Kinetic effects in parallel electron energy transport channels in the scrape-off layer

We present an analysis of parallel electron energy transport in the scrape-off layer (SOL), considering the convective and conductive channels, as well as the radiation and neutral inelastic energy transfer channels involving atomic deuterium. Kinetic effects in both equilibria and conductive transients are explored by utilizing the capability of the SOL-KiT code to treat electrons as either a fluid or kinetically. We find kinetic effects in multiple channels, with an emphasis on those occurring during the investigated conductive transients. Energetic electron effects in the heat flux, as well as a modification of ionization rates of up to 40% compared to Maxwellian rates during perturbations in detached conditions, are reported and discussed.


Introduction
Energy transport in the scrape-off layer (SOL) is of great importance in the design of future fusion reactors [1], and the understanding and management of the divertor target heat load is crucial in the success of the ITER and DEMO projects. Energy transport parallel to the field lines in the SOL can occur via a multitude of channels, and a thorough accounting of the importance of each is a necessary step towards the integrated understanding of energy transport. For example, one parallel energy transport channel is the combination of the conductive and convective heat transport, usually treated using fluid modelling, where the classical theory of Braginskii [2] is used and often modified with flux limiters [3]. Another channel is the plasma-neutral interaction and radiation channel, for which Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. the rates are traditionally calculated using collisional-radiative models assuming Maxwellian electron populations [4]. The parallel direction in the SOL is, however, characterized by variations of collisionality, and this can lead to modifications of both heat flux transport and rate coefficients due to kinetic effects [5].
Modifications to various aspects of parallel transport due to kinetic effects have been reported in the literature, ranging from time dependent flux limiters and target plasma sheath properties [6][7][8][9][10][11][12][13], to the enhancement of ionization rates [14]. Conductive transients have also been identified as potentially heavily influenced by kinetic effects [15].
In order to provide a systematic comparison between kinetic and fluid models of electron transport in the SOL, a new transport code has been developed. SOL kinetic transport (SOL-KiT) [16] is capable of treating electrons as either a fluid or using a kinetic model, enabling one-to-one comparisons. A power scan with medium size tokamak relevant parameters was performed using SOL-KiT, and the modifications to the heat flux and sheath properties have been examined previously [17]. Kinetic effects have been observed during both equilibria and conductive transients, and we continue the investigation here by delving into the electron energy channels, with a special focus on kinetic effects on ionization during conductive transients.
The structure of the text is organised as follows. We briefly go over the basics of the model in order to provide a standalone basis for the discussion ahead, while referring the curious reader to more detailed resources. An overview of the simulation data used in the analysis is given before introducing the energy channels explored in this work. We investigate the differences in energy transport between fluid and kinetic electron models in both established equilibria and during conductive transients. Finally, we discuss the possibilities of refining the model in terms of included energy channels.

