Kinetic modeling of the plasma–wall interaction in the DTT divertor region

A precise estimate of the local energy fluxes and erosion profiles at the divertor monoblocks of a fusion reactor requires a kinetic modeling of the plasma–wall interaction. Here, a two-dimensional Particle-in-Cell code is used to quantify the particle and energy fluxes and ion impact distribution functions across the divertor monoblocks of the ‘Divertor Tokamak Test’ reactor, focusing on poloidal gaps with toroidal bevelling. The considered critical locations are close to the strike points at both Inner and Outer Vertical Targets. A worst-case scenario for particle fluxes corresponding to attached plasma conditions and featuring a single-null magnetic configuration is assumed. The separate and cumulative effects of including electron wall emission and ions/electrons collisions with a background neutral gas (recycled at the walls) are also assessed. It is found that a non-negligible energy flux affects the shadowed regions of the monoblocks, especially when accounting for collisions, and that the ion impact distribution functions are strongly influenced by the considered kinetic effects, with important implications on the induced sputtering yield.


Introduction
In state-of-the-art magnetic fusion experiments in tokamaks, the divertor is one of the most critical elements, whose integrity must be guaranteed throughout operations.A significant fraction of the plasma energy and particle flux reaching the divertor region is channelled towards the strike points of the divertor Inner and Outer Vertical Targets (IVT and OVT).This permits to reduce the contamination of the core plasma (compared to the limiter technology) at the cost of a larger peak energy flux to the targets.The Divertor Tokamak Test facility (DTT) [1,2]) is in construction with the specific aim of providing useful information in plasma reactor relevant conditions for the definition of the optimal divertor configuration.The DTT magnetic system is able to realize the standard Single Null Divertor configuration (SND) and all the alternative divertor configurations (ADCs).The DTT divertor has been specifically designed to be optimized for the SN positive triangularity (SND-PT) configuration and to be compatible with the X-divertor (XD) configuration characterized by flux flaring at the targets, and with the SN Negative Triangularity (SND-NT) configuration, which provides a possible solution to avoid ELMs and provide, at the same time, a high energy confinement.To test technological solutions relevant for a fusion reactor the divertor will be realized with actively cooled tungsten monoblocks able to withstand very high stationary heat fluxes (up to 20 MW m −2 ).Although nominally, both DTT and ITER [3,4] will consider, for nominal operations, a partially-detached regime (the best compromise between core plasma performance and sustained divertor fluxes), the worst case in terms of energy fluxes is represented by the attached divertor case.This features a relatively low neutral pressure, an electron temperature of a few tens of eV, and a plasma density of approx.10 20 m −3 [5].Under these conditions, any exposed edge (due to major radius misalignments) to the parallel energy flux would easily reach melting conditions, and this is avoided in DTT (just like in ITER [6]) by shaping toroidally the monoblocks, through both components tilting and bevelling.
In order to assess these effects, simulations dedicated to plasma-wall interaction are crucial.These typically consider as input conditions the plasma properties at critical divertor positions, as obtained from Scrape-Off-Layer (SOL) fluid codes [5], which, in the near totality of cases, assume axisymmetry thus neglecting the monoblocks toroidal shaping effects.A first approach to account for this is to consider optical methods [7], which assume that both ions and electrons move (without excursion) along the magnetic field lines.At an increased level of complexity, the ion-orbit codes [6,7] propagate the ion trajectories considering their finite Larmor radius but assume constant and non-self-consistent electric fields.Finally, the most accurate approach is represented by the fully kinetic Particle-in-Cell (PIC) method [7][8][9], which permits to self-consistently simulate not only the complex divertor geometry, but also the large variety of physical phenomena occurring close to the divertor walls.These include (i) the non-neutral Debye sheath and quasi-neutral magnetic pre-sheath; (ii) the neutrals recycling from the wall; (iii) the secondary electron emission due to both electron and ion impacts; (iv) the thermionic electron emission, (v) the collisions of ions/electrons with neutral atoms/molecules, and (vi) the ion sputtering.Finally, as suggested by previous research [7], PIC models tend to predict a lower amplification factor of the energy fluxes at the exposed surfaces, e.g. at surfaces oriented almost perpendicularly to the magnetic induction field, a worst case scenario that could occur due to major radius misalignments of successive divertor monoblocks along the toroidal direction.
In this work, we present 2D simulation results obtained with the Particle-in-Cell DESPICCO code [9] ('Divertor Edge Simulator of Plasma-wall Interaction with Consistent Collisions'), on the plasma-wall interaction for the DTT divertor.In particular, the critical locations in terms of energy flux at both the IVT and OVT are considered.Both the toroidal bevel angle and the gap between the monoblocks are included.Moreover, by switching on/off additional kinetic physics effects (e.g.collisions between charged particles and neutrals, and electron wall emission from either thermionic or electron-induced mechanisms), the effects of such phenomena on both the energy fluxes and ion impact distribution functions are finally discussed.

