Two-dimensional collisional particle model of the divertor sheath with electron emissive walls

A novel two-dimensional particle-in-cell (PIC) code, named Divertor Edge Simulator of Plasma-wall Interaction with Consistent COllisions (DESPICCO) and developed at CNR-ISTP, is capable of simulating the thin plasma layer of several millimeters, adjacent to the divertor tiles of a Tokamak fusion reactor. Here, kinetic effects and non-neutral plasma physics in the Debye sheath can be self-consistently captured by the PIC approach. The code is firstly benchmarked against literature one-dimensional codes and additional theoretical predictions for a magnetized sheath. Then, it is applied to a realistic divertor scenario featuring an attached plasma with monoblocks (MBs) radial misalignment and gaps, to compute the energy flux amplification factor at the exposed MB edge. A non-ambipolar local current density close to the leading edge and an average sheath heat transmission coefficient larger than the one predicted by classical sheath theory, are found. The effects of electron wall emission and plasma-gas collisions on the ion Mach number and on particle and energy fluxes to the walls are finally estimated to determine future guidelines for simulations. Ion collisions with recycled neutrals and both secondary and thermionic electron emission from the wall are found to have a relevant impact, with the overall effect of reducing by 25% the average ion impact energy, and by 15%–20% the total heavy particles energy flux to the walls, with relevant implications on the divertor wall erosion.


Introduction
One of the most important open issues in the development of a thermonuclear tokamak fusion reactor is represented by the * Author to whom any correspondence should be addressed.
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. heat and particle fluxes on the divertor plates. The divertor is a system to create a diversion of the magnetic field confining the plasma; it is the most reliable alternative to the limiter because it prevents contamination of the bulk plasma by metal impurities, being placed at a higher distance from the core [1,2]. The divertor wall is made of monoblocks (MBs), whose surfaces intersect the open magnetic field lines of the reactor and are therefore directly exposed to the plasma. The thin plasma layer adjacent to these surfaces is interested by several physical phenomena [2,3]: (i) formation of a non-neutral sheath and a quasi-neutral magnetic pre-sheath; (ii) neutrals recycling from the wall due to recombination of impacting ions; (iii) secondary electron emission (SEE) due to both electron and ion impacts; (iv) thermionic electron emission (TEE), (v) inter-particle collisions, with a major relevance of ion/electron collisions with the recycled neutrals, and (vi) ion sputtering, i.e. the surface erosion with the release of wall atoms due to high-energy ion impacts. Moreover, the gaps existing between the MBs make two-dimensional effects, like peaks in the heat flux at exposed edges, very relevant.
In the last decades, many authors in literature have simulated the plasma-wall transition through a magnetized plasma sheath. Chodura [4], in 1980s, generalized the Bohm's criterion for the stability of the magnetized collisionless presheath, stating that ions must be at least sonic in the direction parallel to the magnetic field at the pre-sheath entrance. With the use of a 1D particle-in-cell (PIC) code, Chodura obtained the dependency of pre-sheath and Debye sheath widths and of the corresponding potential drops on the angle of the magnetic field relative to the surface. To first order, the potential drop across the combined Debye sheath and magnetic pre-sheath only depends on the ion to electron mass ratio, as predicted by simple fluid models [2,4]. However, the simulations [4] have shown that it also depends on the ratio between the electron cyclotron frequency and the plasma frequency, thus highlighting the importance of using kinetic models to accurately reproduce the plasma-wall interaction when magnetized sheaths are involved. Riemann [5] later extended Chodura's model by also accounting for collisions. More recently, Khaziev and Curreli [6] have studied the plasma-wall transition at arbitrary angles between the magnetic field and the material surface with both a 1D fluid model and a 1D PIC code, showing that a correct quantitative prediction of the ion energy and impact angle distributions at the wall requires a kinetic treatment. Such distributions are particularly relevant for the induced wall sputtering and can deviate significantly from the optical (impact along the magnetic field direction) or fluid (impact along the average ion velocity) approximations, even in collisionless scenarios, as also predicted by ion orbit studies [7]. A recent work by Thompson et al [8] includes the effects of ion/neutral collisions in a similar 1D PIC code.
Coming now to divertor scenarios applications, Tskhakaya [3] and Taccogna et al [9] studied the plasma-wall transition using a 1D PIC code, including particle collisions with an Monte Carlo Collisions (MCC) model. The former showed that collisional effects can induce large deviations in potential profiles and energy fluxes to the walls, when compared to classical sheath model predictions. The latter, on the other hand, included a large variety of collisional processes, including molecular assisted recombination (MAR) reactions [10,11] and concluded that MAR reactions can be very effective in high recycling regimes. In this respect, Kukushkin et al [10] suggested that, in deuterium discharges, MAR based on the D − precursor can be particularly effective. Finally, Komm et al [12] recently presented the first two-dimensional PIC study of the plasma-wall interaction at the divertor, considering a collisionless plasma and a realistic MBs geometry. For an International Thermonuclear Experimental Reactor (ITER) divertor geometry, a comparison of the particle and energy fluxes obtained by a PIC code with those predicted by ion orbit and optical methods showed relevant discrepancies, with the former predicting lower peak values at exposed MB edges. These differences are due to the self-consistent electric field considered by the PIC and are expected to grow if collisional effects are also included. Recent PIC studies [13][14][15][16] have also dedicated particular attention to the effect of the electron wall emission on the ion flux to the divertor wall. Thermionic and electron-induced SEE processes contribute to the reduction of the sheath potential drop, and in extreme cases, to the formation of a potential well or virtual cathode and inverse sheath [13,17], when either the wall temperature is very large or when the average SEE yield becomes larger than 1. However, both space charge effects and the prompt re-deposition [14] due to electron gyromotion have the effect of generally inhibiting the virtual cathode formation.
PIC codes have therefore been used extensively in the literature to study the plasma-wall interaction, as they enable a self-consistent calculation of the transport coefficients, which fluid codes necessarily assume as input data. However, to the authors knowledge, no fully two-dimensional kinetic simulation including all relevant physical processes (e.g. collisions, and wall electrons emission) has been carried out so far. This represents the only way to self-consistently assess the total energy fluxes to the divertor MBs, and, more importantly, assess the local erosion profiles, which are very dependent on the ion impact energy and angle distribution function. The goal of this work is then to present a new two-dimensional PIC-MCC code, DESPICCO ('Divertor Edge Simulator of Plasma-wall Interaction with Consistent COllisions'), capable of simulating the above mentioned kinetic effects for a realistic divertor MBs geometry. After validating the code with dedicated benchmarks against theoretical and 1D PIC predictions [2,18,19], the code is used to reproduce a realistic ITER scenario, featuring attached plasma conditions and an exposed edge due to the radial misalignment of the MBs. Finally, the effects of including both collisions and electrons wall emission on plasma properties and particle/energy fluxes are assessed quantitatively.
The paper is organized as follows: section 2 describes the structure of the code and presents its most relevant models. Section 3 describes the benchmark simulations, while section 4 presents the realistic ITER scenario simulation. Section 5 reports the results of a parametric analysis on the effects of electron-wall emission and particle collisions, and, finally, conclusions are drawn in section 6.

