Particle-in-cell simulations of the plasma interaction with poloidal gaps in the ITER divertor outer vertical target

Predictive modelling of the heat flux distribution on ITER tungsten divertor monoblocks is a critical input to the design choice for component front surface shaping and for the understanding of power loading in the case of small-scale exposed edges. This paper presents results of particle-in-cell (PIC) simulations of plasma interaction in the vicinity of poloidal gaps between monoblocks in the high heat flux areas of the ITER outer vertical target. The main objective of the simulations is to assess the role of local electric fields which are accounted for in a related study using the ion orbit approach including only the Lorentz force (Gunn et al 2017 Nucl. Fusion 57 046025). Results of the PIC simulations demonstrate that even if in some cases the electric field plays a distinct role in determining the precise heat flux distribution, when heat diffusion into the bulk material is taken into account, the thermal responses calculated using the PIC or ion orbit approaches are very similar. This is a consequence of the small spatial scales over which the ion orbits distribute the power. The key result of this study is that the computationally much less intensive ion orbit approximation can be used with confidence in monoblock shaping design studies, thus validating the approach used in Gunn et al (2017 Nucl. Fusion 57 046025).


Introduction
The question of monoblock (MB) top surface shaping is the key outstanding design issue for the ITER tungsten (W) divertor and is an important feature for any discrete comp onent (i.e. consisting of edges and gaps) located in regions of high heat flux and glancing magnetic field incidence [2]. The shaping should ensure protection of leading edges (LE) but not at the expense of too much loss of steady state power handling capacity.
An extensive shaping study has been performed for the ITER divertor using optical and 3D ion orbit approx imations to estimate the heat flux distribution along the poloidal and toroidal gaps between MBs [1]. Although these calculations do capture most of the zero order power loading effects, including in particular those governed by the finite ion Larmor radius (r Li ), they do not account for the role of local electric fields which arise in the Debye sheath and magnetic presheath, and which can have a significant impact on particle trajecto ries. In order to assess this effect, a more sophisticated model Predictive modelling of the heat flux distribution on ITER tungsten divertor monoblocks is a critical input to the design choice for component front surface shaping and for the understanding of power loading in the case of smallscale exposed edges. This paper presents results of particleincell (PIC) simulations of plasma interaction in the vicinity of poloidal gaps between monoblocks in the high heat flux areas of the ITER outer vertical target. The main objective of the simulations is to assess the role of local electric fields which are accounted for in a related study using the ion orbit approach including only the Lorentz force (Gunn et al 2017 Nucl. Fusion 57 046025). Results of the PIC simulations demonstrate that even if in some cases the electric field plays a distinct role in determining the precise heat flux distribution, when heat diffusion into the bulk material is taken into account, the thermal responses calculated using the PIC or ion orbit approaches are very similar. This is a consequence of the small spatial scales over which the ion orbits distribute the power. The key result of this study is that the computationally much less intensive ion orbit approximation can be used with confidence in monoblock shaping design studies, thus validating the approach used in Gunn et al (2017 Nucl. Fusion 57 046025).
Keywords: tokamak, plasma, ITER, particleincell, heat loads, monoblock (Some figures may appear in colour only in the online journal) Original content from this work may be used under the terms of the Creative Commons Attribution 3.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. based on the particleincell (PIC) technique is required. The purpose of this paper is to report on an effort to benchmark selected cases from the ion orbit calculations used in [1] with the much more computationally demanding PIC approach and to examine whether or not the simpler approach is adequate as a physics design tool for shaping optimization. The emphasis will be on the case of divertor plasma conditions in between edge localized mode (ELM) transients.
The plasma interaction with ITER MBs has already been a subject of a number of PIC simulation studies [3][4][5]. The results presented here are, however, novel in two main aspects: 1. The simulated plasma parameters target the baseline burning plasma divertor scenario, generated using the SOLPS4.3 plasma boundary code suite [6] for inter ELM conditions. This had not previously been attempted due to the extreme computational demands. 2. The geometry and magnetic field orientation are those appropriate to the ITER divertor design.

