First Monte‐Carlo modelling of global beryllium migration in ITER using ERO2.0

ERO2.0 is a recently developed Monte‐Carlo code for modelling global erosion and redeposition in fusion devices. We report here on the code's application to ITER for studying the erosion of the beryllium (Be) first wall armour under burning plasma steady state diverted conditions. An important goal of the study is to provide synthetic signals for the design of two key diagnostics: the main chamber visible spectroscopy and the laser in‐vessel viewing systems. The simulations are performed using toroidally symmetric plasma backgrounds obtained by combining SOLPS simulations extended to the wall using the OSM‐EIRENE‐DIVIMP edge code package. These are then further combined with a shadowing model using magnetic field line tracing to provide a three‐dimensional correction for the flux patterns. The resulting plasma wetted area, which amounts to ∼10% of the total first wall area, is in excellent agreement with shadowing calculations obtained with the SMITER field line tracing code. The simulations reveal that the main Be erosion zones are located in regions intersected by the secondary separatrix, in particular the upper Be panels, which are close to the secondary X‐point. For the particular high‐density Q = 10 background plasma case studied here, ∼80% of the eroded Be is found to re‐deposit on main chamber surfaces. The rest migrates in almost equal parts to the inner and outer divertor and is deposited close to the strike lines.

• sputtering in the tungsten (W) divertor, [3] • fuel retention by co-deposition, [4,5] • material mixing [6][7][8] and • dust formation. [9,10] Predictive modelling of the Be source strength and distribution is useful to estimate these effects, and thus provide an input to plasma-facing component and diagnostics design. For instance, for in situ monitoring of the Be source and erosion rates, ITER will respectively use the main chamber visible spectroscopy system [11] and the laser in-vessel viewing system. [12] To make final decisions on the diagnostic configuration and assess system performance, predictions for the three-dimensional (3D) Be light emission sources and the erosion rates are required.
A unique tool that can provide such predictions in a realistic 3D wall geometry is the Monte-Carlo code ERO2.0. [13] This code calculates the erosion flux of wall material (in this case Be) by interpolating a database of sputtering yields produced on the basis of binary collision approximation and molecular dynamics (MD) simulations. Subsequently, it simulates the transport, redeposition and radiative emission of eroded material by following test particle trajectories (ions and neutrals) in 3D, including the chain of collisions and atomic processes, and resolving the detailed 3D gyro-orbits of ions taking into account the surface sheath electric field. The trace impurity approximation is used, which means the impurity concentration is assumed to be low enough that the test particles do not influence each other or the plasma background (PBG). The latter is treated as constant and has to be imported from other plasma boundary transport codes as described below.
In the past, the predecessor code ERO was used to predict the Be erosion in ITER and estimate armour lifetime. [1,2] These previous ERO studies focused on only a single first wall panel of interest, but the massively parallel new code ERO2.0 is required to study global Be erosion and transport in the entire main chamber. The new code has been quantitatively validated against experiments at JET with the ITER-like wall. [13,14] In this contribution, we present results of the first application of ERO2.0 to ITER.