PIC model description
The PIC approach is a hybrid Lagrangian-Eulerian approach [20,21]. On one hand, it discretizes the distribution functions of ions and electrons using computational macro-particles (each one representing a large number of elementary particles) subject to Newton's equations. On the other hand, it solves for the self-consistent electric and magnetic fields at the nodes of a mesh, after macro-particles charge and current have been deposited there [18]. The model presented here is a Cartesian electrostatic 2D-3V PIC: it considers a 2D (y, z) space domain representing a thin layer attached to the divertor tiles, as shown in figure 1, and 3 velocity components for the macro-particles, needed to compute correctly the Lorentz force and drift/collisional effects. Furthermore, the magnetic field is externally imposed since the plasma current involved is considered to be small enough to induce a negligible self-consistent magnetic field (compared with the external one).
Referring to figures 1 and 2, macro-particles are first moved knowing the fields at the grid nodes, using a standard leapfrog technique [18] based on the Buneman-Boris algorithm. The particle push also performs the periodic reflection of macro-particles exited through the lateral y min , y max boundaries, and a thermal refluxing for those exiting from the upstream boundary. Then, plasma-wall interaction is simulated, taking into account several phenomena, such as SEE from electron impacts and thermionic emission, whose models are described more in depth in section 2.3. As shown in figure 1, all ion and electron macro-particles hitting the walls or exiting through the bottom boundaries are removed from the simulation (electrons emitted by the divertor surface are in fact treated as new particles). Collisional events are then simulated in each PIC mesh cell, following the algorithms described in section 2.4. After removing all macro-particles lost due to wall impact and/or collisional events, the injection of new macro-particles is carried out either from the upstream boundary surface or from a volumetric region attached to it. Finally, macroparticles are weighted to the nodes of the mesh to obtain their bulk properties and the charge density, used by a Poisson's solver to obtain the electric potential and hence the electric field at the mesh nodes. Finally, after some post-processing operations, the loop advances to the following time step. In the following sub-sections, the main algorithms of the PIC loop are described in detail.