Geometry of the simulations
The ITER divertor targets will be castellated, i.e. split into individual MBs separated by gaps. The nominal MB in the high heat flux (strike point) areas of the vertical targets extends 28 mm toroidally and 12 mm poloidally [1,7]. The MBs are bonded to cooling pipes extending in the poloidal direction, constituting a 'plasmafacing unit' (PFU). Several PFUs comprise a single vertical target, which is in fact a com bination of two 'halves' which are separated to reduce eddy current forces during disruptions. A similar design is used for both the inner vertical target (IVT) and the outer vertical target (OVT). Nominal gaps of 0.4 mm separate poloidally adjacent MBs in the HHF areas (thus between MBs on the same PFU, referred to as 'toroidal gaps' (TG)). Between toroidally adja cent MBs (thus on separate PFUs and referred to as 'poloidal gaps' (PG)), nominal gaps are 0.5 mm. Wider PGs separate the two 'halves' of each vertical target (3 mm) and are even wider between individual divertor units (20 mm) of the 54 divertor 'cassette' assembly. Each vertical target is globally tilted toroidally by ~0.5° to protect intercassette LEs. This paper will focus on the small interPFU PG between toroi dally adjacent PFUs in one half of a vertical target. The ITER design defines allowable tolerances for the divertor assembly. One of the most important with respect to MB heat flux deposition is the relative radial misalignment of toroidally adjacent MBs on neighboring PFUs. The largest allowable misalignment which the manufacturers of the ver tical targets must respect is specified as m rad = 0.3 mm (see table 1 in [1]). In order to simulate the worstcase scenarios, this misalignment has been used in all simulations of PGs reported here.
The local coordinate systems and magnetic field orienta tions are shown in figures 2 and 3 in [1]. The magnetic field is assumed to be constant within the area of interest.

Description of the SPICE2 code
The electric field in the electrostatic magnetic presheaths, as well as the electric potential in the vicinity of the gaps in between MBs, can affect the plasma deposition in such geom etries. The PIC technique selfconsistently generates these regions, describing fully the resulting trajectories of charged particles. SPICE2 (Sheath ParticleInCell) is a 2D3V (2D in space and 3D in velocity) particleincell code developed over several years in collaboration between IPP Prague and CEA Cadarache. It was designed to perform simulations of magnetized plasmas in contact with solid objects and has been successfully used for the study of plasma deposition in the vicinity of castellated plasmafacing components (PFCs) [3,8,9] and probes [10]. SPICE2 calculates the motion of charged particles in a prescribed static magnetic field and a selfconsistent electric field, which is calculated using a direct Poisson solver based on the LU decomposition [11] (using the UMFPACK library [12]). The code is optimized for simula tions of particles in the electrostatic sheath and magnetic pre sheath on PFCs oriented at oblique angles with respect to the magnetic field. A typical geometry for the SPICE2 simulation box is shown in figure 1. Particles are injected from a box located at the upper boundary, while the limiter/divertor tiles are located at the lower boundary of the simulated region. For grazing angle of incidence, a magnetic sheath can develop from the surface of the object towards the plasma with a typical scalelength √ 6c s /Ω i [13], with c s = (kT e + kT i )/m i being the ion sound speed. This is accounted for here by ensuring a min imum distance of 10c s /Ω i between the top of the tile and the upper simulation boundary, thus avoiding any perturbation to the bulk plasma.
We use the result of a collisionless kinetic model of the scrapeoff layer (SOL), which calculates a sheathedge ion parallel velocity distribution function satisfying the Bohm criterion [14]. Ion perpendicular velocity components and all electron velocity components are assumed to be Maxwellian. The yaxis corresponds to the toroidal direction in case of PG simulations and the zaxis is aligned to the target surface normal vector. In the case of castellated PFCs, SPICE2 can implement a variable tile geometry, as long as it can be approx imated by rectangular, triangular or rounded shapes. It can be used to simulate shaped, unshaped or misaligned MBs and the 3D magnetic field orientation can be chosen to simulate PGs or TGs. The tiles are usually set at a prescribed potential, with respect to the plasma potential, which is fixed at the upper boundary. The side boundary conditions are periodic.
SPICE2 simulates collisionless plasma only and all the sim ulated objects are perfectly absorbing, meaning that particles reaching material surfaces are removed from the simulation. The simulation starts with an empty region, which is gradually populated by injected particles. Once the simulation reaches a steadystate level (i.e. the total number of particles, or density, is constant), dedicated diagnostics acquire desired outputs, such as particle and heat fluxes at the surfaces, spatial distributions of electric fields, densities, velocities, etc. Most of these diag nostics involve a large number of temporal or spatial samples (~10 4 ) averaged together. In the work reported here, the time step of the simulations was determined by the electron thermal velocity (electrons must not cross more than one cell during one time step), so an artificially reduced ratio of ion and electron mass = 200 was used to reduce the electron velocities, allowing an increase in time step and a reduce computational time for the simulations. In order to evaluate the effect of reduced mass ratio on the sheath properties, a series of simulations with plasma parameters corresponding to cases PG.1, PG.2 and PG.3 (see table 1) with µ = 200 and 4650 was performed for simple planar surface. The profiles of electric field E z acting normal to surface are plotted in figure 2. The profiles show very good agreement for z > 40 µm (z = 0 mm corresponds to surface location), which corresponds to half of ion Larmor radius. We can therefore conclude that the modification of electric field profile due to use of reduced mass ratio should not affect sig nificantly the ion deposition profiles on MB surface.

