{'fraction_non_alphanumeric': 0.03381607530774801, 'fraction_numerical': 0.01217958001448226, 'mean_word_length': 4.38451341235184, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
SOL-KiT -- fully implicit code for kinetic simulation of parallel electron transport in the tokamak Scrape-Off Layer

Here we present a new code for modelling electron kinetics in the tokamak Scrape-Off Layer (SOL). SOL-KiT (Scrape-Off Layer Kinetic Transport) is a fully implicit 1D code with kinetic (or fluid) electrons, fluid (or stationary) ions, and diffusive neutrals. The code is designed for fundamental exploration of non-local physics in the SOL and utilizes an arbitrary degree Legendre polynomial decomposition of the electron distribution function, treating both electron-ion and electron-atom collisions. We present a novel method for ensuring particle and energy conservation in inelastic and superelastic collisions, as well as the first full treatment of the logical boundary condition in the Legendre polynomial formalism. To our knowledge, SOL-KiT is the first fully implicit arbitrary degree harmonic kinetic code, offering a conservative and self-consistent approach to fluid-kinetic comparison with its integrated fluid electron mode. In this paper we give the model equations and their discretizations, as well as showing the results of a number of verification/benchmarking simulations.


Introduction
The heat flow onto the plasma facing components of both present day and future magnetically confined fusion (MCF) devices is of considerable importance [1,2], as it will greatly affect the lifetime of the material. This is true in both steady state operation and during transients (such as ELMs -Edge Localalized Modes). Understanding the heat flux in the Scrape-Off Layer (SOL) is thus of key importance for the design and operation of future fusion devices.
Classic fluid modelling of the parallel (to the magnetic field lines) energy transport in the edge region of MCF devices has relied on the fluid closure of Braginskii [3], or otherwise on various flux limiter approaches [4]. However, it is now well known that there exist discrepancies between experiments and the widely used fluid codes. These discrepencies can be, at least partly, attributed to the effect of non-local transport in the SOL [5]. In this paper, the term "non-local" is used to describe behaviour that stems from strong departure of the electron distribution function from a Maxwellian, in particular due to the fact that electron-ion mean-free paths in situations of interest are comparable to or greater than the temperature gradient scale lengths.
In the following text we focus mainly on aspects of the divertor SOL. The main feature of the divertor configuration is that the location of the primary plasma-surface interaction is relatively far away from the hot core [6] to specifically designed target plates. We distinguish between the "upstream", closer to the core, and the "downstream", where the plasma near the divertor targets is considerably cooler, and the ionization degree can be well below 100%, rendering plasma-neutral interaction important. As such, a large temperature gradient is present, and plasma collisionality (measured with the electron-ion collision mean free path λ ei ) varies greatly along the magnetic field lines in the SOL. Of critical importance is the ratio of the mean free path to the temperature gradient scale length L ∇T = (∇ || T/T ) −1 . Once the ratio λ ei /L ∇T is no longer much less than unity, the classical transport results are no longer valid [4]. However, another important concept in the understanding of energy transport in the SOL is that of the high energy heat-carrying 2 P R E P R I N T electrons (HCE), which become marginally collisionless before the bulk of the distribution [7,8], and will remain so in most SOL situations. In other words, even if the ratio λ ei /L ∇T might still imply the correctness of classical transport coefficients, the HCE could be collisionless, and thus modify the transport by producing non-Maxwellian distribution functions. A further complication in the understanding of the SOL, as mentioned above, is the importance of electronneutral interactions. This is especially true during detachment [9], when the ionization degree drops considerably, and a neutral cloud is formed between the divertor targets and the upstream plasma. As this regime of operation offers better protection to the divertor plate materials, it becomes important to understand the interplay of kinetic/non-local effects already present in the SOL with detachment.

SOL kinetic modelling
In order to properly capture non-local effects in the SOL it is necessary to treat the plasma using a kinetic approach.
Broadly speaking, the two main approaches in the kinetic modelling of the SOL are the Particle-in-Cell (PIC) and finite difference methods solving the Vlasov-Fokker-Planck equation sometimes combined with the Boltzmann collision integral.
A representative PIC code for SOL simulations is BIT1 [10,11], used in a variety of simulation scenarios corresponding to present day machine conditions. The PIC method is highly parallelizable, and naturally accommodates the addition of many different collision types (e.g. tungsten impurities [12,13]). While detailed simulations covering many aspects of SOL transport are possible, two issues make PIC codes complicated to operate.. These are the usually long run times (compared to finite difference methods and fluid codes), and the fact that the number of particles simulated in a PIC code can never approach reality, requiring smoothing techniques [14], and potentially not resolving the high energy tails of distributions with enough accuracy.
Finite difference methods do not suffer from the noise problems of PIC codes, as they solve for the distribution function directly. A number of finite difference codes with different approaches have been utilized in the modeling of the SOL, with a few examples mentioned here. An early example of a completely kinetic code (treating every species kinetically) was the ALLA code [15]. Another code, utilizing a similar method to what is presented in this paper, albeit with an explicit algorithm, is the FPI code [16,17,18], where electrons are treated kinetically while others species are stationary. More recently, the code KIPP [8,19,20], with kinetic electrons, has been coupled with the 2D fluid code SOLPS, providing the latter with kinetically calculated transport coefficients.

Motivation to develop SOL-KiT
Due to the great mass difference between electrons and other species within a hydrogen plasma, electrons mainly suffer pitch-angle scattering collisions when colliding with those heavier particles. Eigenfunctions of such collision operators are spherical harmonics, and an expansion of the electron distribution function in spherical harmonics becomes natural [21]. This approach has been used in the modeling of Scrape-Off-Layer transport to a limited extent [16,17,18], but has been used both in codes dealing with laser-plasma interactions (KALOS [22]/OSHUN [23], 3 P R E P R I N T and IMPACT [24]), as well as electron swarm transport models [25,26]. The expansion has been used to efficiently model both Coulomb and electron-neutral collisions, and has proven itself to be a powerful tool in treating plasmas of various collisionality. With this in mind, a marriage of the Vlasov-Fokker-Planck (VFP) approach in laser plasmas and the Boltzmann approach of electron swarms in neutral gases seems to be potentially a highly applicable model for the Scrape-Off Layer plasma, where collisionality changes, and neutral-plasma interactions carry a great deal of importance.
As is typical with solutions of differential equations, the boundary conditions tend to define the system behaviour and dictate the approach in the numerical solution. In modeling the SOL, it becomes necessary to incorporate the effect of the plasma sheath formed at the boundary, i.e. at the divertor target. While the traditional approach of Procassini et al. [27] has been used in many kinetic codes, when utilizing the spherical harmonic expansion it becomes necessary to formulate the well known boundary condition in terms of the expansion basis. We present this formulation, and give its implementation in SOL-KiT.
The divertor target plate acts as a sink of particles, and as such generates flows towards the target. Since the divertor boundary condition is formulated in the lab frame, it is then necessary to treat the ion flow in the lab frame as well. As a consequence, the electron-ion collision operator must be extended to account for the moving ions. This is a different strategy for incorporating ion motion in the electron VFP equation, than used elsewhere. There the Vlasov terms are transformed instead into the local rest frame of the ions [28]. We present a simplified treatment of this lab frame operator in the next section, together with fluid ion equations. In order to be able to resolve the low velocity cells for higher harmonics (avoiding the CFL condition) it then becomes necessary to treat the electron kinetic equation implicitly.
Another consequence of the divertor boundary conditions is that it also acts as a source of neutral particles. The inclusion of electron-neutral collisions on a nonuniform velocity grid poses particle and energy conservation problems.
We present a method of mapping inelastic collisions on a nonuniform velocity grid which conserves both energy and particles, as well as obeying a numerically consistant detailed balance condition when calculating superelastic collision cross-sections. In order to self-consistently model the interaction of atomic states and the electrons, we include a diffusive-reactive Collisional Radiative model [29] for the evolution of hydrogenic atomic states.
Finally, in order to provide a modular, self-consistent one-to-one comparison of kinetic and fluid modelling, the code also includes fluid equations for the electrons, which can be solved instead of the kinetic model, while making sure that all of the physical processes are the same, and that the atomic data is used consistently between the two models.
To the authors' knowledge, this is the first fully implicit arbitrary Legendre polynomial/Spherical Harmonic code that has an inbuilt sheath boundary condition, inelastic electron-neutral collisions, and a self-consistent fluid mode for clean comparisons. In the following sections the equations of SOL-KiT and their numerical implementation will be presented, starting with the analytical aspects of the model in section 2, before moving on to the model's numerical implementation in section 3. Finally, details of performed benchmarking runs will be given in section 4. We discuss 4 P R E P R I N T the various aspects of the code in section 5.

Physical model
In this section we will introduce the equations being solved in the SOL-KiT model, giving a condensed overview of the physics before expanding on individual operators.
While the code is capable of handling both fixed and periodic boundary conditions, since the most involved cases utilize the Scrape-Off Layer domain, we start with describing it. The domain is 1D, and is assumed to be along a (straightened) field line. The field line is taken as the x-axis, around which the domain is symmetric. The point x = 0 is taken to be the symmetry plane, representing the "upstream" of the SOL. A sketch of the simplified SOL domain is given in Figure 1.  : the x-axis is the principle axis of the system, with x = 0 being the upstream symmetry plane; the right boundary of the system is at the divertor target, which acts as a sink for the plasma, as well as a source of neutrals via the recycling flux

Γ REC
The equations solved by SOL-KiT are the following: • Electron equations -either fluid (density, parallel velocity, temperature) or kinetic • Ion fluid equations -density and parallel velocity (assuming either T i = 0 or T i = T e ) • Diffusive-reactive Collisional Radiative Model for atomic states where E is the electric field (assuming azimuthal symmetry and straightening out the magnetic field we ignore magnetic effects). The RHS contains all of the collision and source operators. Details on collision operators and the Legendre polynomial decomposition of the electron distribution function are given in Section 2.1. In the electron fluid mode, the continuity, momentum, and temperature equations are solved instead.
The continuity equation is given by where n e ,u e , and T e are the electron density, flow velocity, and temperature, respectively. S = S ion + S rec is the particle source including ionization and recombination, and R ei = R T + R u is the classical Braginskii [3] friction where the τ e is the Braginskii collision time [3]. R en is the total friction from all electron-neutral collisions, assuming a slowly (compared to electron thermal speed) drifting Maxwellian distribution for electrons.
The electron temperature equation is where q e = q T + q u is again the classical Braginskii [3] heat flux where κ is taken to be either the Spitzer-Härm result or the Lorentz result [30,31], depending on whether we treat electron-electron momentum transfer collisions. Q = Q ext + Q en is the combination of the external heating energy source (see 2.1.4.), as well as any inelastic collision energy transfer between the electrons and the atoms.
The ion equations are analogous, with a few differences. Firstly, the ion continuity equation can be solved, or quasi-neutrality can be enforced artificially by setting Zn i = n e . Secondly, there is no ion temperature equation in the current version of the code. Instead, ion temperature is set to either zero, or to the electron temperature. The equations where R ie is obtained by requiring total momentum in electron-ion collisions to be conserved, ie. R ie = −R ei . R CX is a simple charge exchange friction term, given by where the sum is over neutral atomic states, and both the ions and neutrals are approximated as cold, with the simplified constant hydrogenic charge exchange cross sections given by approximate low energy values obtained from Janev [32] To calculate the electric field, we solve Ampère-Maxwell's law, which contains only the displacement current where j e is either given as a moment of the electron distribution function, or simply as j e = −en e u e in the electron fluid case.
Finally, the atomic state distribution of the neutrals must be tracked. This is done using a diffusive-reactive Collisional Radiative model (CRM) to obtain the evolution of the neutral state densities n b (where b here denotes the principal quantum number of the state) where the ionization and (de-)excitation rates K, as well as three-body recombination rates α are calculated using moments of the distribution function (see 2.1.3.). The inelastic cross-sections and radiative de-excitation/recombination rates A and β are all taken from Janev [32] and NIST [33]. Radiative de-excitation is included only up to state number 7 P R E P R I N T b = 20, due to lack of available data. Since higher excited states are primarily collisionally coupled in most situations of interest this should not cause significant discrepancies.
The classical 1D diffusion coefficient is simply where v tn is the thermal speed of neutrals, σ el is the elastic collisions cross-section (see electron-neutral elastic collision operator in 2.1.3.), and n 1 is the ground state density. When charge exchange is used, it is the dominant term in the diffusion coefficient, but the elastic collision diffusion is included for cases when charge exchange is turned off.
Since gas temperature and elastic cross-section are free parameters in SOL-KiT, this operator can be tuned. Ideally, however, a self-consistent neutral dynamics model should be implemented, and this is a planned extension of the model.
The boundary condition at the divertor target is the logical boundary condition [27] when electrons are treated kinetically, while the ions are always assumed to reach sound speed (as per the Bohm criterion). Details of the kinetic boundary condition are given below. In both the fluid and kinetic model, the flow into the sheath is ambipolar Γ e = Γ i . When electrons are treated as a fluid the sheath heat transmission coefficient [6] in q sh = γkT e Γ e is taken to be γ e = 2 − 0.5 ln(2π(1 + T i /T e )m e /m i ). Finally, the atomic neutrals are recycled with a recyling flux where R ≤ 1 is the recycling coefficient. This is simply imposed by setting D 1 ∂n 1 /∂x in (13) at the boundary to Γ REC , whereas it would otherwise be zero.

Electron kinetic equation in Legendre formalism
Spherical harmonics are an orthonormal basis set in the solid angle space (θ, ϕ) ∈ [0, π) × [0, 2π), and can be written as a suitably normalized product of associated Legendre polynomials P m l (cos θ) and the complex phase e imϕ . The spherical harmonic convention used in SOL-KiT is same as the ones in KALOS [22]/OSHUN [23]. As such, the traditional Cartesian coordinate system for velocity (v x , v y , v z ) will be rotated so that the angles θ and ϕ are defined as part of a spherical coordinate system: If then we decompose the distribution function f (v, θ, φ) in spherical harmonics we get: where the complex expansion coefficients f m l satisfy ( f m l ) * = f −m l . This allows writing the kinetic equation for electrons as a set of equations for the amplitudes f m l , whose physical significance becomes obvious when moments of 8 P R E P R I N T the distribution function are expressed using the expansion. If φ is a scalar function of v then and similarly for higher order tensors [21,23]. As can be seen from this, not only does the spherical harmonic expansion have useful properties in relation to collisions, it also provides a physically meaningful decomposition for evaluating transport quantities.
Since the current model is 1D, and azimuthal symmetry around the x-axis is assumed, we do not treat magnetic field effects, and the spherical harmonic decomposition reduces to a Legendre polynomial decomposition (m = 0 always). Accordingly, in the rest of this paper, harmonics will be labeled only by their l-number, i.e. f l (v).
Following the Legendre decomposition of the distribution function as outlined above, the 1D kinetic equation can be written as a set of equations where now all of the operators are moved to the RHS and are functions of l in addition to whatever arguments they naturally have. These will be examined in detail below.

Vlasov terms
The terms on the LHS of equation (1) are usually referred to as the Vlasov terms. The two Vlasov terms in equation (19) are the spatial advection term A l (corresponding to second LHS term in (1)), and the velocity space advection term due to the electric field in the x-direction E l (corresponding to third LHS term in (1)).
Firstly, the spatial advection term (advection in the x-direction), for a given harmonic l is Spatial advection couples harmonics with different l numbers. The physical significance of this coupling is most easily seen in the coupling between f 0 and f 1 . The moments of f 0 are the density and total energy, while f 1 is associated with flows (of particles and energy). Thus gradients in f 0 (density and temperature) drive advection of f 1 (flows), and vice-versa.

P R E P R I N T
The velocity space advection term due to the electric field couples harmonics as well, albeit through velocity space gradients in f l . As only the x component of the electric field is treated, in the remainder of the paper it will simply be written as E for brevity. The velocity space advection operator is given by [22] where As is evident from equation (21), the electric field couples harmonics through the G l and H l functions, which contain velocity space gradients of the coupled harmonics.
Thus, Vlasov terms provide coupling of harmonics through either spatial gradients or the electric field.

Coulomb collision terms
Let us consider the effect of Coulomb collisions on the distribution function f of particles of mass m and charge q = ze colliding with particles of mass M = µm and charge Q = Ze with distribution F. We follow the formalism of Shakorfsky et al. [21], starting from the Rosenbluth coefficient formulation of the Fokker-Planck collision operator for Coulomb collisions where ∇ = ∂/∂ v and Γ zZ = (zZe 2 ) 2 ln Λ/(4π(m 0 ) 2 ), with ln Λ denoting the Coulomb logarithm. The Rosenbluth drag and diffusion coefficients are respectively H and G. We separate the distribution functions into their isotropic and anisotropic components F = F 0 + F a , f = f 0 + f a . The key assumption going forward is that the anisotropic component is small compared to the isotropic one, so that it becomes possible to linearize equation (24). Expanding the distribution function and the Rosenbluth coefficients in harmonics and using the integrals[21] one can derive the expressions for the harmonic components of the Fokker-Planck collision integral for all l. For l = 0 this is For l > 0 the following is obtained [23] P R E P R where the C coefficients are functions of l .
For electron-electron collisions µ = 1. The effect of e-e collisions on the isotropic part of the distribution function is given by where the drag and diffusion coefficients are defined as Note that D is not in the form one would expect from equation (26). Instead, it is written in an analytically equivalent form (see section 3.5.). The electron-electron collision operator for l = 0 is important for the proper relaxation of the electron distribution function to a Maxwellian.
For higher harmonics, from (27) we get As the ions are either assumed cold or with the same temperature as the electrons, in the current version of SOL- KiT the effect of electron-ion collisions on f 0 is not included.
For higher harmonics we distinguish two cases of electron-ion collisions. The first is the classical stationary ion case, where F 0 = n i δ(v)/(4πv 2 ). It can easily be shown that equation (27) reduces to the following eigenfunction form Here it can be seen that this is purely angular scattering, which dampens higher harmonics, helping us truncate the expansion.
However, when ions are not stationary, but are moving at some velocity much smaller than the electron thermal velocity (the case we expect in the SOL), it is necessary to modify the electron-ion collision operator. This is done by first getting rid of all the terms proportional to the inverse mass ratio in equation (27). To calculate the I and J integrals for a cold ion stream we let the ion distribution function be Recasting the Dirac delta into spherical coordinates assuming azimuthal symmetry and expanding in Legendre polynomials gives us Substituting these harmonics into the equations for the I and J integrals gives us where Θ denotes the Heaviside step function. It is now trivial to see that all but the I 0 integrals vanish when u i = 0, recovering the stationary ion collision integral. It can also be easily shown that for small enough u i , the collision integral for f 1 reduces to [21] 1 If f 0 is taken to be Maxwellian, this gives a stationary solution to the f 1 collision operator to be that f 1 which yields a slowly drifting Maxwellian with drift velocity u i .
Taking a closer look at the obtained collision integral, we see that electrons slower than the ions are being carried around by them, as one would expect. This modifies all harmonics for small v, so we lose the convenient property of clean truncation of the distribution function by electron-ion collisions. However, in realistic simulation cases, this happens only for a small handful of velocity space cells, while higher harmonics are normally dampened in the rest of velocity space. This is to be expected, as our electrons see an ion Dirac delta and are collisionally driven towards it.

Boltzmann collision terms
We start from the general form of the Boltzmann collision integral for the effect of collisions on the distribution function of species s colliding with species s where primed velocities denote values before a collision, and σ is the appropriate differential cross-section. The following results are all derived under the assumption of a small mass ratio and stationary (slow compared to the electrons) neutral particles (atoms) [21,25,26]. Furthermore, all differential cross-sections are assumed to be azimuthally symmetric, i.e. are only a function of energy/velocity of the impacting particle (electron), and the deflection angle (here χ). This just means that the cross-sections do not depend on the orientation of the neutral particle.
Elastic electron-neutral collisions. The first term in the small mass expansion of the electron-neutral elastic collision integral cancels for l = 0. It is therefore necessary to include higher order effects of collisional energy transfer. This way, one can allow for long term relaxation of the electron distribution function to a Maxwellian with temperature equal to the gas temperature T g . If M is the neutral mass, and n b the gas density, the collision integral takes the form As is the case with electron-ion Coulomb collisions, due to a great mass difference, the l > 0 integral is considerably simplified, and can be written as where P l is simply the l-th Legendre polynomial.
While implemented and tested in the current version of SOL-KiT, these two processes are rarely used, because of the lack of proper elastic collision cross-section data. Currently the cross-section for elastic collisions of an electron with a hydrogen atom in a given state is obtained using the classical expression for orbit size. Namely, this gives for the integral elastic collision cross-section of electrons with hydrogen atoms in the b-th state where a 0 is the Bohr radius.
Inelastic electron-neutral collisions. We start with inelastic collisions where the total number of particles of each species is conserved (e.g. excitation). A standard procedure exists [21,25] for an inelastic collision for which the pre-collision and post-collision velocities are related as where is the inelastic energy loss.
Here σ T OT = dΩσ(χ, v) is the integral cross section, while Using equation (43) for l = 0 one can easilly show that particle number is conserved.
On the other hand, collisional processes such as ionization do not conserve the total number of particles. The main difficulty in treating ionization is the fact that it is, ultimately, a 3-body process, and ideally one would like to know the triply differential cross section for such a process. However, such an approach is not only complicated, but cross-section data (to the author's knowledge) are not systematically available. Because of this, the approach taken here will be the simplest possible [15], where all electrons produced in ionization are put in the lowest velocity cell, while the original electron experiences a standard energy loss/deflection as in the case of excitation. The collisional operator for ionization takes the following form where δ f l δt ex (σ ion b ) is equation 43, but with σ ex replaced with σ ion , and with the collisional ionization rate coefficient defined (unconventionally) as P R E P R I N T Particle sources are then computed using S ion = b K ion b n b , and similarly for recombination, while the inelastic collision contribution to Q is similarly calculated by taking the product of each transition rate and the associated transition energy.
Of course, one operator of the above kinds exists for each possible process of the given kind, i.e. one operator for each excitation process, and one for each ionization process, taking in the appropriate neutral state densities and cross-sections for the given process. The rate coefficients are the same ones used in equation (13).
Inverse processes (deexcitation and 3-body recombination) are treated using the principle of detailed balance [34,35] to obtain cross-sections. For deexcitation, this gives where i and j are atomic states ( j < i), and g i and g j their statistical weights. For hydrogen these are simply g n = 2n 2 .
Equation (42) defines the velocities, but we use a negative .
For 3-body recombination, using the statistical weights of a free electron gas we get for the cross-section where h is the Planck constant, n e and T e are the electron density and temperature, respectively, and g + 1 is the ion ground state statistical weight (for hydrogen simply g + 1 = 1). To calculate the electron fluid mode R en in equation (3) we use terms of the form , and similarly for other neutral processes.

Electron heating operator
The implemented diffusive heating operator has the form where Θ(L h − x) is the step function designating the heating region. It is easy to check that this operator conserves particle number if ∂ f /∂v = 0 on the system boundaries. If we assume a spatially uniform heating, it is easy to show that where W h (t) is the heat flux entering the SOL over length L h . This is related to the fluid model heating Q ext via

. Particle source operator
In order to treat upstream density perturbations, the following electron and ion particle sources (in kinetic mode only) are implemented where F R is the source rate coefficient with Γ in being the effective upstream flux. The particles are injected over a length of L s and with temperature T source , which can be the background temperature.
The ion particle source is simply

Divertor target boundary condition with Legendre polynomials
The boundary condition at the divertor target is calculated using the standard logical boundary condition [27], setting the ion and electron fluxes to be equal at the sheath entrance. The logical boundary condition assumes that all electrons with a parallel velocity above some v c moving towards the target are lost, while all others are reflected.
This translates to having a sharp cut-off in the electron distribution function. The challenge when formulating this condition in a Legendre polynomial formalism is the extreme anisotropy that results from it, which would require a high number of harmonics to resolve to a satisfactory level. Fortunately, this number is usually not prohibitively high (see 4.4.). The harmonic content of the "cut-off" distribution f cl can be written as a linear combination of known For details on the calculation of the transformation matrix P ll see Appendix A. Knowing the form of the distribution function, one can solve the ambipolarity condition where n i,sh is the extrapolated density at the sheath boundary (see 3.7. below), and u i,sh is the ion velocity at the boundary, given by the Bohm condition where T e is the electron temperature in the last simulation cell, and T i is the ion temperature. Solving the ambipolarity condition gives the value of v c , and with it the value of the sheath potential drop ∆Φ = m e v 2 c /(2e). The electron distribution harmonics with the correct cut-off velocity can then be used as the dynamically updated boundary condition.

Numerical Methods
In this section we present numerical details of the SOL-KiT algorithm, starting with the definitions of the normalization scheme and the grids used. An overview of discretization schemes for the various operators in the code follows.

Normalization
The temperature is normalized to some reference value (in eV), while the reference density is assumed to be given in m −3 . These normalization constants will be refered to as T 0 and n 0 , respectively. The velocity is normalized , and the length to the thermal electron-ion collision mean free path x 0 = v th t 0 . Here T 0 [J] denotes the normalization temperature converted from eV to Joules, with and where ln Λ ei (T 0 , n 0 ) is the Coulomb logarithm for electron-ion collisions calculated for the normalization temperature and density (taken from [36]). All normalized quantities arẽ where σ 0 = a 2 0 π (where a 0 is the Bohr radius) . In the following sections the normalized quantities will be written without the tilde in order to lighten the notation.

Grids
The velocity grid is a uniform or geometric grid of cells with (starting) width ∆v 1 and width multiplier c v , with N v cell centres distributed as while the spatial grid is staggered, i.e. consists of cell centres and boundaries. N c is the number of cells (cell centres), while N x is used to denote the total number of spatial points (cells and boundaries), which depends on the boundary conditions of the grid (while N c is an input parameter).
In the following text, spatial points with an odd index (x 1 , x 3 , etc.) will denote cell centres and those with an even index will denote cell boundaries. The values of these points are determined by where m = k if k is odd, or m = k − 1 if k is even. ∆x c m denotes the cell width of the cell whose centre is at m. For a uniform grid, this is constant, while for a "logarithmic" grid it is an exponential function that starts at a prescribed width for the first cell dx, and drops to a prescribed width of the last cell ∆x L (with a fixed number of cells N c ).
Variables are positioned on the staggered grid in the following manner: • In cell centres: Variables not evolved in cell centres are linearly interpolated from neighbouring cell boundaries, and vice versa.

Timesteps, nonlinear iteration, and vectorization of variables
SOL-KiT uses either uniform timesteps of prescribed length ∆t, or rescales ∆t with min(T e (x) 3/2 /n T OT (x)) (where n T OT is the total density of heavy particles). This is a conservative estimate for the shortest Coulomb collision time. In the following text the timestep will be assumed uniform, but generalization to adaptive timesteps is straightforward.
Timestepping is done using a first order implicit backwards Euler method. Let F be the vector containing all the evolved quantities, and M(F) the evolution matrix. Then Within a given timestep, we use fixed point iteration to solve the non-linear system (57), with M(F i * ) being evaluated with the solution at the previous iteration. For the first iteration,F i * = F i . (57) is iterated until convergence is established is established within a prescribed tolerance. Implicit variables (those that appear at time i + 1 on RHS of equation (56)) in the most important nonlinear terms are given in Tables 1-3. b If using kinetic model this is replaced by 4π v 2 f 0 dv.  Term The structure of the variable vector (with all variables present) is the following: where F loc (x k ) is the (spatially) local subvector for quantites at point x k given as where N n is the total number of neutral states tracked, l max is the highest resolved harmonic, and f l (x k ) is the l-th harmonics subvector . . .
Note that when running in kinetic mode, the vector does not contain electron fluid quantities, while when the code is running with fluid electrons the distribution function harmonics are not evolved, but are updated after each timestep to be a slowly drifting Maxwellian with current temperature, density, and electron fluid velocity.
From knowing the input parameters N c (and its derived parameters 2N c − 1 ≤ N x ≤ 2(N c + 2) − 1), l max , N n , and N v we can calculate the total vector length (for a kinetic run) as For a representative system with l max = 5, N v = 80, N n = 30, and N x = 128 this gives N total = 55424.
To solve the above matrix system, we use the MPI and PETSc libraries [37,38,39]. Domain decomposition is done in the spatial dimension with as close to even distibution of grid points between processors as possible.

20
P R E P R I N T

Velocity and spatial derivative discretization
The velocity space derivatives appearing in the various kinetic operators are all implemented using a central difference scheme: where F is any velocity space function.
Spatial derivatives are mostly discretized using an analogous central difference scheme to the above velocity space one. The only exceptions are the advection terms in the momentum and temperature equations (eq. (3), (6), (10)), where an upwind scheme can be used instead of the central difference. The upwind scheme used is simply and if u(x k ) < 0 where if F is not evolved on x k±1 we use use x k±2 instead (this is the case for temperature advection in (6)).
Finally, the diffusion type derivatives are given by where the velocity grid boundaries v n+1/2 are given as v n+1/2 = n l=1 ∆v l . The only simple derivatives that do not obey the above are in the velocity advection terms, where in equations (22) and (23) we use the conservative form where F(x k , v n+1/2 ) is obtained through linear interpolation.
The following velocity space boundary conditions are assumed for the distribution function harmonics i.e. f 0 at v = 0 is quadratically extrapolated (this is used whenever f 0 at v = 0 is required) and higher harmonics are set to 0. At the velocity space boundary after v N v all f l and their derivatives are assumed to be zero.

Coulomb collision operators
Here, Coulomb collision operator discretization will briefly be presented to supplement the material in previous sections. The method used for the electron-electron collisions for l = 0 is the Chang-Cooper-Langdon [40,41,42] method. The implementation in SOL-KiT is very similar to that in IMPACT [24], and as such we will leave out some of the details. We write the collision operator as a divergence of fluxes F where A 0 ee = Γ 0 ee n 0 t 0 /v 3 t = 1/(Z ln Λ ei (T 0 , n 0 )), and the flux is given by f i+1 0 (x k , v n+1/2 ) is then calculated using a special weighted interpolation, which ensures relaxation to a Maxwellian The friction coefficient is where v * n+1/2 = (v n + v n+1 )/2. Boundary conditions for the C and D coefficients are taken in such a way to conserve particle density (see [24]). The resulting submatrix from this operator is tridiagonal. Thus we have an iterative method which conserves particles and energy (up to nonlinear iteration tolerance), and which relaxes the distribution to a Maxwellian in the absence of other driving forces.
The electron-electron collision terms for higher l (equation (31) have the same normalization constant as that for l = 0. Here we use a discretization method similar to that in OSHUN [23], and will hence again, for the sake of brevity go over only the most important elements. The first three terms in (31) produce a tridiagonal matrix when discretized, while the remaining terms produce an upper and lower triangular submatrix due to the the I and J integrals, which are discretized as For l = 1 the discretized e-e collision operator does not numerically conserve momentum, unfortunately, as the discertized form loses the partial integration properties that would analytically conserve momentum. However, the numerically lost momentum is transfered to the ions, and total momentum is thus conserved.

P R E P R I N T
The electron-ion operator for stationary cold ions is trivial, and is discretized straightforwardly, while the moving ion operator is discretized similarly to the e-e operator for higher l, with the exception of having mainly a tridiagonal component (with terms containing F 0 ), and a part (terms containing F l ) where ion density is the implicit variable (if working on a cell boundary, it is implicitly interpolated) . For the second part we need the previously presented I(F l ) and J(F l ) integrals for cold ions ( (35) and (36)). For a given lagged ion velocity u i * i (x k ) and density n i+1 i (x k ) these integrals in discrete form are

Electron-neutral collision numerics
In this section we briefly go over the basic properties of the elastic electron-neutral collision operator before moving on to the important conservative discretization of inelastic electron-neutral collisions.
Elastic collisions. The discretization of elastic collisions borrows greatly from the Coulomb collision operators.

Equation (39) is discretized in a way similar to the Chang-Cooper-Langdon scheme. However, since the integral
is linear in f 0 , there is no need for special interpolation ( f 0 and σ are simply linearly interpolated), and just the flux formalism was used, with the C and D coefficients (here left unnormalized) being where Since no differential cross section data is available, this is simply set to the constant cross-section (41).
The above discretization produces an operator that conserves particle number and relaxes the electron f 0 to a Maxwellian with temperature T g (within finite grid effects). For higher l, equation (40) is discretized in a straightforward way.

24
P R E P R I N T

Conservative discretization scheme for inelastic collisions
In order to streamline the arguments in the next sections, let us introduce the following notation. If we label inelastic processes with transition energy as Π, we can label the corresponding superelastic (inverse) processes as Π −1 (with transition energy − ). Electron energy and particle conservation are governed by the l = 0 harmonic, which we will denote simply as f in the following derivation. Then the value of f in the n-th velocity cell centre can be written as f n , while the total cross-section notation for any given process will be labeled as σ n .
Equation (43) contains two terms, a loss/emission term and a gain/absorption term. While the first term is defined on the used velocity grid, the second term requires evaluation of the distribution function on (most likely) non-existant points α(v)v (where α = (1 ± 2 /mv 2 ) 1/2 ). A straightforward interpolation here fails to produce a discretization that conserves particles and energy. In order to develop such a scheme we start by writing the collision operator in the emission/absorption form where E n = n neut v n f n σ n , and where W nm are weights determining the contribution of the emission from cell m to the absorption in cell n. For particle conservation, we want the number of particles emitted by a single cell m -4πE m v 2 m ∆v m to be absorbed exactly by some other cells n -4πW nm E m v 2 n ∆v n . After tidying up, the resulting particle conservation condition is v 2 m ∆v m = n W nm v 2 n ∆v n .
For energy conservation, we note that all energy being emitted via collisions from one cell needs to either be lost to the energy of internal atomic states or be absorbed by some cells n. In a similar way to the above, one can cancel distribution functions and cross-sections to obtain the energy conservation condition where n l and n u denote the densities of the lower and upper excited states, respectively. Using a Maxwellian for f , relating n u and n l via the Boltzmann distribution, and utilizing equation (84) we obtain Note that this depends on the temperature and the excitation weights, contrary to the analytical version from equation (46). This is due to the discrete nature of the grid, where energy differences between absorbing and emitting cells do not have to be equal to the analytical value .

Two-absorber mapping
Note that in the above conservation condition we haven't specified the summation range for n. The simplest way to do this is the following. For each emitter m we choose exactly two (consecutive) absorber cells, n 1 and n 2 such that v 2 m ∓ lies between v n 1 and v n 2 . We will refer to these two points as the ideal absorber pair. This way we do not need any further instructions on partitioning the emitted particles and energy, and can proceed to calculate weights. Since for any given m there are only two absorbing cells, we can denote the weights simply as W m 1 and W m 2 , and solving the conservation conditions we get .
(89) Figure 2 shows an illustration of this mapping for an excitation collision with both particle and energy transfer obeying conservation conditions.
To establish that this scheme can reduce to the analytical result, we note two properties we expect in the analytic limit. Firstly n 1 → n 2 = n, i.e. each emitting point has one and only one absorbing point. This can be done by letting one of the weights tend to 0. For the sake of this illustrative argument, let that be W m 1 , from which we see that the second fraction in (89) tends to unity. Then, we note that the first fraction will tend to α(v n ) 2 ∆v m /∆v n (since v m = α(v n )v n in analytical limit). Taking the differential of (42) points us toward a grid that satisfies Noting that we have chosen to write the gain term as a sum of weighted loss terms, another α factor can be extracted from the single emitter cell velocity in the absorption term, and the analytical α 2 form from (43) is recovered.
At first glance, this mapping looks good, we have found pairs of absorbers and weighted their absorption terms to conserve particles and energy, but after a closer look at eq. (87) a potential problem reveals itself. Suppose that for some n and every m the weights W nm = 0 in an excitation process. In this case σ deex n = 0 for the corresponding deexcitation process. In other words, if cell n does not absorb in process Π, it will not emit in Π −1 , and gaps are left on our velocity grids where cells that should emit do not. This is not physical, and the mapping needs to be refined.

Two-absorber mapping with pair partition
In the previous section we have associated every cell (emitter) m with an ideal absorber pair (n 1 (m), n 2 (m)). We define the absorber list A Π m as a list of absorber pairs of point m in process Π. The two-absorber mapping implies A Π m = {(n 1 (m), n 2 (m))}, where the only pair is the ideal absorber pair. As noted above, this produces gaps in the the velocity grid where cells do not emit in the inverse process Π −1 .
In order to fill out the aforementioned gaps, we must potentially include absorbers other than those in the ideal pair. This can be done by looking at all cells n that aren't in the ideal absorber pair, but whose v 2 n ± (note the change in sign!) falls into cell m. We then refer to m as the ideal emitter for cells n.
We then move to generate absorber pairs that would include the non-ideal absorbers n, while satisfying the original constraint, namely that v 2 m ∓ is between the points of each new pair. In SOL-KiT this is done by always pairing P R E P R I N T non-ideal absorbers with one of the ideal absorbers. If we denote P m ≥ 1 as the number of (both ideal and non-ideal) absorber pairs of cell m, the absorber list becomes A Π m = (n p 1 (m), n p 2 (m)), p = 1, ..., P m , where we label the first and second cell in pair p as n p 1 (m) and n p 2 (m), respectively. This allows us to use the previous solution, but we require a way to partition the emitted energy/particles among pairs. We do this by defining a total energy width of the entire absorber list and normalizing the energy widths of each cell within the list to the total energy width of the list: 2v n ∆v n , β n = 2v n ∆v n /β m T OT .
Then we define δ n as the number of pairs in A Π m which contain point n. We can then define to be the fraction of the emission given to pair p. It is trivial to check that the sum of all γ p for a given absorber list is equal to unity, as one would expect. An example with an absorber list with three distinct cells is given in Figure 3.  Then we can go through the same process as the one for the two-absorber case, except the final results for W m 1 and W m 2 will now be functions of the pair p for which they have been calculated, and will have an extra γ p factor multiplying the previous simple results (see below). To then calculate the total absorption term for a given cell n, it where we scale the harmonics by the ratio of the electron densities at spatial cells N x and N x − 1. The choice to perform this sort of extrapolation, and not a linear one, comes from the danger of large density gradients (which are often present at the divertor target) to produce negative extrapolated densities when extrapolated linearly.
Knowing the forward going distribution allows us to calculate the cut-off distribution using equation (53), while equation (54) imposes v c . However, due to the discretization of velocity space, this will not be one of the resolved cell centres, and we are required to interpolate in order to calculate the precise cut-off. A diagram of the interpolation region is given in Figure 4.
The interpolation replaces cell K containing the cut-off and its preceding cell K −1 with two new cells with centres The electron flux to the boundary (given in equation (54)

Benchmarking SOL-KiT operators
A number of verification tests have been performed using SOL-KiT, aimed at checking the properties of various implemented operators. The test details are presented in the following sections, while an overview of the tests and their scope is given in Table 4. In all tests we consider hydrogenic species, ie. Z = 1.
P R E P R I N T

Heat conduction tests
In order to test a number of operators we performed both local tests with different scenarios, as well as non-local perturbation tests.
Local heat conduction tests. Two types of tests were used to verify the local limit of heat conduction in SOL-KiT.
The first type is a small perturbation test on a periodic grid, with the initial temperature given by where L is the total length of the system, given by L = N c ∆x. The density was set to n 0 = 10 19 m −3 , and the rest of the simulation parameters are with ∆v,∆t, and the total number of timesteps being able to vary between runs. The total time simulated was always set to 30 t 0 , i.e. to 30 electron-ion collision times. In all runs the full set of Coulomb collision operators was active, and the highest harmonic number was kept at l max = 1. The calculated conduction heat flux was compared to the Spitzer-Härm (SH) heat flux q S H (equation (7)), with the heat conductivity for Z = 1 in SOL-KiT normalized units being κ = 0.6 √ π.
The ratio of the calculated heat flux and the reference SH flux at the end of each run would be averaged along the simulation domain, and these are the results presented below. Three pure kinetic runs (with initial f 1 set to 0, and no fluid ion operators active) were performed. The reference run (Run 1) with ∆t = 0.1 t 0 , ∆v = 0.1 v th , N t = 300 produces an average heat flux of q = (0.988296 ± 1 × 10 −6 )q S H , i.e. reproduces the reference results with less than 1.2% error. Run 2 used a smaller timestep ∆t = 0.05 t 0 , N t = 600, but the ratio obtained was the same as in Run 1.
Run 3, on the other hand, tested the velocity resolution dependence by setting ∆v = 0.05 v th , and the obtained heat flux was q = (0.997488 ± 1 × 10 −6 )q S H , reducing the relative error below 0.3%. At this point we believe the relative error is smaller than the precision of the reference SH flux value, and should thus be taken with a grain of salt.
The next set of small perturbation tests aimed to compare the results of kinetic simulations to those performed using the fluid mode of SOL-KiT. For this purpose, since the fluid model assumes the reference SH heat flux described above from the very start of the simulation, we initialized f 1 in the kinetic simulation (Run 4) to a local solution which would give the same heat flux as the one in the fluid run (Run 5). Other than this, the parameters of Run 4 were identical to those of Run 1, and the final heat flux ratio was the same as the one there. Run 5 was performed with the same parameters as Run 4, but instead of solving the kinetic equation for the electrons, fluid equations for both electrons and ions were solved. The total changes in the temperature value at its maximum for both runs were compared. It was P R E P R I N T found that the relative diference in the changes was less than 0.3%, showing good agreement between the kinetic and fluid model in the small perturbation test.
The final heat conduction test run (Run 6) was done for the fluid model, where the plasma was initialized on a logarithmic spatial grid (∆x = 8.5 x 0 and ∆x L = 0.5 x 0 ) with the Two-Point Model [6] profile where n u = 0.368 × 10 19 m −3 and T u = 18eV are the upstream density and temperature, while T d = 5eV is the downstream temperature. Boundary conditions were fixed. We expect that this profile should not evolve, and after 30000 collision times (for a reference plasma of n 0 = 2.5 × 10 19 m −3 and T 0 = 10eV), the largest deviation from the initial profile was 0.12 %, which is due to linear interpolation close to the downstream boundary, with most points having less than 0.01% deviation.
Non-local heat conduction tests. In order to test how SOL-KiT handles a more non-local regime, we have performed a set of runs (Set 1) with a setup similar to Run 1 and related runs, featuring a sinusoidal temperature perturbation.
For completeness, we give the common run parameters for Set 1 The various runs in the set were performed with different system lengths and with varying number of resolved harmonics l max = 1, 3, 5, 7. In Figure 5 we plot the κ/κ (B) = q/q S H for each of the runs, comparing it with values obtained with the code KIPP [8,19,20] and reported by Brodrick et al. [43]. The ratio q/q S H was obtained by running each run in the set until the ratio equilibriated. The different values in Figure 5 are plotted as a function of kλ (B) ei , where λ (B) ei = 3(π/2) 1/2 x 0 /4 is the Braginskii electron-ion collisional mean freepath, as used by Brodrick et al., and k = 2π/L is the perturbation wavenumber. KIPP results were fitted using the following function based on equation (25) in [43] κ where the values obtained for the parameters a and b are 51.409789 and 4.4682314, respectively. We show that the increase in number of harmonics leads to an increase in agreement between SOL-KiT and KIPP. Already at l = 5 SOL-KiT results appear to be very close to the fit values, with the increase to l = 7 having a negligible effect on the result. We also note that the diffusive approximation (l = 1) appears to break down around kλ (B) ei = 0.075, while l = 3 seems to hold up further into the non-local regime.

Collisional-Radiative model and inelastic collision tests
The discretization scheme for inelastic collisions presented in this paper is designed with three properties in mindparticle and energy conservation, as well as numerically consistent detailed balance. In this section we present several tests aimed at confirming the desired properties.
Conservation property tests. We start with two runs differing only in the utilized velocity space grid. The common run parameters are n 0 = 10 19 m −3 , n e = n 0 , n 1 = 0.1n 0 , N n = 30, where N n is the total number of resolved hydrogen states. Both runs included only the inelastic electron-atom processes, and were performed in quasi-0D (low number of spatial cells). The first of the two runs (Run 7) used a uniform grid with ∆v = 0.05v th , while the second (Run 8) used a geometric grid with ∆v 1 = 0.01 and c v = 1.025. Figure 6 shows the average heavy particle density (n i + n b ) in the two runs. As can be seen the total error is of the order of 10 −14 , which is consistent with round-off errors. Similar results are obtained for the relative deviation of the total energy density E T OT = 3n e kT e /2 + n e ion + n b b from initial value, where b is the 1 → b transition energy.
This result is shown in Figure 7, where we see that the geometric grid (Run 8) is performing somewhat better in the later part of the simulation, likely due to a finer velocity grid near the origin.
Detailed balance test. In order to test the numerical detailed balance condition for the inverse processes we treat the electrons as a fluid, thus forcing the distribution function to a Maxwellian. This simulation (Run 9) was performed 1e−14+1.1 (b) Geometric grid -Run 8 Figure 6: Evolution of total density in the two inelastic collision discretization test runs.
with the same grid, normalization, and initial conditions as the uniform grid run discussed above. The only difference was the timestep parameters being N t = 1000 and ∆t = 100t 0 . By the end of the simulation the atomic states settled close to a Saha-Boltzmann equilibrium, as would be expected from the detailed balance condition. The ionization degree relative error computed against the analytical solution for hydrogen was δX = 1.134 × 10 −8 , while the total density and energy errors at the end of the simulation were ≈ 6 × 10 −14 and ≈ 5 × 10 −14 , respectively.
Note that in all three runs discussed above we disabled all radiative processes, as those would mask any energy conservation processes in the first two runs, and would not allow for equilibrum in the detailed balance test.
Integrated test. Finally, an integrated test (Run 10) was performed in order to test the full interplay of electon and neutral processes. We have attempted to replicate as closely as possible the simulation performed by Allais et al. [18] using the code FPI. Parameters used in this run were T 0 = 10eV, n 0 = 2.5 × 10 19 m −3 , N n = 30, n T OT = n 0 , amounting to a domain of length L = 20.31m and total simulation time t T OT = 49.96µs. At x = 0 the temperature is initialized to 25 eV, and the plasma is 100% ionized. From x = 2.5m to x = 4.06m the temperature drops exponentially to 1eV, while the ionization degree drops to 10%. All inelastic collisions were enabled, as well as radiative processes, while ions and neutrals were left stationary as in the original paper. Boundary conditions were fixed to initial values. Figure 8 shows the evolution of electron temperature and density corresponding to Fig. 1. and Fig. 2. in [18]. Qualitative behaviour is recovered, while discrepancies of less than 10% are most likely caused by potential differences in the intialization, as well as likely different spatial and velocity space resolutions (not reported in [18]). Another potential cause of discrepency is the use different databases in SOL-KiT compared to FPI.

Ion flow tests
Acoustic wave test. To test the ion advection, we performed the following isothermal ion acoustic wave test (Run 11).
Parameters of the simulation were  N v = 120, ∆v = 0.1v th , N c = 128, dx = 0.0078125x 0 , with the density and ion velocity initialized on a periodic grid of length L = x 0 as Electron-electron collisions for l = 0 and electron-ion collisions were turned on in the simulation. The density and ion velocity are presented at several points during the evolution of the acoustic wave in Figure 9.

Divertor boundary condition tests
Analytic limit convergence test. The first condition a cut-off logical boundary condition would need to satisfy is the analytical Maxwellian limit for the sheath properties. In order to test this aspect of the operator, the cut-off procedure and calculation of the sheath heat transmission factor and potential drop were performed without evolving the initially Maxwellian distribution. A set of single step ∆t = 0 simulations (Set 3) with different velocity space resolutions (constant v max = 6v th ) was used to evaluate convergence of the cut-off calculation without evolving the distribution.
The results of this test are shown in Figure 11, where we see that both the sheath heat transmission coefficient and the potential are within a few percent of the analytical values [6] ( γ e = 2 − 0.5 ln(2π(1 + T i /T e )m e /m i ) and ∆Φ = γ e − 2) even for the coarsest grid used.  Figure 12. Simulations were performed up to l = 25, and the relative changes in the sheath potential drop and heat transmission coefficient were tracked. Very quickly the change drops below 1%, before nonmonotonic behaviour is observed around l = 9. While convergence appears slow, this is to be expected with such a sharply cut-off distribution being expanded in Legendre polynomials. Fortunately, even in this highly stressed scenario a relatively low number of harmonics l ≤ 5 seems to be enough to capture the physically important behaviour.

Discussion
In The novel discretization scheme for the Boltzman collision integral for inelastic collisions on a velocity grid has been shown to conserve particles and energy up to round-off error. This is to be expected as with N n = 30 hydrogen states, the number of transitions treated, not including radiative transitions, is 930. Given the stiff nature of the Collisional-Radiative model matrix, as well as the corresponding terms in the electron kinetic equations and the large number of transitions, we expected relative errors of order ≈ 10 −13 . Runs 7-8 in the previous section agree with this estimate, with the geometric grid notably having better conservation properties. This is likely due to the fact that the high excited state transitions are better resolved with a finer grid near the origin. Detailed balance has also been shown to be observed within a high accuracy, owing to the numerically adapted detailed balance cross-sections derived in Section 3.6.1.
At the sheath boundary, we implement the logical boundary condition [27]. While this is a standard boundary condition in SOL kinetic codes, this is the first time it has been implemented in a harmonic decomposition code. A reasonable concern in this case would be that the sheath properties would be hard to pinpoint with accuracy, given that the sharp cut-off is poorly approximated with Legendre polynomials. While we do notice the effects of the Gibbs phenomenon in the reconstructed distribution function when high l low collisionality systems are treated, we show in Set 4 that a relatively low number of harmonics is sufficient to provide an estimate of the relevant sheath properties. Furthermore, the Maxwellian analytic limit of the sheath properties requires only resolving up to l = 1. This further supports the notion that, while important in general, high harmonic approximations of the cut-off are not necessary for the transport calculations at the sheath boundary.
SOL-KiT is a fully implicit code, and is thus not limited by the CFL condition. It is, therefore, able to resolve the lowest velocity cells for higher harmonics, unlike explicit codes [23]. This is one of the key features allowing the implementation of both lab frame ion flow, as well as inelastic electron-neutral collisions, with the former being the main beneficiary of the implicit scheme. The reason for this is that the effect of flowing ions in the lab frame is visible primarily in the low velocity cells, where electrons are highly collisional and are dragged with the ions. A scheme that cannot resolve those cells must resort to calculations in the ion centre of mass frame [44,28]. Given the complicated sheath boundary condition necessary in the simulation of the SOL, the lab frame is the natural choice. The effect of this choice is that the complexity normally arising in various LHS operators (in the Vlasov part) of the electron kinetic equation is transferred into the electron-ion collision operator. In this work the ions are treated as a cold beam for the sake of this collisional computation. While the addition of ion temperature effects would in theory increase the fidelity of the operator, it would strongly depend on the velocity grid resolution around the ion flow velocity. It is likely that the effect would be negligible as the ion thermal velocity is much smaller than the electron thermal velocity, and that the main collisional effect of the ion thermal population is well captured by the delta function of a cold beam.
However, refinement of the electron-ion collision operator is a potential path in the development of the code.
Momentum is not explicitly conserved in the electron-electron collision operator for l > 0 as implemented in the code. As suggested by Tzoufras et al. [23], we ensure total momentum conservation by transfering the momentum lost in electron-electron collisions to the ions. This produces a spurious electric field in the regions of high flow P R E P R I N T speeds. However, this electric field is negligible compared to the physical electric fields occuring in realistic SOL simulations, as the region of non-zero ion flow tends to be in front of the divertor target, where the pre-sheath electric field is generated through pressure gradients.
The way we calculate the electric field also ensures that quasineutrality is observed up to displacement current effects. This allows SOL-KiT to treat collisionless phenomena on the Debye length scale. However, as the logical sheath boundary condition does not require resolving the Debye sheath, this regime has not been fully explored. While SOL-KiT is implicit, in practice the timestep is limited by the capabilities of the solver. We currently use the Bi-CGSTAB iterative matrix solver with Block Jacobi preconditioning in the PETSc package. This solver tends to struggle when the timestep is many (≈ 50) times the electron-ion collisional time, as well as when higher harmonics are included due to the added stiffness of the matrix due to the fast evolution of f l for high l's. However, in most situations we are aiming to resolve the collisional phenomena, and the timestep-limiting effect of higher harmonics becomes evident at a number of harmonics already high enough where the dense matrix effects already cause the majority of the slowdown.
SOL-KiT is currently parallelized using MPI, with domain decomposition along the x-axis. Basic optimization was performed on the code, however, many avenues of potential improvement remain. Besides the already mentioned grouping of excited neutral states, parallelization in l number is also being considered, though the degree of speedup attainable would depend heavily on the inter-processor communication. Different preconditioners and solvers could also be employed depending on the problem, though this would require more in-depth tests of the code's performance.
One of the code's design goals was to enable consistent and convenient one-to-one comparisons between a fluid P R E P R I N T and kinetic model of electron transport. This was accomplished by ensuring that the physics implementated in the kinetic and fluid modes use the same data (e.g. cross-sections etc.). A possible extension of this approach would be the inclusion of various non-local fluid models as comparison options. The easiest to implement would possibly be the flux limiter method [4], and more complicated models like the SNB model and others [43] could be implemented self-consistently with the current framework.
The elastic electron-neutral collision operator currently implemented assumes a simple cross-section (see 2.1.3.).
As we believe the cross-section could be improved with a more detailed model based on experimental data, while tested, this operator is not in regular use. Improvements of the neutral model are being planned, including potential additions of molecules, the implementation of a fluid as opposed to a diffusive neutral model, and the option to include high Z impurities. To accompany the improved neutral model, the next step in SOL-KiT's development will be the addition of an ion temperature equation, as well as an electron-ion collision operator for l = 0. This will allow us to probe more varied SOL regimes, where the ion and electron temperatures are decoupled.
Compared to PIC codes, the finite difference approach taken with SOL-KiT provides a speed-up in computation, as well as avoiding the noise issues present in PIC codes. The closest comparison to SOL-KiT, however, is the code FPI [16,17,18]. While a form of the logical boundary condition appears to have been implemented in FPI [16,17], to our knowledge no in-depth discussion of its Legendre formalism form is available. Similarly, electron-neutral inelastic collisions have been implemented and reported in both ALLA [15] and FPI, but the conservative properties of the operators were not explored in detail. In this paper we report in detail both of these numerical aspects in the context of a Legendre decomposition model. Furthermore, the combination of a fully implicit algorithm together with the lab frame ion treatment and the aforementioned numerical facets has been realized for the first time in SOL-KiT.
The addition of a self-consistent fluid mode for the electrons adds another aspect to the code, with work on coupling the two modes under way.

Conclusion
We have presented the model and the numerics behind the newly developed fully implicit arbitrary harmonic Here it's helpful to visualize the two-dimensional velocity space of the cut-off. This is shown in Fig A. The integral in (A.9) can be found in the literature (see for example [46] for general case) and well known recurrence formulae for the derivatives of Legendre polynomials can be used to get x min l + l + 1 P l (x min )P l (x min ) + 1 (l − l )(l + l + 1) lP l−1 (x min )P l (x min ) − l P l (x min )P l −1 (x min ) , and For P 0 ll we get δ ll if l , l odd and l = l 0 if l , l odd and l l 2l+1 (l−l )(l+l +1) (lP l−1 (0)P l (0) − l P l (0)P l −1 (0)) if l odd and l even It should be noted that for v < v c the above equations give P ll = 0 for all odd l, which is the consequence of reflection symmetry, from which one can also see that there is no effective flux contribution for v < v c .