Particles injection
In order to maintain stationary plasma conditions, macroparticles must be injected into the simulation domain. Here, two injection approaches have been considered. The first one, 'type 1', consists of injecting a number of ion/electron pairs equal to the number of impacting ions on the divertor tiles and on the bottom boundary (see figure 1). Therefore, it requires to initially load the simulation domain with a macro-particles population resembling the targeted average densities. The second injection approach, 'type 2', consists of the injection of a constant prescribed flux of ion/electron pairs, and enables simulations starting from a vacuum. Both approaches have been considered in the literature [22]. While injection type 1 permits to obtain self-consistent fluxes to the wall, knowing bulk properties upstream, type 2 permits to study cases with a known fixed particles inflow. As mentioned earlier, all ion and electron macro-particles that cross the upstream boundary are refluxed, i.e. re-injected downstream from a prescribed velocity distribution function. This is the same for both injection and refluxing and consists in a drifting Maxwellian flux distribution function, as used in [19] and as shown in figure 1. In particular, the source probability distribution function S i,e for injecting an ion/electron macro-particle with velocity v is: where T i,e is the ion/electron temperature and u i,e is the ion/electron injection fluid velocity vector. In the simulations The proportionality of the source distribution on the z component of the velocity (or normal component to the upstream boundary) is necessary to reproduce the correct volumetric distribution function, as explained in [19]. For the sake of clarity, figure 3, compares the resulting distribution function of electrons next to the injection surface on the top z boundary, in a simulation in which electric fields are set to zero on purpose, and assuming either the source distribution of equation (1) (with u e = 0), or the simple Maxwellian distribution function (without proportionality with v z ). Clearly, the latter injection is not capable of reproducing the desired half-Maxwellian distribution and actually tends to a Dirac delta distribution, as the simulation advances in time: particles with zero z velocity are, in fact, the most probable and, in the absence of electric field, they remain in the same location and accumulate there.

Particles weighting and Poisson's equation solver
Macro-particles are weighted to the PIC mesh nodes through a standard area weighting technique based on the Cloud-In-Cell (CIC) method [20], i.e. using a 2D bilinear interpolation shape function. This means that the charge density at a generic node g is obtained as:  (1), and blue lines to the distributions (at successive time instants as shown by the arrow) obtained using a simple Maxwellian source function.
where the summation extends to all N macro-particles located within a cell size distance along both y and z. ∆y, ∆z are the uniform mesh spacings along y, z, e is the elementary charge, W p is the macro-particle weight (per unit length along x), Z p is the charge number of the considered p th macroparticle (−1 for electrons, +1 for singly charged ions), and η ′ p = (y p − y g ) /∆y, ζ ′ p = (z p − z g ) /∆z are the computational coordinates of the p th macro-particle relative to the considered grid node. The same shape function S(η ′ p , ζ ′ p ) is used also for depositing other plasma bulk properties (like the mean quadratic species velocity or the plasma species densities and vector fluxes). With the knowledge of the charge density at the mesh nodes, a Poisson's equation solver, based on PETSc libraries and using the HYPRE preconditioner [23,24], obtains the electrostatic potential ϕ and field E = −∇ϕ. The boundary conditions for Poisson's equation are illustrated in figure 1, and consist in homogeneous Dirichlet conditions on the divertor MB nodes (ϕ = 0), and homogeneous Neumann conditions (∂ϕ /∂z = 0) at the upstream top boundary. Periodic conditions are finally imposed on the left and right domain boundaries. Central finite differences schemes are used for the discretization of the equation.

Electron wall emission model
The current code version implements two electron emission processes; TEE and electron-induced SEE.
For thermionic electron emission, the Richardson-Dushman formula is considered for the emitted electron current density j e,th , which depends on the wall temperature T w and on the material work function W f : For an atomic polycrystalline clean tungsten [14], W f = 4.55 eV and A eff = 60 A cm −2 K −2 . Thermionic electrons are emitted following a Maxwellian flux distribution (equation (1) with v z now referring to the velocity component normal to the wall) with temperature T w .
Regarding SEE, the same model implemented in the SPICE2 code by Tolias et al [15,16] is here considered. This distinguishes between back-scattered electrons, which feature a probability η bse , and true secondary electrons, produced with a yield δ tse . For each impacting electron, featuring an impact energy E p and angle θ p (relative to the surface normal direction), an average number σ see = δ tse (E p , θ p ) + η bse (E p , θ p ) of electrons is emitted from the surface. Back-scattered electrons are emitted according to the Jablonsky-McAfee energy distribution [25,26], with a backscattering direction given by the specular reflection direction in 1/3 of cases and a diffuse Lambertian direction in the remaining 2/3 of cases. True secondaries, on the other hand, are emitted according to the Chung-Everhart energy distribution [27], with an average emission energy equal to 2W f and a diffuse Lambertian direction. The reader is encouraged to read [15] for the details on the computation of η bse and δ tse . These are assumed to tend to 0 as E p → 0, although other references in literature [28,29] seem to suggest larger yield values at low impact energies. Although these discrepancies clearly deserve to be further investigated in the future, here Tolias's approach is preferred as it permits a  [31] more direct comparison with the already existing SPICE2 PIC code.