The simulation model
DESPICCO, described in detail in [9], has undergone an extensive verification and benchmarking campaign.It consists of a 2D-3V Particle-in-Cell code, which is 2D in position (e.g.covering the toroidal direction y and the direction z normal to an axisymmetric divertor surface) and 3D in velocity in order to correctly capture the charged particles gyromotion and collisions.The model is electrostatic, meaning that only Poisson's equation is self-consistently solved: where ϵ 0 is the vacuum dielectric constant, e is the electron charge, ϕ is the electrostatic potential, n e is the electron density, and n s is the number density of the s th heavy particle species featuring a charge number Z s .In this study, however, only one species of singly charged atomic deuterium ions (Z i = 1) is considered.The number densities of ions/electrons are obtained by weighting the macro-particles to the nodes of a PIC mesh, featuring square cells, with a bilinear interpolation function [9,10].Regarding the boundary conditions for the potential, this is assumed to be fixed at the divertor monoblocks (grounded surfaces), while periodic conditions and homogeneous Neumann conditions (∂ϕ /∂z = 0) are imposed on respectively the lateral (y min , y max ) and upstream (z max ) boundaries.Equation ( 1) is solved with the Petsc solver, using the HYPRE pre-conditioner [11].The magnetic induction field is considered to be uniform and constant in time (no self-magnetic field induced by the plasma current inside the domain).
Charged particles are moved according to the local electric and magnetic fields with the non-relativistic Newton's equation, following the well known leap-frog method, based on the Buneman-Boris explicit algorithm [10].This means that both drifts and gyromotion effects are implicitly included.In order to maintain the requested plasma conditions, both deuterium ions and electrons are injected from a volumetric region adjacent to the upper boundary of the simulation domain, with a fixed influx obtained from SOLEDGE2D-EIRENE fluid code [12].The injection probability distribution function S i,e (v) corresponds to a drifting Maxwellian flux distribution [9]: where m i,e is the elementary ion/electron mass, v is the generated stochastic velocity vector, u i,e is the ion/electron drift velocity, and T i,e is the ion/electron injection temperature.While electrons are injected with no drift velocity, ions are injected with sonic conditions along the magnetic field line (i.e. with conditions corresponding to the magnetic presheath edge).The Maxwellian distribution is well justified for electrons because of their reflection within the sheath, although it may yield to an underestimation of high energy tail effects, as seen in [8].For ions, on the other hand, the Maxwellian-character is not guaranteed by a similar physical reasoning, and other kinetic studies have considered different distributions [13,14].In order to avoid plasma cooling, ions/electrons that cross the upper boundary of the domain are finally re-fluxed, which means that they are re-injected towards the wall with the prescribed injection distribution function S i,e .Periodic reflections along y are finally imposed on particles crossing the y min /y max boundaries along the toroidal direction.
The effects of collisions are simulated by using Monte Carlo Collisions algorithms [9,15], featuring the 'Null Collisions' method.A neutral background of deuterium atoms is considered, fixed in time and uniform in space, in order to reduce the computational cost.In fact, coupling plasma and neutral dynamics would require a large physical time to reach steady-state (given the lower neutrals velocity).A possible solution, which will be considered for future work, would be to separate neutrals and charged species time scales, as done in [16].The effects of the assumed uniform neutral background are expected to be relevant only in regions very close to exposed edges and between monoblocks, and over distances from the wall much larger than the ones considered here.In fact, neutrals are emitted according to the local ion impacting flux (non-uniform along the monoblock) and to the energy reflection coefficient, but are essentially a thermal population, so that their density tends to become uniform quickly.Secondly, neutrals are consumed by ionization collisions as they get farther from the wall, so that significant variations of their density are expected over distances from the wall in the order of their ionization mean free path, which, in the considered scenario, is quite larger than the domain size along z.Regarding the considered collisional processes, electron macro-particles collide with the neutral background through elastic, electronic excitation and ionization collisions, while deuterium ions suffer momentum transfer and charge exchange collisions (cross sections are taken from [17,18]).No Coulomb collisions between charged species are accounted for because these can be considered negligible at the energies of interest (tens of eVs).As an example, at an ion energy of 50 eV, the ion-neutral elastic collision presents a cross section of about 10 −18 m 2 , while the ion-ion Coulomb cross section is 4×10 −20 m 2 .Finally, again for computational reasons, no deuterium molecules are simulated, with their corresponding vibrational and excitation kinetics, which may be relevant very close to the divertor monoblocks, where dissociation has not consumed completely the emitted molecules (see treatment of Molecular Assisted Recombination in a 1D PIC model [19]).
Regarding plasma-wall interaction, every charged particle hitting the divertor walls is eliminated from the simulation, and recorded in terms of particle and energy fluxes deposited on the walls.Additionally, impacting electrons can induce secondary electron emission (SEE) electrons, which are differentiated between back-scattered electrons and true secondaries, following the approach introduced in [20].This model, contrary to some experimental characterizations [21], predicts a zero yield in the limit of zero electron impact energy, and may therefore underestimate the effects of low-energy reflected electrons.No ion-induced SEE is simulated, as this is expected to play a marginal role at the impact energies of interest (hundreds of eVs).In fact, existing experimental data for protons impacting on metals [22] suggests yield values well below the expected ones for electron-induced SEE.Thermionic electron emission is simulated through the Richardson-Dushman's formula [23], which provides the emitted electron current per unit surface as a function of the wall temperature T w and material work function W f : Since the wall temperature is not known a priori, a reasonable engineering value (kept constant throughout the simulation) is assumed.
Finally, volumetric radiation effects are not included, and the plasma is treated as a transparent non-emitting medium with no radiation effects on the walls.This is justified by the scenarios considered (see section 3): a worst-case in terms of particle energy fluxes, for which the relative contribution of volumetric radiation to the total energy fluxes is only around 2%.