Boundary conditions for simulations with radial misalignment
Although, as mentioned earlier, the simulations take into account the radial misalignment between toroidally adjacent MBs (m rad ), this is not compatible with the periodic boundary conditions (BCs) imposed on the sides of the simulated region. As shown in figure 3, this problem can be solved by extending the simu lated region to the right in such a way that the misaligned tile is contained together with a region where the particle fluxes are perturbed by the magnetic shadow and the modified electrostatic potential generated by the preceding misaligned tile.
This approach leads to correct results, but due to small angles of incidence in the cases under study, the extension of the simulated region must be rather significant, increasing the computational demands of such simulations. An alternative solution consists of a modification of the boundary conditions such that whenever a particle passes through one of the side boundaries, it is displaced vertically by ±m rad (figure 4). In this way, no magnetic shadow is cast and the extension of the simulated region can be significantly reduced. However, to obtain the correct potential profile, the Poisson solver must also be modified so that stencils which extend through the periodic boundary are composed from the grid points which are also vertically displaced. This leads to a strip of grid points near the injection region where the potential is not properly defined. Since the potential profile in this region should be flat and close to zero (equal to the plasma potential), the poten tial on such grid points is fixed at zero. Even if it is observed though in some of the simulations that this strip leads to per turbation of the potential, the impact on the tile heat flux dis tribution is not significant.

Electron heat flux
SPICE calculates the motion of both ions and electrons and therefore it is in principle capable of providing the heat flux delivered to the tile by both species. However, in the case of computationally demanding simulations, where a reduced iontoelectron mass ratio is used (µ = 200), the electron heat flux is artificially increased due to finite Larmor effects [15]. Moreover, for large simulated regions, numerical fluctua tions of the potential can induce artificial heating of electrons, which can be of the order of 20% in the most demanding cases. At the same time, the electric fields in the sheath are insufficiently strong to cause demagnetization of electrons and so they follow the Larmor motion. For these reasons, in the simulations presented here the optical approximation is used to compute the electron heat flux. When calculating the electron heat flux magnitude, standard sheath theory is used and the following formula [16] is adopted: where q i is the ion heat flux received by the planar top sur face, τ = T i /T e and V s is the potential drop in the sheath nor malized to T e (V s = 3 in all simulations).

Local coordinate system
In order to systematically represent the heat fluxes received by the MB surfaces, the s tor coordinate adopted in [1] is also introduced here (see figure 13 in [1] and figure 5), with s tor = 0 located at the corner of the plasma wetted gap side and being positive on the MB top surface.