Inter-particle collisions model
DESPICCO is capable of simulating collisions between ions/electrons and the neutral atoms or molecules produced at the walls due to recycling of the main ionized species (deuterium or hydrogen). Currently, neutrals are considered to be a uniform background of density n n , with a non-drifting Maxwellian velocity distribution at a given temperature. The collisional processes listed in table 1 are considered, with the corresponding cross sections reported in figure 4. An MCC algorithm, with the null collisions method to reduce the number of required random numbers, is here implemented [32,33]. For each PIC cell with a non-zero number of colliding macro-particles (e.g. ions or electrons), this consists in two main stages: a sampling phase, and a postcollisional velocities generation phase.
Regarding the sampling phase, a total cross section σ tot summing all available processes for the considered species (e.g. elastic, electronic state excitations and ionization cross sections for electron-atom, momentum and charge exchange cross sections for ion-atom) is firstly computed for the centerof-mass energy range of figure 4. Then, in each PIC cell and every N coll time steps, a number of macro-particles is selected to check for collisions, considering a maximum collision probability computed as: where (σ tot v rel ) max is the maximum value of the product between the total cross section and the relative velocity over the considered energy range, n n is the density of the background neutrals, and ∆t is the PIC time step. N coll is chosen so that p coll,max < 0.01, which, for the considered neutral background densities and collisional processes, yields N coll,i = 2000 and N coll,e = 1000 for respectively ions and electrons. Then, the collisional reaction effectively taking place for each selected macro-particle (including the null-collision outcome) is evaluated randomly. Firstly, a random velocity is generated for the neutral background macro-particle, thus permitting to compute the relative velocity v rel = v i/e − v n and the center-of-mass kinetic energy E k,com = 1/2 m red v 2 rel , with m red = m i/e m n /(m i/e + m n ) representing the reduced mass for the considered collision. Secondly, the probability of each available reaction r is evaluated as: Note that p null = 1 − r p r > 0, so that there is a non-zero probability that the selected macro-particle does not suffer any collision (the null collision outcome). Once the collisional reaction has been randomly selected (based on the probabilities of equation (5)), the new velocity of the colliding charged macro-particle is computed. For elastic processes, the new velocity is computed assuming an isotropic scattering of the relative velocity vector. If the center-of-mass velocity vector is v com = m i/e v i/e + m n v n / m i/e + m n , and v ′ rel is the newly scattered relative velocity vector (isotropically distributed in direction with |v ′ rel | = |v rel |)), then, the new ion/electron macro-particle velocity is obtained as: For the atomic electronic excitation reactions between electrons and neutrals, the same algorithm is applied, but the initial electron velocity v e (used to compute both v com and v ′ rel ) is firstly reduced in magnitude taking into account the threshold energy loss of the reaction: where E thr is the threshold energy of the reaction (see table 1) and E k,com is the center-of-mass kinetic energy prior to collision. For the ionization reaction, no ion/electron macroparticles pair is generated (which would be inconsistent with the injection strategy '1', see section 2.1), and the same approach of excitation reactions is followed. Therefore, ionization is here modeled as a simple energy sink for electrons with no distributed volumetric source of ions/electrons. Finally, concerning the charge-exchange reaction between ions and neutrals, the ion macro-particle is directly given the neutral macro-particle velocity, thus considering a pure charge-exchange without any momentum exchange.