The SOL-KiT model
A more detailed description of the SOL-KiT model is given elsewhere [16,17]. Here a short overview is given in order to capture the general physical aspects of the model, while omitting implementation details.
We focus on a 1D SOL, straightening out the magnetic field, and taking x = 0 to denote the location of the upstream symmetry plane. The domain ends at the divertor target sheath entrance. The equations solved are the fluid or kinetic electron equations, ion fluid equations, an atomic Diffusive-Reactive collisional-radiative model, and an equation for the parallel electric field. Care was taken to ensure consistency between the fluid and kinetic electron models in order to provide adequate comparison opportunities.
The fluid electron equations are where E is the electric field in the parallel direction and S = S ion + S rec is the ionization and recombination particle source. n e and T e are the electron density and temperature, respectively, and u e is the parallel flow velocity of the electron fluid. We use the Braginskii closure [2,18] R ei = R T + R u , R u = −m e n e 0.51(u e − u i )/τ e , R T = −0.71n e ∂(kT e )/∂x, q e = q T + q u , q u = 0.71n e kT e (u e − u i ), where τ e is the electron-ion collision time [2] and κ e = 3.2n e kT e τ e /m e ∝ T 5/2 e is the classical Spitzer-Härm heat conductivity. R en is the total electron-neutral friction, calculated using a slowly drifting Maxwellian for the electrons. The external heating and the energy exchange with atoms in the temperature equation is given by Q = Q ext + Q en .
For the deuterium ions we take n i = n e , and set T i = T e . Thus the only ion equation that is needed is the momentum equation where R ie = −R ei . Charge exchange friction R CX is given by where the sum is over neutral atomic states, with the simplification of cold ions and cold and stationary neutrals. The constant charge exchange cross sections σ CX,b are approximated by the low energy values from Janev [19].
The electric field is obtained using Ampère-Maxwell's law For a discussion on this form of the electric field equation the reader is directed to our previous paper [17]. The Diffusive-Reactive collisional-radiative model for the neutrals is given by where we use moments of the electron distribution function (Maxwellian in fluid model) to calculate the ionization and (de)excitation rates K, as well as three-body recombination rates α. Hydrogen atomic data are all taken from Janev [19] and NIST [20]. Spontaneous emission rate data A is included up to state b = 20. This truncation should not introduce a substantial error for higher states, as those are primarily collisionally dominated. Finally, we include radiative recombination β as a function of temperature [19]. The 1D diffusion coefficient is simply for the sake of which we treat neutrals as having a thermal velocity v tn = √ 2kT n /m n , where the neutral temperature T n is a free parameter in our model and m n = m i . We assume that the diffusion is due to elastic collisions with ions and ground state neutrals and charge-exchange. Usually negligible compared to charge-exchange, elastic collisions are approximately modelled using σ el = πa 2 0 . At the sheath boundary where we assume ambipolar flux Γ i = Γ e , ions reach the sound speed. The electron sheath heat transmission coefficient [21] in the transmitted sheath heat flux is set to γ e = 2 − 0.5 ln(2π(1 + T i /T e )m e /m i ) when using the fluid model. Neutrals are recycled with flux Γ REC = −RΓ i , where R ≤ 1. For more details on the boundary condition in both the fluid and kinetic models the reader is directed to our code paper [16]. When treating electrons kinetically, we decompose the classical 1D kinetic equation using Legendre polynomials, where the RHS contains all of the collision and source operators. The details of the formalism, based on the KALOS/OSHUN [22,23] spherical harmonic decomposition, are presented in more detail in both the previous study [17] and the dedicated code paper [16]. Here we note only that equation (9) becomes a system of equations for the various distribution function harmonics f l where A l and E l are the advection terms, and C l are all other operators. C l includes Coulomb collision operators for both electron-electron and electron-ion collisions, as well as all of the inelastic electron-atom collision operators corresponding to the various processes in (6). Finally, an isotropic heating source (acting only on f 0 ), as well as an implementation of the logical boundary condition [24] are included.