Sensitivity to the shape of the ion velocity distribution function
In order to assess the differences in power deposition profiles caused by the presence of the local electric fields, the first step is to ensure that when electric fields are not taken into account in SPICE2, the results of the simpler ion orbit approach are recovered for the same geometry and plasma conditions. Even in absence of Efields, the PIC and ion orbit codes differ in the choice of the ion parallel velocity distribution function (DF). SPICE2 uses as input a DF calculated using the Chung and Hutchinson (C&H) model [14]. Since the ion orbit code used in [1] does not include ion acceleration in the sheath, simulating particles with parallel velocity close to zero (which are included in the kinetic distribution) can require long simulation times. To eliminate these low velocity particles, the ion orbit code normally employs a simple top hat DF (see figure 11 in [1]) having the same mean velocity and thermal width as that resulting from the C&H model. Moreover, since the ion orbit model is based on an analytic solution of the equations of motion, which are helices, it is not subject to cumulative gyrophase error, which arises in the PIC method due to discrete time steps [17].
To benchmark the PIC and ion orbit approaches, a PG PIC simulation has been performed with the Efield turned off (referred to henceforth as a 'ballistic' simulation) and the resulting ion heat flux compared to that obtained from a 2D ion orbit code run for the same geometry and plasma   conditions (corresponding to case PG.1 in table 1). Both codes used an identical DF generated from the C&H model. The results, shown in figure 6, demonstrate very good agree ment between the ion orbit and PIC simulations. Note that the heat fluxes are normalized to q top , which corresponds to the heat flux received by the plasmafacing surface of a MB in the absence of any tilting or surface shaping (perfectly axisymmetric surface).
It should be noted that the C&H distribution itself, although accurately describing the sheath edge distribution function on the basis of a first principles physical model, corresponds to a collisionless, fully ionized plasma. At the ITER divertor tar gets, the interELM plasma in high performance phases, of most interest here, will generally be highly collisional due to the elevated plasma density. The motion of particles within the sheath under these conditions is likely to be modified by ion-ion collisions, which can occur at frequencies exceeding the ion Larmor frequency. In addition, the real distribution function at the sheath edge should ideally be calculated by a fully collisional kinetic model (such as the BIT1 code [18]), including interactions with neutrals and excited impurity spe cies. In this sense, the SPICE2 calculations described here represent only a first order approximation to the complete solution of the problem of plasma interaction with MBs under ITERrelevant conditions.