Code benchmarking
The code has been benchmarked against an existent 1D unmagnetized PIC code [19] that uses the same injection strategy, and against both theoretical and 1D PIC predictions for magnetized sheaths and pre-sheaths [4].
The first benchmark features an unmagnetized collisionless plasma scenario with a floating wall (at z = 0.1 mm), and input parameters listed in table 2 (BM-1 column). The goal is to match the results obtained with a 1D PIC code by Procassini et al [19], considering the same plasma and numerical conditions. In particular, the injection strategy N.1 is considered, so that any ion hitting the downstream wall is re-injected as an ion-electron pair from a volumetric injection region, located next to the upstream boundary and extending 5 mm in height. Furthermore, an average plasma density of 1.5×10 16 m −3 (producing a peak density of 2.2×10 16 m −3 inside the injection region), equal ion and electron temperatures (50 eV), a zero injection drift velocity for both particles species, and an ion to electron mass ratio of 100, are further assumed. The resulting z-profile (direction normal to the wall) of the timeaveraged electric potential, at steady state conditions, is reported in figure 5 for both codes, showing a very good agreement. Both the Debye sheath drop close to the wall, and the potential drop across the injection region, which develops to make ions reach a sonic condition and hence a stable sheath, can be clearly appreciated.
The second benchmark, BM-2, consists of a set of 5 collisionless simulations at varying B field angle relative to a horizontal non-emitting divertor surface (refer to the α angle of figure 1). Input parameters are summarized in the BM-2 Table 2. Input parameters for simulations: (BM-1) injection/refluxing algorithms, (BM-2) B-field angle parametric simulations, (DIV) realistic ITER divertor scenario, and (PWI-Colls) plasma wall interaction and collisions parametric simulations. The simulated elementary masses for ions and electrons are expressed respectively in units of real hydrogen ion (m H ) and electron (me) masses and determine the simulated particles trajectories. The reported macro-particle weight Wp is the number of elementary particles represented by a computational macro-particle (per unit length along x) and only affects the statistical PIC noise.  Finally, the magnetic induction field intensity, 2.94 T, has been chosen to reproduce a cyclotron to peak plasma frequency ratio equal to 1 (necessary to compare the obtained results with [4]).
The z-profile of the electric potential for the various α values is shown in figure 6(a). As expected [2], the total potential drop (from top boundary to wall) is only dimly dependent on the angle between the wall surface and the magnetic induction field. However, as α decreases, the size of the presheath clearly grows. The z-profile of the parallel (to B) ion Mach number, M ∥,i , is reported in figure 6(b). Given the considered ion sonic injection, this is equal to 1 at the upstream boundary, thus satisfying the Bohm-Chodura condition [4] and avoiding the formation of an artificial source potential drop. The quasineutral magnetic pre-sheath edge can be here identified as the point at which M ∥,i begins to increase significantly from 1. This occurs earlier (i.e. at larger z), the lower α, meaning that the magnetic pre-sheath becomes wider. The corresponding profile of the wall-perpendicular ion Mach number, M z,i , is shown in figure 6(c). At the Debye sheath edge entrance, on the other hand, M z,i reaches -1, as requested by a stable monotonic plasma sheath. This point, as expected, becomes closer to the wall as the magnetic field becomes more parallel to it: as a consequence, also the potential drop inside the non-neutral Debye sheath decreases with α. In order to validate the simulation results, the potential drops obtained with DESPICCO are finally compared with those of [4] in figure 6(d). From ion and electron continuity equations, and assuming Boltzmann electrons, it is straightforward to obtain a theoretical evolution of the magnetic pre-sheath potential drop [2,4] as: T e e ln(sin(α)).
This theoretical evolution is shown by the dashed blue line in subplot (d). On the other hand, the total potential drop across pre-sheath and Debye sheath is, to first order, only dimly dependent on the angle α, and can be approximated as ∆ϕ tot,f ≈ −(T e /2e) ln [2π m e /m i (1 + T i /T e )]. Chodura [4], however, provides the exact evolution of the total potential drop (red dashed lines), obtained using a 1D PIC code for the considered plasma conditions (hydrogen ion plasma with cold ions and a cyclotron to plasma frequencies ratio of 1). Finally, the Debye sheath potential drop is obtained as the difference between the total potential drop and the magnetic pre-sheath drop: ∆ϕ DS = ∆ϕ tot,f − ∆ϕ MPS . An excellent agreement can be clearly appreciated between the DESPICCO simulation results and [4].