ERO2.0 MODELLING
The simulations are performed in a 20 • toroidal sector and assume periodic boundary conditions. Wall models, which correspond to the geometry in the sector with the diagnostic equatorial port, were provided by the ITER organization (IO) in the form of 3D triangular surface meshes. These are based on an analytic approximation for the surface shaping; fine structures such as castellations are neglected. The meshes were provided at different resolutions (triangle side lengths). For the simulations, a wall model with ∼50 mm edge length was chosen, since tests showed that using finer meshes significantly increases the computational time, but does not notably affect the results. As an approximation in the present initial study, we treat Be as the only impurity in the simulations (the burning PBG also contains neon impurity, used for extrinsic seeding, and He from the fusion reactions). Since the focus is on the main chamber, we also neglect the build-up and re-erosion of Be deposited layers in the divertor, which consists of pure W for the entire simulation. However, it should be noted that taking deposited layers into account should have some influence on main chamber emission due to increased stray light from the divertor. For the plasma parameters, we focus on a single reference scenario: the steady-state phase of the baseline D-T H-mode burning plasma, with power gain Q = 10 and scrape-off layer (SOL) input power P SOL = 100 MW, toroidal plasma current I = 15 MA, and central toroidal field B = 5.3 T. [15] As an approximation, we assume deuterium only for the D-T plasma. Furthermore, ELMs are neglected.
ERO2.0 requires an input PBG to calculate the impinging particle fluxes and energies with which to obtain the erosion source. The PBG is also required for the impurity transport, where the atomic processes and collisions depend strongly on plasma density and temperature. The input PBG used here, provided by IO, was calculated using the two-dimensional (2D) OSM-EIRENE-DIVIMP edge (OEDGE) code as described in Lisgo et al. [16] The OEDGE solution is based on SOLPS case No. #1514, [17] and extends the grid up to the ITER first wall contour using an onion-skin model. In the far-SOL, the plasma parameters are prescribed by the ITER heat and nuclear load specifications. [18] For example, at the entire outer wall the plasma temperatures are fixed at T e = 10 eV and T i = 20 eV.
However, even the extended OEDGE grid does not provide the required plasma parameters for ERO2.0 on all surface regions; the Be first wall panels are toroidally shaped and the OEDGE grid boundary extends only to the plasma-facing F I G U R E 1 Outer midplane (Z = 50 cm) profiles for electron density, electron temperature and ion temperature. The vertical lines indicate the separatrix (dashed) and the wall contour (solid), the latter being also the OEDGE radial grid boundary. Outside this boundary ("mini-SOL"), constant extrapolation is used F I G U R E 2 (a) Connection length L calculated with ERO2.0 by numerical field line tracing. The numbers 1-18 of the panels are indicated. (b) Shadowing patterns obtained with ERO2.0 for a threshold L thr = 8 m. Blue regions are shadowed, yellow ones are plasma-wetted. (c) Shadowing patterns obtained with SMITER, using the same L thr , wall geometry model and magnetic equilibrium ridges. For the retracted surface regions ("mini-SOLs"), extrapolation is required. Here, we assume a constant extension of the profile. This is illustrated in Figure 1, which shows the outer midplane profiles for density and temperature.
Additionally, a 3D correction is applied to account for the fact that some surface regions are magnetically shadowed, and thus experience lower parallel ion fluxes. This correction is based on the connection lengths L, which are obtained in ERO2.0 by numerically following a magnetic field line starting from each surface cell, until either another surface cell is hit or a maximum length (here: 60 m) is exceeded. The resulting connection length pattern is shown in Figure 2a. We select a threshold L thr and define all surface cells with L < L thr as magnetically shadowed. For these cells, we set the impacting ion flux (and hence the eroded Be flux) to zero. Figure 1b shows the shadowing patterns obtained for a threshold L thr = 8 m. This particular value was chosen here since it was also used recently for ITER power load estimations with the codes PFCFLUX [19] and SMITER, [20,21] which rely on a similar shadowing model as that used in ERO2.0. As a benchmark, a SMITER calculation for the same conditions (wall geometry model and magnetic equilibrium) is shown in Figure 1c. The very good agreement with Figure 1b confirms the correct implementation of field line tracing, wall geometry and magnetic equilibrium in ERO2.0.
The shaping of the Be panels and resulting shadowing leads to a wetted area that is only about 1/10th of the entire Be surface area. On the other hand, due to the higher magnetic field angles resulting from the toroidal first wall panel shaping, the parallel flux in the wetted areas increases by about a factor 10 compared to a toroidally symmetric wall as used in OEDGE. Integrated over the entire main chamber, OEDGE gives a D + impact rate of R +in D = 5.4 × 10 23 ∕s, while ERO2.0 gives R +in D = 5.0 × 10 23 ∕s. This mismatch is an artefact resulting from the interface between the 2D and 3D codes. To maintain particle balance between outgoing D + in OEDGE and incoming D + in ERO2.0, we multiply the D + flux for all surface cells in ERO2.0 by the constant factor 1.08 so that the total rate of 5.4 × 10 23 /s is obtained.
To calculate Be erosion, we take into account physical sputtering by D + ions, D 0 energetic charge-exchange neutrals (CXN), and self-impact from Be Z+ . Chemically assisted physical sputtering with release of BeD [2,8] is neglected here.
For the physical sputtering yields, we proceed as described in Borodin et al. [2] and use the Eckstein fit formula, [22] where the yield Y (E, ) is a function of impact energy E and angle . The fit parameters [23] are based on MD runs and SDTrimSP [24] runs with 50% D surface content ("ERO-min"). Figure 3a,b shows the electron density and temperature mapped to the 3D surface. Figure 3c shows the magnetic field line angle (measured from the surface normals), which is typically very shallow with angles close to 85 • -90 • . For sputtering by D + , the impact velocity distributions are obtained from preliminary runs of D + sheath trajectories, as described in Borodkina et al. [25] In this way, effective sputtering yields are obtained for each surface cell which depend on the local plasma density, temperature, and magnetic field. Figure 3d shows the resulting gross flux of Be eroded in the main chamber by D + impact. The surface-integrated gross erosion rate is R eroded Be←D + = 3.6 × 10 22 ∕s. The peak erosion flux, with about Γ eroded Be←D + = 1.7 × 10 21 ∕m 2 s, is found on the top of the machine on first wall panel number 8 (counting clockwise from inner to outer wall, as shown in Figure 2a). Other regions with high erosion are panel 5 (inner wall) and panels 11 and 18 (outer wall). These are the regions that are touched or intersected by the secondary separatrix [1] .
The erosion by impact of CXN is treated somewhat differently. Their flux is directly derived from OEDGE, since they are not affected by shadowing. For that reason, the CXN impact affects the entire surface area including the shadowed regions. Therefore, the CXN total impact rate of 3.8 × 10 23 /s is comparable to that of D + , even though the flux is much smaller than that of D + in plasma-wetted regions. We assume mono-energetic impact for each surface cell, with energy equal to the neutral temperature for that surface cell, which is provided by the EIRENE part of OEDGE. Furthermore, we assume that the ions are moving along the magnetic field line in the direction of the surface before the charge-exchange reaction takes place, and therefore the CXN impact angle is equal to the magnetic angle. The Eckstein fit parameters are the same as for D + . Figure 3e shows the resulting gross erosion flux. The most affected panels are numbers 4, 11, and 18. The surface-integrated erosion rate is R eroded Be←D 0 = 2.7 × 10 22 ∕s. Finally, for Be erosion by self-sputtering, the impact fluxes and velocity distributions are directly obtained as a result of Be trajectory tracing in ERO2.0. The corresponding "ERO-min" Eckstein fit parameters are again taken from Borodin et al. [23] . For sufficient Monte-Carlo statistics, an ensemble of 10 6 Be test particles was traced. Figure 3f shows the self-sputtered Be flux. It is strongly concentrated at the top of the machine at panels 8-9. The surface-integrated self-sputtered erosion rate is R eroded Be←Be = 4.8 × 10 22 ∕s. Thus, the total gross erosion rate from all three contributions is R eroded Be,tot. = 1.1 × 10 23 ∕s, to which D + , D 0 and Be Z+ impact contribute 33, 24, and 43%, respectively. Figure 3g shows the total gross erosion flux from all three contributions. The peak gross erosion is on first wall panel 8, where particularly strong self-sputtering occurs, with a flux Γ eroded Be,tot. = 1.1 × 10 22 ∕m 2 s. Figure 3h shows the Be deposition flux obtained by trajectory tracing. To account for reflection of Be on Be or Be on W, SDTrimSP reflection coefficients are used. About ∼20% of the Be eroded in the main chamber is deposited in the inner or outer divertor, the remaining ∼80% are redeposited in the main chamber. The inner and outer divertor each receive about 50% of the deposited Be. In the inner divertor, the deposition occurs mainly within a ∼40 cm broad stripe above the inner strikeline. In the outer divertor it is more homogeneously distributed. In the main chamber, the deposition pattern mostly resembles the erosion pattern. Figure 3i shows the net Be flux obtained by subtracting the gross erosion from the deposition flux. One can see that in the main chamber, most wetted areas are net erosion zones. An exception is panel 8, which is the zone of highest gross erosion, but becomes a net deposition zone due to the very strong Be transport to and redeposition in that region. In the other wetted areas, net erosion is also significantly reduced compared to gross erosion. The strongest net erosion flux is found on panel 5 with Γ eroded Be,tot. = 1.2 × 10 21 ∕m 2 s, corresponding to a layer erosion rate of 0.036 mm/hr. Apart from the erosion and deposition patterns, ERO2.0 also provides the spatially resolved impurity particle and line emission densities. These are critical quantities required for assessment of the ITER main chamber visible spectroscopy diagnostic performance. As an example, Figure 4a shows the Be particle density (all charge states), averaged in the toroidal direction. The dashed lines show the separatrix (green) and the inner simulation boundary (white). The latter is chosen at a normalized plasma radius = 0.9 ( = 1 is the separatrix), and test particles crossing this boundary are reflected. Figure 4b,c shows the Be I 457 nm and Be II 527 nm emission density. Due to the short ionization lengths, the emission occurs close to the regions of strong erosion or deposition. While in the toroidally averaged view the Be I and Be II lines look very similar, the patterns are more complex in 3D and a detailed analysis of synthetic spectroscopy using ray-tracing techniques is required. This analysis is now ongoing at the NRC Kurchatov Institute by the ITER Russian Federation Domestic Agency team responsible for the design and supply of the main chamber visible spectroscopy diagnostic.