Plasma scenario specifications
The key plasma scenarios of interest are fully described in [1] and consist of the burning plasma, steady state, dissipa tive (partially detached) baseline case and the 'slow transient' case, specified to correspond to roughly twice the baseline heat flux density and intended to capture the condition in which there is a temporary loss of partial detachment control. The timescales (steady state, or slow transient on the order of seconds) have no meaning in the PIC simulations since only times shorter than 1 µs can be followed.
The ion orbit simulations [1] show that in steady state (when divertor plasma temperatures are very low), r Li m rad and as such finite orbit effects cannot provide significant power spreading over an ITERrelevant LE (such spreading is only observed for m rad < 2 r Li at almost normal incidence to the edge [19]). It is possible, however, that modifications of the particle orbits due to local electric fields may play a role in reducing the heat flux on such misaligned edges. This is the principal motivation for the application of a more realistic PIC approach.
The plasma parameters required for the PIC simulations are identical to those assumed in [1] for the ion orbit and optical approximation calculations and extracted from SOLPS 4.3 boundary code simulations of the ITER W divertor. All three simulation approaches (PIC, ion orbit and optical) thus employ the same plasma background and can thus be reliably compared. Calculations are restricted here to the OVT since the steady state heat flux densities are predicted to be highest there and the angles of incidence lowest, giving the most con servative situation for misaligned edge loading.
As already mentioned, the computational demands of the PIC simulations restrict the size of the simulated area and the number of time steps. Using the geometry shown in figure 3, the width of the box must be sufficient to compensate for magnetic shadowing. The new scheme employing a boundary condition with a jump (figure 4) does not require this and so the simulation box can be shorter, however at the risk of intro ducing a perturbation in the potential field. To examine this issue, some of the computationally less demanding cases (see table 2 for the list of all PG simulations attempted) were run with both schemes. The more computationally difficult cases are calculated only with the new periodic boundary condition to render the problem more tractable.
The dimensions of the simulation (without the new boundary condition scheme) have already been noted in figure 3. Apart from the length covering the magnetic shadow (y 1 = m rad /tanθ ⊥ ), an additional length corresponding to at least 2 Larmor gyrations along the field (Larmor helix length) should be allowed on each side of the gap (y 2 = 4π c s /ω ci , with ω ci the Larmor frequency). The height of the simulation box above the misaligned tile should be sufficient to allow the magnetic presheath to develop naturally in the simulation and should include a gap depth z gap of at least 2r Li below the optical wetted length on the gap side. In the case of the new boundary condition scheme, the length y 1 is set to zero and the box height is incremented by m rad . The local heat flux calcu lated by each model is normalized to the the nominal heat flux that would strike an ideal, axisymmetric divertor surface for a given launched parallel heat flux, q ⊥,pk = q sinθ, where θ is the total angle of incidence with respect to the axisymmetric divertor surface.
Unfortunately, since PIC code cell sizes must be at most equal to the Debye length λ D = (ε 0 T e /n e e) 1/2 to guarantee numerical stability [17], simulations which simultane ously respect the SOLPS4.3 generated background plasma parameters and satisfy the geometrical constraints described in figure 3 would require enormous simulation boxes. This would render simulations with any reasonable particle num bers impracticably long. In particular, the high divertor densi ties cannot be matched.
Since the heat flux density and the scale size of the mis alignment in comparison with r Li are the primary quantities of interest, the two key parameters to be preserved in the simula tions are r Li and q ⊥ ,pk . This can be achieved by increasing T e (thus decreasing τ) and decreasing n e , which will increase λ D and thus decrease the number of cells required in the simu lation box, while allowing the heat flux to remain constant. However, this procedure is not without penalty, since λ D is modified which in turn modifies the electric field and the heat flux distribution. The assumption here is that as long as the ordering λ D r Li m rad is preserved, the modification of the heat flux distribution remains sufficiently small. In order to verify this hypothesis, a series of cases have been run in which the real SOLPS4.3 divertor plasma parameters are gradually approached (i.e. PG.1 to PG.4 in table 1) and the differences in the heat flux distributions studied.
Since the SOLPS4.3 simulations take no account of target tilting or MB surface shaping, the q ⊥ ,pk values are for an axisymmetric surface only. In fact, the SOLPS4.3 code output is in terms of T e , T i , n e , etc in front of the target plates, from which the power density is derived through simple sheath theory (i.e. the application of a sheath heat transmission factor). Thus, for example, in the steady state (interELM), dissipative (partially detached) case (PG.4 in table 1), SOLPS 4.3 yields q ⊥ ,pk = 9.3 MW m −2 at the OVT, corresponding to q || = q ⊥ ,pk /sinθ = 197 MW m −2 for a total angle of incidence θ = 2.7°. The plasma parameters yielding this peak power density are n e = 5.8 · 10 20 m −3 , T e = 12 eV, T i = 4.7 (τ ~ 0.4). The local magnetic field is B = 6.0 T. In reality, global target tilting increases this angle to 3.2°, in turn increasing the top surface heat flux density to ~11 MW m −2 .
With these points in mind, a series of PIC simulations for a misaligned PG separating unshaped monoblocks has been performed, beginning with an artificially low value of τ and increasing the latter towards more realistic values obtained from the SOLPS4.3 simulation (τ ~ 0.4) for both the steady state and 'slow transient' conditions. The parameters of the simulations are summarized in table 1. Each has a corre sponding ion orbit code calculation with which the results can be compared. Cases PG.1-4 represent a series of runs approaching the target conditions for steadystate (PG.4), while PG.5-6 target the slowtransient conditions (PG.6).
In order to illustrate the computational demands of the simulations, table 2 contains grid sizes and run times for each simulation which was performed. The simulations were com puted at the HELIOS supercomputer on 1632 cores depending on the grid size. Suffix S in tables 1 and 2 denotes simulations in which the periodic boundary condition scheme with ver tical jump is used ( figure 4). All other cases use the traditional BC scheme. Note that the computational demands depend not only on the grid size but also on how the simulations could be distributed between computational nodes (depending mainly on the RAM requirements) and on the chosen Poisson solver. For case PG.4S, the grid size was already too large to be han dled by the UMFPACK Poisson solver employed as standard in SPICE2. Instead, the distributed, but slower solver using the MUMPS library had to be employed, further increasing the run time. It can be seen that use of the BC scheme with jump reduces the computational demands by a factor of 2 in case of run PG.1, which is already significant, given the time scale of simulations. Figure 7 compiles results for case PG.1 obtained from 4 models: optical approximation (no ion orbits), ion orbit model [1], 'ballistic' PIC code (Efield turned off) and PIC code (Efield switched on). The plasma wetted area of the PG and the misaligned edge is exposed to nearly perpendicular heat flux, reaching more than 20× the optical flux impacting the top surface of the MB. The ion orbit and ballistic PIC models show only a minor reduction of this flux, which is due to the very small ion Larmor radius compared to m rad (r Li = 80 µm = 0.25 m rad ). The ion and electron heat fluxes to the top surface (far away from the gap) are quite similar in the full PIC simulation results, in contrast to the expectation from simple sheath theory of q i /q e = 5/2 [16]. This is due to the fact that T i /T e = 0.1, so the energy of ions hitting the surface is practically only due to energy gained in the sheath, making the ratio q i /q e = 3/2 (the ratio 5/2 is valid for T i /T e = 1).