Simulation setup
For the analysis presented here we use the same runs as those in our previous work [17]. Here a quick overview of the simulation setup is given, and the reader is referred to the code paper [16] for further numerical details.
In the velocity space, 80 cells are used, with cells near the origin being finer. The total length of the velocity space domain is ≈12v th,0 , where v th,0 is the electron thermal velocity for a reference temperature of 10 eV. In this work the diffusive approximation in the kinetic model is used, solving only equations for f 0 and f 1 [25].
The simulation domain is spatially decomposed into 64 cells, with cells closer to the sheath boundary being smaller, to allow for better resolution near the target. The system length is L = 10.18 m, with the external step function heating source localised over a length of L h = 3.75 m upstream. 100% recycling is used to maintain a total line-averaged density This is due to the fluid and kinetic models predicting detachment at different input powers, with the kinetic model predicting a higher input power for the rollover.
Recycled deuterium atoms are assumed to have a temperature of T n = 3 eV (mimicking Franck-Condon enhancement [21]). A total of 30 neutral states are tracked in the simulations. Equilibria used have effective power fluxes ranging from 1 to 6 MW m −2 , with a step of 0.5 MW m −2 . The choice to do an atypical parameter scan, in input power instead of density, comes from the original computational constraints with the code. Recent optimization improvements enable us to perform density scans efficiently as well, and this will be the topic of future work. The conductive transients were simulated by starting from the obtained equilibria, and increasing the input power flux to 45 MW m −2 for ≈10 µs. The perturbation was allowed to relax for a further ≈10 µs by returning the input power to its original value.
In order to aid the following section, we present the comparisons of temperature and density profiles between the two electron models for two of the input powers in figures 1 and 2. We found that the fluid electron runs detach below an input power of 2.5 MW m −2 , while this threshold is moved up to 3.5 MW m −2 when electrons are treated kinetically. This was done by investigating particle flux rollover as well as the pressure drop near/at the target (see figures 6 and 8 in our previous paper [17]). The difference in detachment onset is visible in figures 1 and 2, as the kinetic model has higher densities and lower temperatures near the target when the input power flux is 3 MW m −2 . Note, however, that we did not observe particle loss detachment (hence the higher divertor plasma densities in our runs), likely due to low overall densities and missing molecular physics.   (13) when electrons are treated as a fluid. The channels are: sheath heat flux (q sh ), energy lost to electric field (q E ), energy spent on maintaining ionization and excitation states (qn), and energy lost to atoms that is subsequently radiated (q rad ). The position of the maximum q rad is shown as a dashed black line.

Electron energy transport channels
In order to introduce the electron energy channels considered in this work, we begin with the electron fluid energy equation, derived by taking the energy (second) moment of equation (9), where W e is the electron energy density. The inelastic energy transfer rate Q en can be decomposed into the radiated power density −Q rad and the power density going into ionization and excitation of atomic states −Q n (minus signs introduced for convenience). If we then integrate the energy equation over the entire domain, we get where we have simply used the fact that the upstream boundary is reflective (a symmetry plane), and set the outgoing electron heat flux to q sh . In equilibrium this can be reduced to where the RHS represents the effective input energy flux, and the LHS gives the four major channels considered in this work. These channels are: • The conductive and convective channel q sh -delivered directly to the target by the electrons through the sheath. • The electric field channel q E -energy taken from the electrons by the electric field; transferred to the ions and delivered by them to the target. • The hydrogenic radiation channel q rad -energy radiated by the atomic neutrals; this is delivered to the target and the walls over a potentially larger area than the considered flux tube. • The neutral transfer channel q n -energy transferred to the atoms via ionization and excitation, used to maintain plasma density in equilibrium; this is delivered to the target by the ions and atoms as they impact it.