CONCLUSIONS
The Monte-Carlo code ERO2.0 was applied for the first time to ITER for predictions of Be erosion, redeposition, and emission in the main chamber. The predicted overall peak net erosion rate is 0.036 mm/hr on panel 5. On panel 11, which was in the focus of earlier ERO studies, [1,2] the erosion is 0.024 mm/hr. Note that in ERO, the concentration of background Be plasma impurities for self-sputtering and deposition was a free parameter that led to a large uncertainty in the predicted Be layer evolution: in Borodin et al. [1] 0% Be leads to net erosion on panel 11, while 3% Be leads to net deposition. In contrast, ERO2.0 calculates the spatially varying Be concentration in a self-consistent way, which removes the uncertainty. Note also that in the ERO2.0 simulation about 80% of the Be is redeposited in the main chamber and self-sputtering is the strongest source of erosion, which explains the importance of taking into account global migration modelling for net erosion/deposition predictions. It is important to note that the results presented here have been obtained from the first ever attempt to apply the new ERO2.0 code to ITER. While the 3D geometry of the first wall is fully accounted for, the simulation is for a single input PBG and a set of fixed modelling parameters. In particular, the very high far SOL ion temperature, T i = 20 eV, is the most conservative value in the ITER Heat Load Specifications. Together with the assumptions that this high temperature propagates out to all wall surfaces at a constant plasma density (fixed by the last density point on the PBG extending to the nominal 2D surface), this leads to very high Be sputtering yields which are unlikely to be found in reality. Such rather pessimistic plasma specifications are fixed by the requirement to provide the worst case thermal plasma heat loads on the most exposed first wall panels for the purposes of thermohydraulic engineering design. They necessarily lead to rather extreme erosion, which would imply correspondingly short rather unrealistic component lifetimes.
For a meaningful discussion, the results must be put into context of a wider parameter study addressing these uncertainties regarding far-SOL conditions. Additionally, different plasma scenarios should be compared, since changing parameters such as the magnetic equilibrium, the heating power, and the plasma species (e.g. hydrogen and helium plasmas, which will be used in the ITER pre-fusion power operation phase [26] ) will lead to very different Be erosion and migration. Finally, one should also address uncertainties related to free parameters in the ERO2.0 modelling itself: • Constant extrapolation of radial profiles in the "mini-SOLs" was used. Choosing a decaying profile would significantly reduce the erosion.
• A constant diffusion coefficient of D ⊥ = 1 m 2 /s was used to describe anomalous cross-field transport of Be. Changing this value should modify the redeposition patterns and the emission.
• Different assumptions could be used for the impact velocity distributions of D + ions and D 0 CXN, for example, using impact angles other than = B for CXN.
Studies in which such parametric variation is performed are underway and will be the subject of future publications.