Simulation of a realistic divertor scenario
A realistic ITER divertor geometry is now considered, featuring an exposed divertor MB edge to the incoming particle flux. This is a worst case scenario in terms of heat flux, that could occur due to an excessive major radius misalignment between successive MBs in the toroidal direction. Referring to table 2 (DIV column), the simulation box extends 9.5 mm in the toroidal direction y and 2.1 mm along z, with an MBs height difference of 0.3 mm and a 0.5 mm gap. Particles are injected with a fixed flux within a rectangular box adjacent to the upper boundary with 0.3 mm height. Half-Maxwellian electrons at T e0 =47 eV and sonic ions (parallel to the magnetic field) at T i0 =4.7 eV are considered. Inside this injection region, electrons are further thermalized (i.e. their velocity is re-generated according to the prescribed distribution) with a frequency of 10 GHz. This artifact has negligible effects on the plasma properties outside of the injection region and allows to increase the numerical stability of the simulations, enabling a slightly larger cell size and hence a reduced computational cost. The magnetic field is inclined by 3.2 deg with respect to the horizontal and goes from right to left towards the MBs (α = 176.8 deg). Neither collisions nor electron wall emission are included. Finally, an artificially increased electron mass, equal to 9.1 m e , and the real hydrogen ion mass have been considered. The spatial gradients are thus captured correctly (since they depend on the ion Larmor radius), while the total potential drop and hence the total sheath heat transmission coefficient are underestimated (due to the lower ion to electron mass ratio). The plasma conditions described above are representative of the worst point of an attached divertor scenario, which can be considered as a sizing conservative case, as it features the largest particle energy fluxes to the walls [34].
The obtained electric potential 2D map, normalized with the bulk electron temperature T e0 , is shown in figure 7(a). The particle injection region can be appreciated in the figure, with the potential peaking at its end, i.e. at z = 1.8 mm. The corresponding electron density and electron temperature maps are shown in figures 7(b) and (c). In both maps, the physical magnetic shadow region casted by the right MB corner in the gap region is evident. In fact, in the absence of collisions, no electrons can traverse the magnetic induction field. Figure 7(d) finally shows the ion Mach number and the ion streamlines: ions become sonic in the direction parallel to B soon after the injection region and further accelerate towards the wall due to the potential gradient within the quasineutral magnetic presheath and non-neutral Debye sheath. Only less than 0.5 mm from the walls, do ions streamlines deviate from the magnetic induction field. Figures 8(a) and (b) show respectively the normalized total energy flux to the walls (with the total energy flux on the undisturbed divertor surfaces, away from the corners), and the corresponding particle flux to the divertor, versus the wall curvilinear coordinate s, distinguishing the contribution of ions and electrons. The normalization heat flux is around 11 MW m −2 , and the ion and total energy fluxes are compared with those of the PG.1 case reported by Komm et al [12] (using the SPICE-2 code). The amplification factors (around 17) predicted by DESPICCO and SPICE-2 clearly match very well. Along the vertical surface of the left MB, the ion energy flux varies almost linearly with the distance from the leading edge as a result of their finite Larmor radius orbits, while electrons have a nearly constant energy flux. Regarding the particle fluxes to the walls (subplot (b)), it is interesting to notice that, close to the exposed leading edge, the current flow is locally not ambipolar, with ions having a generally larger flux on the horizontal surface and electrons dominating on the vertical surface. This is probably a filtering effect of the divertor walls due to the ion gyration radius and the abrupt changes of the electric field map on the two sides of the leading edge. In any case, at a distance of around 1 mm from the exposed edges, ambipolarity applies again. The average ion and electron impact energies are then shown in figure 8(c): ions have a larger impact energy than electrons, since they are accelerated within the non-neutral Debye sheath, and they reach peak values (of around 4 T e0 ) on the vertical MB surface, close to the point where the magnetic shadow intersects the wall (within 1 ion Larmor radius from it). This point can be identified as the location where the electron mean impact energy suddenly drops to zero (no electrons can traverse the magnetic shadow region in a collisionless plasma). The profile of the total sheath heat transmission coefficient γ along the divertor wall is finally shown in figure 8(d), and it is computed as Note that, in equation (9), we are assuming that, in a steady state, the energy transferred to the walls equals the energy transferred to the sheath, and we are neglecting the recombination heat in the evaluation of the energy fluxes. The variation of γ along the MB surface is mainly due to the non-ampipolar character of the impacting particle flux: it shows peaks where the ion flux contribution to the total flux is small, and minimum values where it is dominant. The average sheath heat transmission coefficient (weighted with the impacting ion flux) is equal to 4.8. This is larger than the one predicted by classical sheath theory [2] for a collisionless plasma with a totally absorbing cold surface, and for the considered ion to electron mass ratio of the simulations (m sim i /m sim e = 200) : This result clearly highlights the need of using a kinetic model to obtain appropriate boundary conditions for scrape-off-layer fluid codes.

Effects of electron emission and collisions
In this section, the effects of wall electron emission, due to secondary and thermionic electrons, and of collisions between ions/electrons and a neutral atoms background, are analyzed to identify relevant trends and guidelines. In particular, a background of hydrogen atoms with a uniform prescribed density n n = 5 × 10 19 m −3 is considered. An approximate value for this neutral density can be obtained by assuming a perfect thermal accommodation of impacting ions with the divertor walls (so that emitted neutrals are in thermal equilibrium with the walls) and equal ion influx and neutral outflux, which is a reasonable assumption in steady state operations with a negligible ion implantation rate. For emitted neutrals at the wall temperature T n = T w = 3000 K and for an incoming ion flux Γ i = 3.6 × 10 23 m −2 s −1 (see table 2), this yields  Table 4. Parametric simulations with wall-emission and collisions, showing the normalized plasma potential ϕp (relative to the wall at the end of the injection region, z = z inj ), the parallel ion Mach number M ∥i at z = z inj , and the peak and average energy flux deposited to the wall for ions and electrons.