Simulated heat flux profiles
At the top of the LE the flux is determined by the unper turbed parallel flux, whereas with increasing depth into the  2 3 × 10 6 -PG. 2S 1.3 × 10 6 64 PG. 3 6.5 × 10 6 -PG. 3S 2.4 × 10 6 138 PG. 4 10.6 × 10 6 -PG. 4S 3.7 × 10 6 300 PG. 5 1 × 10 6 45 PG. 5S 7 × 10 5 -PG. 6 3.4 × 10 6 -PG. 6S 1.8 × 10 6 100 gap, the flux diminishes due to ion orbits being scraped off by the preceding MB further upstream. In previous simula tions of PGs between misaligned tiles [19], it was observed that due to the finite Larmor radius, the LE in some cases can be protected from the ion component of the heat flux. This is, however, only significant when the misalignment is less than ~2r Li . In figure 7, the LE power loading due to ions is lower on the gap edge in the full PIC simulation by nearly 20%. This difference is due to the sheath and presheath electric field, as will be explained below. Energy conservation dictates that the total power injected into the simulation domain must be recovered on the surface, whatever its shape, and whatever model is employed. That is, the integral of the surface heat load should be identical for all cases. The reduced power load on the exposed edge seen in figure 7 is compensated by a zone of elevated heat flux onto the top surface close to the gap. This enhancement with respect to an infinite flat surface is due to ion orbits which extend below the surface into the volume preceding the LE before passing over it on their last gyration; such orbits are unpopulated in the case of an infinite surface. This qualita tive explanation applies to both ion orbit and PIC models. The difference between the spatial distributions is an effect of the electric field, as will be discussed below. Integrating the total heat flux received by the elevated edge (between s tor = −0.5 and s tor = 2.3 mm), the heat fluxes received by the LE in the different models agree within ~3%. These small discrepancies are within the error of numerical fluctuations.
Simulations of steady state plasma conditions are the most difficult cases to calculate due to the very short Debye length. The target case PG.4 (table 1) has λ D = 1 µm, so that in order to cover the entire MB surface in the toroidal direction (28 mm), N y = 28 000 cells on the PIC grid would be required in the y direction alone. This is many orders of magnitude beyond the capabilities of the SPICE2 code. Even when the simulated region is restricted to a smaller area around the PG and the periodic boundary condition with a jump is applied, the estimated runtime for PG.4S is ~300 d (as shown in table 2), which is beyond the computational resources avail able for a single run.
In order to assess the heat flux deposition for the steady state case, the series of proxy cases PG.1, PG.2S and PG.3S, as shown in figure 8, are examined to determine if the profiles converge. Evidently, there is no significant variation between the cases, suggesting that the mechanism of ion transport to the LE is dominated by gyromotion perpendicular to the field lines since r Li is maintained constant for all simulations (by means of adjusted T i /T e ). The deposition on the top surface is located slightly closer to the MB edge for cases PG.2S and PG.3S, which is expected due to the smaller specified  T e . Significant differences are therefore not expected for case PG.4S, which closely approaches the real SOLPS4.3 param eter range. Note the perturbation of the heat flux profile in case PG.3S at s tor = −0.25 mm. This is an undesired result of the potential perturbation caused by the BC scheme with ver tical jump. However, the magnitude of the effect is sufficiently small for it to be considered insignificant for the analysis.
For the case of slow transient, or more attached divertor con ditions, two PG cases have been attempted (table 1), one (PG.6S) which matches the SOLPS value of τ (=0.2) and a second (PG.5) which is a computationally less demanding proxy, following the same philosophy as for cases of PG.1-3 for the more detached plasma conditions. The results of the simulations are shown in figure 9 and are qualitatively similar to those for the steadystate conditions. The perturbation caused by the BC scheme is again visible on the ion heat flux (red curve) deposited on the exposed edge of the MB (s tor < 0 mm). The deposition on the top surface is limited to ~1.5 mm in the vicinity of the MB edge, as for the PG.13 simulations in figure 8.