General fluid vs kinetic channel behaviour
In equilibrium, as stated in (13), the considered channels sum up to the input power, and the behaviour of the channels reflects the SOL regime. The contributions to the various channels as a function of input power when electrons are treated as a fluid are shown in figure 3. The same channel contributions obtained with the kinetic electron model are given in figure 4. The different detachment behaviour observed previously [17] and inferred from figures 1 and 2 is also evident from the radiation channels in figures 3 and 4. Namely, figure 4 shows a later drop in radiated power when electrons are treated kinetically compared to the fluid case. A selection of percentage contributions to the four channels is given in tables 1 and 2 for three of the input powers, including 3 MW m −2 , at which the kinetic runs are detached, whereas the fluid ones are not. This can be seen from the percentage of input power radiated, which is 20.85% in the kinetic case, and 16.07% in the fluid case. Another interesting feature is the greater amount of power lost to the electric field q E in higher input power runs with fluid electrons. However, this is most likely the consequence of higher temperatures in the fluid runs at high powers, as the electric field acts to balance the plasma pressure, and an increase of temperature would lead to an increase in pressure, and thus a greater electric field.
These results show a presence of kinetic effects in equilibrium simulations. However, as noted in the previous analysis of the simulation data [17], kinetic effects are not dominant in equilibria, but manifest themselves strongly during conductive transients, where the behaviour change is more striking. To illustrate this in energy channel terms, we plot similar stacked area plots as figures 3 and 4 for the case of conductive transients. In this case, the ∂´W e dx/∂t term in (12) is not zero, and is not represented in the figures. The simple interpretation is that the energy density of the plasma increases or decreases during the transient, and as such it is just the outgoing energy flux that is captured in the stacked area plots below.
The fluid vs kinetic differences in the perturbation behaviour on the low power background are evident from figures 5  and 6. Firstly, one must note that the peak outgoing power in the fluid case is greater than in the kinetic case, but the outgoing power relaxes more slowly in the relaxation stage of the kinetic perturbation (second 10 µs period). This is likely due to upstream flux suppression, reducing the rate at which the perturbation energy is dissipated when kinetic effects are included. The next kinetic effect visible is the preheat period, showing an increase in the sheath heat flux, that is not accompanied by an increase in radiation. This can be explained by the fact that the fast particles reaching the divertor target after ≈2.5 µs are too energetic to interact with the neutral cloud, and seem to ignore the detached front. In the fluid case of figure 5 the increase in sheath heat flux is preceded by strong radiation  and ionization, while the kinetic case shows a population of fast particles contributing to the sheath flux channel while the radiation and ionization channels slowly increase in strength. The radiation peak occurs later in the kinetic case, and is also lower than in the fluid electron case. This is again a likely consequence of the heat flux suppression, which leads to lower target temperatures in the case of kinetic perturbations [17], as well as a decreased overall ionization during the perturbation, leading to reduced/delayed detachment burnthrough. However, kinetic effects in electron-atom interactions are also present, as will be discussed below.
The high power, attached kinetic background  with respect to its response to a 45 MW m −2 perturbation. Besides the preheat, the kinetic case (figure 8) shows a slightly increased overall outgoing power flux. However, the attached case shows that kinetic effects did not have a significant impact. A likely cause is that the background profile is already quite flat, and that a tenfold increase in input power was not enough to produce gradients steep enough for kinetic effects to dominate.

Conductive and convective electron heat flux channel
The outgoing sheath heat flux contains two components, the conductive and convective heat fluxes. Figures 9 and 10 show the decomposition of the total electron heat flux throughout the simulation domain for two kinetic equilibria. The lower power equilibrium in figure 9 shows a considerable convection contribution around the target, as would be expected in a detached regime. The attached equilibrium (figure 10) shows a much greater conduction dominance up until very close to the target. This is, of course, a consequence of a much narrower ionization region. In both cases, the uniform heating operator causes a spatial ramp up to the input power flux.
During the conductive transients investigated the sheath heat flux is the primary energy channel responsible for the divertor heat load, as can be seen from the total area associated with it in figures 5 to 8. The evolution of this channel is isolated in figure 11. Here the aspects discussed above become even more evident, with the preheat starting from ≈2.5 µs, and the kinetic perturbation showing a lower sheath heat flux peak in the detached case and a higher peak in the attached case. The slower decay of the kinetic perturbation is also evident on the 1 MW m −2 background. Furthermore, the peak heat flux occurs almost 10 µs after the heating is switched off in the kinetic case, further illustrating the heat load 'smearing' effect of  Evolution of the sheath heat flux channel corresponding to the three considered equilibria-warmer colors correspond to higher input power equilibria. The heat load in the lowest power kinetic case is spread out over a longer period of time, with a very fast preheat and a second peak almost 10 µs after the source was switched off. flux suppression. Interestingly, the borderline detached kinetic equilibrium with input power flux of 3 MW m −2 behaves like the fully attached case. This might be caused by the plasma attaching during the preheat phase of this particular transient, as could be inferred from target temperatures in figure 10 of the sister paper [17].