Considered effects
where the neutral fluid velocity has been approximated with the thermal velocity. The chosen neutral density (5 × 10 19 m −3 ) is lower than this prediction to account for a non-perfect accommodation of recombination neutrals, which yields a larger neutral outflux velocity, and is consistent with the expected range of neutral densities at the divertor walls, for attached plasma conditions [34]. The real masses of both electrons and ions are considered, with a divertor geometry similar to that of section 4, but with a lower MBs height difference (0.1 mm) to yield a lower computational cost (in fact, the extension of the magnetic shadow region along the y direction depends linearly on this). Ions are now injected with half of their sonic velocity, to clearly observe the effects of collisions and wall emission on their required velocity for a stable sheath and pre-sheath. The simulation parameters are listed in table 2, column 'PWI-Colls'. Table 3 provides the estimated mean free path for the considered collisions against the neutral background. These values must be compared against the length covered by particles along the magnetic induction field while traversing the simulation domain (from the end of the injection region). Neglecting the effect of the gyromotion on the distance travelled (good approximation for ions), this is approximately 20 mm. Table 4 reports the plasma potential and the parallel ion Mach number at the end of the injection region, and the peak and average energy fluxes of ions and electrons, for all the considered simulation cases.
The 2D maps of electron density and ion Mach number with streamlines are shown in figures 9(a) and (b), for the case with all effects included. If we compare these results with those of figures 7(b) and (d), the most visible effects are the increase of the peak plasma density upstream (for a nearly equal particle flux towards the walls), the reduction of the ion Mach number, and the existence of a non-zero electron density within the magnetic shadow in the gap region. Figures 10(a)-(d) show the evolution, along the vertical dashed line of figure 9, of the parallel ion Mach number, the electric potential, and the electron and ion temperatures. From subplot (a), it can be seen that the required ion Mach number for a stable pre-sheath diminishes due to the effect of the ionneutral collisions, while electron collisions and wall-emission have negligible effects, except very close to the wall. This ion Mach number reduction has been observed, with quantitative similarities, by Thompson et al [8] in a recent publication, and predicted by collisional pre-sheath theory [5]. Regarding the electric potential, on the other hand, collisions have the effect of increasing the plasma potential by around 15 V. Electronwall emission due to both SEE and TEE, however, reduces this potential drop by a larger amount, so that, when taking all effects into consideration, the potential drop decreases with respect to the case with no effects. This has important implications on the ion sputtering of the divertor MBs, as this depends on the ion impact energy. The electron temperature (subplot (c)) is only dimly affected by collisions (given the large mean free path for the considered electron collision processes, see table 3). Nevertheless, simulation results show that the electron temperature maintains relatively large values to inhibit efficient electron-ion radiative and three-body recombination reactions (which would require electron temperatures in the order of a few eV). Finally, the ion temperature (subplot (d)) is mainly affected by ion-neutral collisions. Away from the wall, these act as an energy sink, yielding to a reduction of the ion thermal energy. In the very proximity of the walls, on the other hand, their scattering effect on the ion velocities produces a  small temperature raise, and hence a more isotropic ion distribution at wall impact. Figures 11(a)-(f ) show the profiles along the wall (zoomed around the exposed leading edge on the left MB) of electron and ion total energy and particle fluxes, and electron and ion average impact energies. Clearly, the inclusion of electron wall-emission has the effect of increasing the electron energy flux everywhere, because of the reduction of the plasma potential relative to the wall (see figure 10(b)). The average and peak ion energy flux, on the other hand, is reduced by ion-neutral collisions and electron wall emission. In particular, ion-neutral collisions have the effect of reducing by approx. 25% the peak ion energy flux on the vertical surface of the divertor, because of the scattering of ion velocities and the creation of slow ions. This reduction of ion energy flux, however, is partially compensated by an increase of high energy neutral atoms (generated by collisions with ions), which is not displayed here since neutrals are treated as a frozen background, and not as Figure 11. Evolutions along the wall coordinate s of (a) electron energy flux, (b) ion energy flux, (c) electron particle flux, (d) ion particle flux, (e) average electron impact energy, and (f ) average ion impact energy. The s coordinate grows from left to right and assumes value 1.6 mm at the exposed leading edge. Results are time-averaged over 50 000 PIC time steps. a particle species. Nevertheless, these neutrals are not accelerated by the non-neutral Debye sheath; therefore, the mean impact energy of heavy particles and hence the related sputtering is expected to reduce. From energy conservation, the increase of average electron energy flux should be very similar to the decrease in average heavy particles energy flux. When including all effects (see table 4), the electron average energy flux increases by approx. 2 MW m −2 (relative to the case with no effects), so that the total heavy particles energy flux (including fast neutrals) should reduce by a similar amount.
The effects of electron wall emission are clearly visible in figure 11(c), where, on the vertical surface, the net electron flux to the wall (which includes also the negative contribution of emitted electrons) reduces by a large fraction (more than 30% in the case with all active effects). From figure 11(d), it is found that the electron wall emission produces a local raise of the ion flux along the vertical surface of the MB, probably as a consequence of the reduced potential drop in the presheath and Debye sheath. Finally figures 11(e) and (f ) report the mean impact energies of electrons and ions (i.e. referring to impacting particles only). Regarding the electrons, their average impact energy reduces by a factor of 2, along the horizontal surface, when both SEE and TEE effects are included. For the ions, on the other hand, the reduction is around 25% on the horizontal surface. Furthermore, when collisions are considered, ions penetrate farther inside the magnetic shadow region, with energies of approximately 100 eV there.