Scrape-off-layer simulations
Many different DTT plasma scenarios are presently foreseen both in terms of toroidal magnetic field (B Φ = 3, 6 T), plasma current (I p = 2-5.5 MA) and density (⟨n⟩ = 0.5-2 × 10 20 m −3 ) [24] but also in terms of the divertor magnetic configuration, like the SND in positive or negative triangularity or the ADCs [25], with and without different impurity seeding [25] (like neon or argon), and in terms of divertor conditions spanning from attached, partially-detached to fully detached.For the simulations of this study, an attached plasma scenario has been considered, featuring an SND positive triangularity (see figure 1(a)), I p = 5.5 MA, B Φ = 6 T, an electron density at the separatrix, n e,sep = 0.7 × 10 20 m −3 , a power crossing the separatrix P SOL = 14 MW and without any impurity seeding to enhance the radiation.This plasma condition was selected because it provides one of the worst divertor target conditions in terms of plasma temperature and heat flux but it is still compatible with the actively cooled divertor tungsten monoblocks implementation.Being in pure deuterium, it also simplifies the PIC modeling while providing useful information on possible divertor extreme conditions.
Edge mode has been done with the 2D fluid code SOLEDGE2D [12] using the Monte Carlo code EIRENE [26] for neutrals.Transport values have been estimated for DTT by modeling JET and C-mod experiments and extrapolating  to match the heat flux decay length with the Eich's scaling law [27] applied to DTT, as described in [28].Drifts are not included in the modeling, because they often cause a considerable slow-down of the code convergence, without adding serious physical details, at least at the accuracy level needed for this analysis.The pumping speed has been assumed equal to 50 m 3 s −1 for deuterium, at the 2 pumping slots shown in figure 1(a), in agreement with the present design of the DTT divertor pumping system, featuring a total pumping speed of 100 m 3 s −1 .Deuterium has been injected by gas puffing from the outer mid-plane (OMP), and the particle flux has been adjusted to achieve the requested condition in terms of electron density at the separatrix.The EIRENE code considers the most important neutrals-ions (atomic and molecular) reactions, while the neutral-neutral collisions are neglected, since they are expected to play a minor role in this attached plasma scenario characterized by a low neutrals density.In terms of parameters uncertainty the transport coefficients are those affecting the most the edge modeling results.Indeed a recent numerical study with the XGC1 gyrokinetic code [29] has shown a large deviation of the heat flux decay length from the Eich's scaling for ITER at high plasma current.The newly proposed scaling law predicts a deviation also at the DTT highest plasma current, although much lower than in ITER.This means that the absolute plasma properties at the targets considered here could also be slightly affected, but this does not affect the work performed, whose goal is to describe poorly known physical phenomena at the monoblocks surfaces, rather than to estimate exact properties.