The radiation and neutral transfer channels
The two channels involving atomic neutrals in the model are the radiation channel q rad and the neutral transfer channel q n . When decomposed into their fundamental subchannels, they can be written as where q deex rad and q recomb rad are the radiative deexcitation (line radiation) and radiative recombination (band radiation) contributions, while q ex , q ion , q deex , and q recomb are energy fluxes associated with the collisional processes of excitation, ionization, deexcitation, and three-body recombination, respectively. To illustrate the subchannel decomposition, the 1 MW m −2 kinetic equilibrium is shown in figure 12. As would be expected for an atomic hydrogenic plasma with temperature over 1 eV, the radiation channel is dominated by line radiation (purple line in figure 12), while recombination radiation is negligible. Ionization and excitation together dominate the neutral transfer channel in most of the domain. However, superelastic collisions (collisional deexcitation) become non-negligible close to the target.
The fluid case of the subchannel decomposition (not pictured here) for the background in figure 12 does not differ greatly compared to the kinetic case, supporting the finding by Zhao et al [26]. that kinetic effects on atomic rates in equilibria are negligible.
In order to contrast equilibrium and transient kinetic effects on atomic rates, figure 13 shows the total ionization source S ion at the start and two thirds through the heating phase of the conductive transient on the 1 MW m −2 equilibrium. Plotted is a comparison between the sources calculated with the distribution obtained with SOL-KiT and sources that would result from a Maxwellian electron distribution with the same density and temperature. Here we see again that at t = 0 the  Maxwellian (fluid) model rates agree well with the kinetically obtained ones. At t = 6.68 µs, the kinetically calculated rates are modified compared to what a fluid model would give using the same density and temperature profiles. Note that this is not a comparison between a fluid and a kinetic perturbation. Instead, comparing different methods of rate calculation on the same profile highlights the pure impact of kinetic effects. The distributions of atomic states used for both the SOL-KiT and Maxwellian ionization source calculations are the same, further isolating the effect of rate modification. However, the logarithmic y-axis and the choice of plotted times in figure 13, while enough to demonstrate a deviation from the fluid model, obfuscate the extent of the rate modification. In order to better grasp the evolution of rate modification and its extent, we plot the maximum relative deviation of the ionization source from fluid value at each timestep throughout both the heating and cooling phases of the conductive transient in figure 14. The ionization source deviation during the transient reaches 40%, but quickly relaxes back toward lower values.
For completeness, we plot the same sort of ionization source deviation for the transient launched on the 4.5 MW m −2 background in figure 15. The general trend that kinetic effects during perturbations are weaker on hotter backgrounds continues, as the ionization source never deviates much more than 5% from the Maxwellian value.