Role of electric fields
2.9.1. Potential structure in the gap. An effect which origi nates from the presence of electric fields in the sheath is the formation of a potential structure inside the gap, which under certain conditions affects the ion motion in the gap. This positive potential structure is caused by ions, which, unlike electrons, can enter the magnetic shadow inside the gap due to their finite Larmor radius. Previous simulations [9] have shown that when the potential inside the structure exceeds the plasma potential, it can prevent ions from entering the gap and force some of them to land on the top surface of the neighbor ing MB, contributing to the peaking of heat flux there. In the simulations for ITER conditions described here, this effect is not observed for two reasons: 1. The potential structure in all cases remains well below the plasma potential (as shown in figure 10 for selected cases). 2. Since m rad ~ 4r Li , the ions influenced by the potential structure cannot reach the top surface of the LE.

The role of polarization drift in the magnetic pre-
sheath. The particle flux in the Debye sheath and magnetic presheath is modified by the polarization drift [20, 21]. Con sider a simple planar surface located at z = 0 with grazing incidence of the magnetic field and in the absence of an elec tric field. As the ions approach the surface, they are gradu ally scraped off from the plasma as a result of their Larmor motion. This effect creates a local decrease of ion density which is the principle mechanism behind the formation of the magnetic presheath. It also gives rise to an apparent velocity component pointing towards the surface, simply due to the fact that in this region some particles moving out of the sur face (due to their Larmor rotation) are missing, because they have been previously scraped off by the surface.  When the sheath electric field is present, the ions passing through the sheath experience a temporally evolving local Efield in their frame of reference, which drives a polariza tion drift. This drift is a real guiding center velocity pointing towards the surface. Due to the strong Efield gradient above the surface, the polarization drift velocity can easily exceed the normal projection of the parallel velocity of the ion, leading to substantial distortion of its trajectory. The processes can be best illustrated by comparing ion streamlines (contours tan gential to average velocity vectors at every grid point of the PIC grid) in all 3 cases, as shown in figure 11.
The fundamental mechanism of LE heat loading for both ion orbit and PIC models, as well as the differences between the two, can be understood with the simplified geometry of a step on the surface, shown schematically in figure 12. The gap is not considered since it does not influence the heat loading at the top of the LE in these cases.
In the case of the ion orbit model, the particle deposition profile at the LE is identical to the parallel particle flux Γ || profile above an infinite planar surface at any position fur ther upstream. This is a trivial result which arises because the magn etic field lines are nearly parallel to the surface, and the ions have no interaction with the LE until the instant they strike it. The LE flux profile is determined uniquely by the scraping off of ions on the preceding surface, without any modification of the velocity of those reaching the LE.
In the case of the PIC model including interaction with the surface via the selfconsistent electric field, while the poten tial profile above a planar surface scales approximately as ~exp(z/r Li ), the potential in front of the leading edge scales only as ~exp(y/λ D ), since the LE is almost perpendicular to the magnetic field (note that in all studied cases, r Li is more than an order of magnitude larger than λ D ). As a consequence of the supersonic motion due to acceleration in the magnetic presheath, particles approaching the LE do not know what is ahead of them until they reach the LE Debye sheath. However, this zone is so narrow that it does not allow for substantial modification of the flow pattern. Since information about the presence of the LE does not propagate upstream against the supersonic inflow, similar reasoning as for the ion orbit model applies, and the deposition profile on the LE is very similar to the parallel flux profile calculated by a 1Dlike PIC simulation representing flow to an infinite surface. Figure 13 gives a com parison of the PIC simulated ion flux profiles on the LE for the case PG. 1 with Efields on (full PIC) and off ('ballistic' case). In both cases, the LE flux profiles calculated on a 2D grid are nearly identical to the result of a 1D calculation. The PIC heat flux is reduced compared to the ballistic case because the energy gained in the sheath is shared between the parallel and perpendicular velocity components. 2.9.3. Power flux on the top surface. Since the total flux arriving at the surface must be conserved, higher deposition than expected to an infinite planar surface (which is denoted as q i⊥ ) is found on the top surface immediately behind the LE, compensating the smaller integral flux to the LE predicted by both the ion orbit and the PIC models. The LE ion flux pro files in figure 13 show that the power flux profile immediately downstream of the LE increases and broadens significantly in the presence of the electric field. In some cases, strong focus ing of the flux is found in the vicinity of the edge (due to strong electric fields). However, in the case of steadystate plasma conditions, the deposition profile can be well approxi mated by a combination of constant flux due to parallel flow and a half Maxwellian profile as shown in figure 14 for case PG.1 (full PIC and ballistic). The sharp peak near s tor = 0 is caused by strong electric fields in the vicinity of the LE corner and cannot be captured by this simple fitting function q i (s) = q i,⊥ + Ae For the ballistic case v || = c s , however, additional acceleration in the sheath must be considered in the presence of sheath electric fields. As an upper limit, it can be assumed that all the energy gained from the potential drop across the sheath (V s = 3 kT e in all the simulations) is converted into parallel accelera tion: v || = c 2 s + 2eVs mi . This yields λ top,model = 0.6 mm (bal listic) and 1.2 mm (with sheath electric field) for the plasma conditions of case PG.4, close to λ top,fit obtained by fitting the deposition profiles (0.7 and 1.4 mm respectively), as shown in figure 14. It is also possible to determine the sensitivity of λ top,model to τ by introducing the ratio of decay lengths for bal listic and PIC cases: which varies between 2.5 (PG.1) to 2.3 (PG.4). This theor etical ratio is higher than that observed when comparing the fits to simulation results (~2 for PG.1) because the energy gained by ions in the sheath is shared between the parallel and perpend icular velocity components.