Conclusions
A two dimensional PIC simulator has been developed to study the thin plasma layer adjacent to the divertor tiles of a fusion reactor. The code has been benchmarked against existing codes in literature and theoretical predictions in unmagnetized and magnetized collisionless cases with no electron wall emission, finding a very good agreement.
The code has been applied to a realistic worst case scenario, considering the ITER divertor geometry with a monoblocks major radius misalignment of 0.3 mm. Considering an artificially increased electron mass (by a factor of 10) to reduce the computational cost, results have shown that, for a magnetic field inclination of 3.2 • relative to the divertor surface, an amplification factor (relative to the energy flux found in undisturbed regions) of around 17 can be observed on the exposed MB surface, confirming previous PIC findings by Komm et al [12]. Moreover, a non-ambipolar profile has been observed along the divertor surface, especially close to the exposed MB edge. This has relevant effects on the local sheath heat transmission coefficient, which presents large fluctuations along the wall and features an average value that is 25% higher than the one predicted by classical sheath theory. Kinetic PIC simulations are therefore key to properly assess this fundamental boundary condition for scrape-off-layer fluid codes.
Finally, the effects of plasma-gas collisions and electron wall emission have been parametrically evaluated in a scenario similar to the above one, but with real ion/electron masses and a reduced radial misalignment of the MBs equal to 0.1 mm. Regarding collisions, it has been found that ion-neutral collisions play a much more important role than electron-neutral collisions, given the smaller (by two orders of magnitude) mean free paths. In particular, both elastic and charge-exchange collisions between hydrogen ions and atoms have the effect of lowering the required ion drift velocity for a stable magnetic pre-sheath and hence the average ion impact energy on the walls, with important implications for ion sputtering. Electron-wall emission, on the other hand, has a strong impact on the total plasma potential drop from magnetic presheath to the wall, which reduces significantly, thus yielding to an increase of the electron energy flux. Due to energy conservation considerations, this means that the heavy particles energy flux (due to both ions and neutrals) must decrease by a similar amount. Summing all effects together (collisions and electron wall emission), it is finally found that the average ion impact energy reduces by 25% and the total heavy particles energy flux by 15%-20%, compared to a case featuring a collisionless plasma with totally absorbing divertor walls. Although these results have been obtained considering hydrogen atoms, it is expected that collisional effects will be even more relevant when considering the more realistic case of deuterium gas, since the corresponding ion/neutral cross sections (for elastic and charge exchange collisions) grow by approximately 30% at the expected center of mass kinetic energies. Moreover, the correspondingly higher ion to electron mass ratio will also induce larger potential drops from plasma bulk to divertor walls.
Future work shall focus mainly on: (i) simulating real divertor MB geometries, including beveling effects to properly protect exposed edges, (ii) performing extended parametric analyses on the effects of different MBs geometries, neutral background densities and wall temperatures, and (iii) including different neutral species (both atoms and molecules), with their relevant collisional processes, as macro-particles of a dedicated Test Particle Monte Carlo (TPMC) module, as done in [9,35]. These particles will be generated at the walls according to known reflection and energy accommodation coefficients [36], and the dedicated TPMC loop will feature a much larger time step (compared to the PIC time step) to account for the much slower neutral dynamics. This latter update will permit to have a fully consistent model of the plasma-wall and the plasma-gas interactions in the divertor region.