Discussion
The current limitations of the SOL-KiT model have been discussed elsewhere in detail [16,17], but we will quickly go over them in order to give context to the following discussion on energy channels. Firstly, the model is 1D, and all magnetic geometry effects and cross-field transport are ignored. Secondly, the T i = T e assumption weakens the model and potentially produces artefacts [17]. This will be improved in a future version of the SOL-KiT code. The neutral model is pending updates as well, and we will discuss this in more detail below. Finally, higher harmonic simulations are planned, as the presented data was obtained using only l = 0, 1. However, preliminary results show that higher harmonics are mainly localised around the target sheath, with the exception of l = 2.
In this work we have presented a number of kinetic effects in the various electron energy transport channels in the SOL. Before discussing the results, it is worth noting the channels that are not present in the model. Firstly, electron-ion energy exchange has been ignored by setting T i = T e , and this channel will be automatically included with the introduction of the ion temperature equation in a future version of SOL-KiT. We expect this to have an effect on high power equilibria and on all transients, as the ions and electrons will be allowed to decouple. It would be interesting to compare the electron-ion channel to the channels presented here, especially if hot neutrals are allowed to deposit their energy into the cold plasma close to the target, through either electronneutral or ion-neutral elastic (or charge-exchange) collisions. Secondly, molecules and impurities are not yet included in the SOL-KiT model, and these would provide further radiation and neutral transfer channels. Other channels might include bremsstrahlung radiation, as well as effective cross-field transport. Nonetheless, the presented comparisons between fluid and kinetic models with a restricted set of channels should still serve well to highlight kinetic effects.
We have presented the different channel contributions for both the fluid and the kinetic electron model throughout the performed power scan. Figures 3 and 4 and tables 1 and 2 highlight the differences in energy channels, and these were associated with a difference in detachment behaviour between fluid and kinetic models [17]. Namely, at an input power flux of 3 MW m −2 , the kinetic equilibrium is detached, and the fluid equilibrium is attached, which is reflected in the drop from ≈21% to ≈16% radiated power when moving from a kinetic to a fluid description. However, differences between the kinetic and fluid model become much more striking in the case of the examined conductive transients, when the input power flux was increased to 45 MW m −2 for ≈10 µs. Multiple kinetic effects have been observed during these transients, and their dominance increases as the input power of the perturbed background is decreased, being most pronounced at the lowest power equilibrium considered. Preheat is evident in both low and high power cases, as shown in figure 11, as hot electrons start depositing energy through the sheath heat flux channel from ≈2.5 µs. Another important effect that was observed was a spreading out of the heat flux over a much longer period in the kinetic case compared to the fluid case, while the peak value was almost halved. Further research is necessary to explore this effect in different regimes.
While no significant kinetic modification of atomic rates was observed in equilibria, consistent with what is claimed in the literature [26], we have found ionization rate modification during transients. Figure 14 shows that, during the transient evolution on the lowest input power background, ionization rates deviate up to 40% from those predicted by using a Maxwellian distribution. This effect is much weaker when the transient is launched on an attached background. We note that the overall ionization rates in the kinetic case are lower, due to flux suppression, and this leads to more of the detachment front surviving the immediate transient. Non-local effects on atomic rates during transients have been previously explored in the literature [14], but the simulated problem was a conductive relaxation of a ramp in temperature and ionization degree, where a much stronger modification was observed. We believe that in order to observe similar strong modifications of hydrogenic ionization rates it would be necessary to explore convective transients, or to use a greater power increase in our conductive transient simulations, while utilizing a detached background. Further study is needed in order to push the non-local atomic rate effects in hydrogen to their limit in situations relevant to the SOL.
Another place to look for kinetic effects in both equilibria and transients is in impurity ionization and excitation rates. With the higher energy thresholds for collisional processes of electrons and heavier atoms, we can expect that energetic electrons in the tail of the distribution (that would normally ignore detachment due to their low collisionality with hydrogen) can interact with impurities and molecules, causing non-local ionization effects. An ongoing research project is being devoted to the inclusion of these effects into SOL-KiT in order to utilize the model comparison capabilities of the code in the systematic exploration of kinetic effects in atomic and molecular processes in the SOL.

Conclusion
We have presented an analysis of parallel electron energy transport in the SOL, with a focus on kinetic effects in equilibria and transients. Decomposing the energy transport into channels, we explored the importance of each channel, as well as the differences between a kinetic and fluid description of electrons.
While some kinetic effects can be observed in equilibria, namely flux suppression and different detachment behaviour, we report that kinetic effects become dominant in conductive transients investigated in this study. Energetic electrons ignore the detachment front and directly contribute to the sheath heat flux in the early phases of kinetic simulations, while an overall 'smearing' of the heat flux over a longer period compared to the fluid model is observed. We also report up to a 40% increase in ionization rates during transients, compared to fluid predictions with the same density and temperatures. However, due to flux suppression and 'smearing', overall ionization is not as high in the kinetic transient simulations as it is in the unsuppressed fluid case, leading to a smaller/slower loss of detachment.
We have shown that kinetic effects during transients dominate energy transport in multiple channels, and have discussed potential new avenues along which to expand the presented work in the future.

Acknowledgments
This work was supported through the Imperial President's PhD Scholarship scheme. This work was partially funded by the RCUK Energy Programme [grant number: EP/T012250/1]. Simulation results were obtained with the use of the Imperial College Research Computing Service [27].