Thermal modelling using PIC and ion orbit heat flux
profiles. To assess the impact of the difference in PG LE loading found by the ion orbit and PIC simulations, thermal calculations identical to those used in [1] were made using as input the heat flux profiles obtained from the two codes and running the calculation for a fixed time over which the heat flux is applied. A selected scenario corresponding to the steadystate conditions was performed, with PIC heat flux profile taken from case PG.3S (black curve in figure 8) and q ⊥ ,pk = 9.3 MW m −2 (note that this particular loading case differs from that considered in [2], where the MB is irradi ated by thermal plasma flux following magnetic field lines and perpendicular fluxes due to radiation), which was compared   to the corresponding profile obtained using the ion orbit code. The time evolutions of the hottest point at the MB are shown in figure 15. The difference between ion orbit and PIC cases, which are differentiated by the extent of spreading of the heat flux over the LE, is about 300 °C, representing about 10% of the peak temperature increase.

Summary
Predictions for the loading of misaligned monoblock edges in the ITER tungsten divertor based on the optical (or guiding center) approach or ion orbit modeling (including Larmor gyration but not electric forces) [1] show that surface temper atures during steady state, burning plasma operation would be sufficient to produce recrystallization of the tungsten mat erial in the vicinity of the edge, and melting during slow tran sient increases in the power loading. These calculations lead to the conclusion that MB shaping, namely a 0.5 mm toroidal bevel, is required to protect leading edges across PGs [2]. The increased manufacturing complexity of this solution moti vated the present study. The objective of the work presented in this paper was to perform predictive simulations of the heat flux distribution in the vicinity of PGs between unshaped ITER MBs with exposed leading edges, while taking into account the presence of sheath electric fields, to evaluate the significance of such fields on edge loading and to clarify the underlying physical processes.
Results of PIC simulations for steadystate and slow tran sient plasma conditions show (in agreement with study per formed in [1]) that due to finite Larmor effects, the leading edge created by the radially misaligned tile receives lower heat flux than predicted by the optical model. In addition, the sheathelectric fields drive polarization drift in the magn etic presheath, which accelerates ions and drives the flow supersonic before the Debye sheath is reached. As a result, the particle flux deposition profile on the leading edge can be well approximated by the parallel flux profile above a simple planar surface.
The heat flux that is missing on the leading edge (when compared to the optical model) is redistributed on the top sur face of the MB within 1-2 mm from the edge. The character istic deposition length on top surface is about a factor 2 longer in PIC simulations than predicted by the ion orbit code. This is due to the additional acceleration of particles in the sheath and magnetic presheath. However, this difference is not sig nificant due to the small spatial scale of the profiles. Thermal calculations which make use of ion orbit heat flux profiles instead of PIC profiles result in a 10% difference in the peak temperature increase of the LE. The influence of electric fields in the sheath on the heat flux distributions is thus insignificant for all cases studied and the ion orbit code can therefore be used with confidence to study this problem.