DESPICCO simulation settings
The worst points in terms of total energy flux have been identified on both the IVT and OVT.A zoom on the OVT and IVT locations is shown in figures 1(b) and (c), where the local Cartesian coordinates system used for the DESPICCO simulations is shown in blue (the toroidal direction y enters the page).Note that these worst points do not coincide exactly with the strike points of the magnetic induction field through the null point.The most relevant physical and geometrical parameters for the simulations are summarized in table 1.The higher power and fluxes at the IVT can be attributed to the 1/B Φ dependence in the transport parameters, which narrows the heat flux channel to the IVT on the high field side.This effect disappears in the most useful detached scenario.Finally, both IVT and OVT cases were run considering the real deuterium ion and electron mass (no tricks on the mass to reduce the computational cost), and with 4 different of modeling conditions: no collisions and no electron wall emission ('no effects' scenario), only collisional effects, only electron wall emission effects, both collisional and electron wall emission effects ('all effects' scenario).

Simulation results
The maps of electric potential, ion density and ion Mach number (computed as 2, for both the OVT (top) and IVT (bottom), and accounting for all kinetic effects (e.g.collisions and electron wall emission).The electric potential peaks at around 3-4 T e0 /e, slightly higher than the expected value for a floating wall, due to the braking effect of ion collisions with the neutrals, compensated by a slightly larger potential drop to guarantee the fulfilment of the ion sonic condition at the magnetic presheath entrance [9].Regarding the ions, their upstream density is above 10 20 m −3 in both cases (with the IVT featuring a larger value), as shown in subplots (b),(e).They also tend to penetrate within the magnetic shadow, where, however, their density is below 10 18 m −3 .Regarding the ion fluid velocity, shown in subplots (c),(f), this aligns with the magnetic field direction, reaching sonic conditions at least 0.5 mm from the divertor surface.The monoblock edge finally distorts the ion streamlines, which bend towards the material surface to satisfy the classical Bohm's condition (perpendicular Mach number greater or equal to 1).
Figure 3 then shows the profile of the energy fluxes along the divertor wall, for both the OVT (top) and IVT (bottom) critical points, with a zoom on the exposed monoblock edge (b),(d).The curvilinear coordinate s is 0 at the left-most simulation domain node (at the plasma-wall interface) and grows from left to right along the monoblock profile covering also the vertical shadowed surfaces (i.e.along z) and the gap region: after the left monoblock edge, it grows downwards along its vertical surface, then along the z min boundary, and upwards along the vertical surface of the right monoblock.
While the electron energy flux (blue lines) is almost constant along the exposed monoblock surface for both OVT and IVT and in all modeling scenarios, meaning that the optical approximation is relatively good, the ion energy flux (red lines) tends to grow closer to the exposed edge because of the finite ion gyro-radius effect and, more importantly, because of collisions (dashed and dash-dot lines).In fact, ions reaching the monoblocks farther away from the edge travel a larger distance and feature a higher probability of suffering a collision with a neutral atom.This produces a slow ion, and the corresponding energy carried to the walls is therefore lower.The overall effect of including collisions is to reduce the average ion energy flux by up to 10 MW m −2 .However, this reduction is partially compensated by the neutrals energy flux, which has not been computed here (neutrals are not followed kinetically).Although future work is required, it is reasonable to expect a neutral energy flux lower than this energy flux reduction (10 MW m −2 ), because neutrals are not accelerated inside the Debye sheath.This also yields a lower impact energy and a significantly lower sputtering induced by neutrals.
Finally, from figures 3(b) and (d), it can be appreciated that collisional effects produce a significant raise in the ion energy flux on the shadowed vertical surface of the monoblock.In the OVT scenario, this increases from fractions of MW m −2 to even more than 5 MW m −2 just around the edge (on the exposed surface it is nearly 20 MW m −2 ).Non-zero electron energy fluxes can also be observed when accounting for electron wall emission (dotted and dashed lines), and they are due to emitted thermionic electrons from other shadowed monoblock surfaces (see figure 2).Overall, SEE electrons have a relatively weak effect on the total energy fluxes, as shown in (a) and (c), and this may be explained by the prompt electron re-deposition in the presence of a small angle between the magnetic field and the monoblock surface [30].
In order to correctly model the ion-induced sputtering, the ion impact distribution functions have finally been computed, as a function of both the ion impact energy E imp and the ion impact angle θ imp .Figure 4 shows the distribution for both the OVT (top) and IVT (bottom), at two different locations: the first (location 1) in the magnetic shadow and the second (location 2) on the exposed surface, both next to the exposed edge (magenta circle marker in figure 2).
Results are also compared among the different modeling scenarios, with solid lines referring to the 'no effects' case, and dashed lines to the 'all effects' case.Non-negligible ion fluxes are observed for impact energies larger than 50 eV (due to the acceleration within the Debye sheath) and angles ranging from a few degrees (nearly normal impact) to more than 80 deg (grazing incidence).At the shadowed locations (subplots (a),(c)), the inclusion of collisions substantially increases the ion flux.At the exposed location (subplots (b),(d)), including collisions broadens significantly the distribution in both energy and angle: energies of up to 500 eV are observed with angles ranging from a minimum of around 10 degrees to a maximum close to 80 degrees.In all cases, the effects of collisions is that of lowering both the average impact angle and energy, as shown by the circle markers.

Conclusions
An innovative study for a single null attached plasma scenario of the DTT experiment has been carried out with the DESPICCO code, for two critical locations at the Outer and Inner Vertical targets.Simulations have considered the real geometry of the DTT monoblocks, accounting for the bevelling effect along the toroidal direction.The effects of collisions and electron wall emission physics have been evaluated in detail.
Results show that ion-neutrals collisions can play an important role in determining both the peak value and the wall profile of the energy flux to the monoblocks.In particular, a peaking effect at the exposed edge is observed for the ions.Moreover, non-zero energy fluxes are also found at the shadowed monoblock surfaces close to the edge, with values of several MW m −2 .Finally, the ion impact distribution function around this exposed edge have been computed, finding that collisions broaden significantly the impact energy and angle ranges, while lowering their average values.
Future work shall focus on (i) assessing the sputtering yield corresponding to the obtained ion impact distribution functions, (ii) simulating the effects of additional species (e.g.due to gas/impurities seeding), (iii) analysing a different critical location of the divertor (e.g. the dome, for an XD magnetic field configuration), (iv) comparing PIC results with a faster ion orbit code, and (v) performing coupled TPMC/PIC simulations to capture self-consistently the neutrals dynamics and the plasma-gas coupling [16,19,[31][32][33], thus avoiding to assume a fixed and uniform neutrals background.This last topic will also allow to estimate precisely the contribution of neutrals to the total energy fluxes to the walls, that has not been estimated in this study.Finally, it is of great importance to carry out experiments on the characterization of the plasma-wall interaction at the monoblocks (estimation of sputtering profiles, heat fluxes, floating potentials, etc), as this would enable a validation of the numerical tools, thus increasing the confidence level in their predictions.

Data availability statement
The data cannot be made publicly available upon publication because they are not available in a format that is sufficiently accessible or reusable by other researchers.The data that support the findings of this study are available upon reasonable request from the authors.

Figure 1 .
Figure 1.(a) Poloidal view of the DTT first wall and vessel outer wall, showing the considered SND magnetic configuration and the SOLEDGE2D simulation domain (light blue region).The location of the deuterium puffing point and of the pumps entrance is also shown by respectively a red dot and thick blue lines.Zooms on the critical (b) OVT and (c) IVT points considered in PIC simulations, where the local DESPICCO simulation axes x, y, z are shown in blue, and the magnetic field direction in the poloidal plane by a red arrow (the dominant component is perpendicular to the page, along the toroidal y direction).White arrows indicate the wall normal direction.

Figure 2 .
Figure 2. Results for the (a)-(c) OVT and (d)-(f) IVT critical points, under the 'all effects' modeling scenario, showing (a), (d) electric potential, (b), (e) ion number density, and (c), (f) ion Mach number.The divertor monoblocks are shown in grey, the magnetic field direction with a white arrow, while on subplots (c), (f), the black lines with an arrow refer to the ion streamlines.The exposed monoblock edge is shown by a filled magenta circle marker.The gap considered here is a poloidal gap, i.e. belonging to the poloidal plane (normal to the considered simulation domain).

Figure 3 .
Figure 3. (a) Energy fluxes along the divertor monoblock curvilinear coordinate s at the (a), (b) OVT and (c), (d) IVT critical locations.Magenta dashed lines show the location of the exposed monoblock edge (shown by a magenta circle in figure 2).Subplots (b), (d) show a zoom around this edge.Black, blue and red lines correspond respectively to total, electron and ion energy fluxes.Results are shown for the 'no-effects' (solid), 'collisions only' (dash-dot), 'wall emission only' (dotted), and 'all effects' (dashed) modeling scenarios.In subplots (b), (d), a logarithmic scale is considered, and the total contributions are omitted for the sake of clarity.

Figure 4 .
Figure 4. Ion impact distribution functions at (a), (c) shadowed and (b), (d) exposed locations next to the monoblock edge, for the (solid) 'no effects', and (dashed) 'all effects' scenarios.Triangle/circle markers refer to the distribution peak/average values, for the 'no effects' (filled markers) and 'all effects' (empty markers) modeling scenarios.θ imp = 0 deg corresponds to normal to surface incidence.The integral of the distributions over the energy and angle 2D space gives the total impacting ion flux.Coloured contour maps of the distributions are available as supplementary material.

Table 1 .
Simulation parameters for the critical points on the IVT and OVT.