Beryllium erosion and redeposition in ITER H, He and D–T discharges

The Monte-Carlo code ERO2.0 was used to simulate steady-state erosion and transport of beryllium (Be) in the ITER main chamber. Various plasma scenarios were tested, including a variation of the main species (hydrogen, deuterium, helium), plasma conditions (density, temperature, flow velocity) and magnetic configurations. The study provides valuable predictions for the Be transport to the divertor, where it is expected to be an important contributor to dust formation and fuel retention due to build-up of co-deposited layers. The Be gross and net erosion rates provided by this study can help identifying first wall regions with potentially critical armour lifetime.


Introduction
ITER will be constructed with a tungsten (W) divertor and a main chamber first wall (FW) armour made of beryllium (Be) [1]. It is expected that an important contributor to the FW lifetime will be the sputter erosion of Be due to steady state plasma particle fluxes under long pulse, burning plasma conditions [2][3][4]. Migration of the eroded Be also contributes to retention due to co-deposition as well as to dust formation [5]. High fidelity modelling of the Be erosion, transport and redeposition (labelled as Be migration) is required to quantify the wall lifetime and to provide refined estimates of the growth of * 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. in-vessel dust and tritium inventories. Such modelling is also useful for generating synthetic signals to test ITER diagnostic systems intended for the monitoring of erosion and migration, such as the main chamber visible spectroscopy system [6] and the in-vessel viewing system [7]. Numerical predictions must account for the 3D structure of the FW panels, which are shaped to protect plasma impact on leading edges arising from radial misalignments between components [1,8].
The present work continues the previously reported modelling of global Be migration in ITER using the threedimensional (3D) plasma-surface interaction (PSI) and impurity transport code ERO2.0 [9]. This code simulates the steady-state Be erosion flux, taking into account the subsequent kinetic transport of eroded Be impurities, which are represented by so-called test particles in the trace impurity approximation [10]. For each Be test particle, a trajectory is calculated until the particle is deposited. The simulation takes into account reflection of and sputtering due to test particles at the surfaces, ionisation and recombination, collisions with the background plasma (described by a Fokker-Planck term) and anomalous cross-field diffusion using an assumed constant anomalous diffusion coefficient (here: D ⊥ = 1 m 2 s −1 ).
The earlier study is significantly extended in scope by taking into account different plasma scenarios, corresponding to various operating conditions foreseen in the ITER D-T Hmode fusion power operation (FPO) and H or He plasma L-mode pre-fusion power operation phases (PFPO) [11].

Overview of the simulation cases
Eight different plasma cases are used in this work, as described in table 1. Importantly, all cases are assuming steady-state plasma conditions-transient events like ELMs or VDEs are not considered.
For each case, the Be erosion is calculated and the subsequent Be transport is simulated using an ensemble of 10 6 Be test particles, which is a number found to produce satisfactorily low levels of Monte-Carlo noise (e.g. below 1% of the total gross erosion). The necessary computing resources are provided by the JURECA supercomputer [12].
The main difference between the ERO2.0 simulations for cases #1-8 are the input plasma and neutral backgrounds, as discussed in more detail in section 2.4. Other simulation parameters are kept constant unless explicitly stated otherwise below. The first five cases #1-5 correspond to FPO D-T H-mode discharges (treated here as pure D plasmas, see section 2.4), 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 [13]. Note that cases #1 and #2 were presented earlier in references [9,14], respectively, and are now put into the context of a wider parameter study. The cases #6 and 7 correspond to PFPO L-mode discharges in H plasma with P SOL = 20 MW, I = 7.5 MA, B = 2.65 T. Case #8 is a He plasma case for which all plasma parameters were taken from case #7, thus H was replaced by He keeping everything else unchanged. Figure 1(a) shows the Be sputtering yields used in the simulations. For simplicity, only physical sputtering of Be is considered in the present study. It should be mentioned that the so-called chemically assisted physical sputtering, which leads to release of molecules such as BeH or BeH 2 (or featuring the respective isotopes D and T), can be an important contributor to the erosion [15,16] and affect the Be migration characteristics [17].

Sputtering and reflection coefficients
The physical sputtering yields are calculated with the aid of the Eckstein fit formula Y(E, θ) dependent on the particle impact energy E, impact angle θ and a set of fit parameters [18]. The derivation of E and θ used to calculate the erosion is discussed below in section 2.4. For D on Be impact, the Eckstein fit parameters are taken from reference [19]. The fits are based on simulations described in reference [20], which use the codes SDTrimSP [21] (binary collision approximation) and PARCAS [22] (molecular dynamics). Here, we use a particular fit (the so-called 'ERO-min' fit) which assumes a high fuel content of 50% in the Be surface. It was shown in recent studies that this fit leads to a better agreement with JET spectroscopy data, compared to the other opposite extreme assumption of a clean Be surface which leads to generally higher sputtering yields [17]. Since sputtering yields with a similar assumption were not available for H impact, the same yields for D impact were also used for the H plasma cases to avoid inconsistent surface assumptions. This likely leads to a somewhat overestimated erosion in these cases, since according to data from reference [23, p 42], H should lead to lower sputtering yields than D.
For the single He plasma case considered in this work, it is assumed that there is significantly less implantation in the Be than for H or D plasmas. Fit parameters given in the Eckstein 2007 report [23] are therefore used which are based on SDTrimSP simulations for a clean Be surface.
The Be test particles simulated in ERO2.0 can be reflected from the Be or W surfaces. Figure 1(b) shows the particle reflection coefficients obtained from dedicated SDTrimSP simulations on a clean Be surface.
Finally, it should be noted that SDTrimSP and PARCAS typically use several different assumptions regarding the Be surface, as discussed e.g. in reference [24]. For example, different interatomic potentials are used for the Be-Be and Be-D interaction. Furthermore, SDTrimSP assumes an amorphous target [21], whereas the PARCAS simulations used for the 'ERO-min' fit assume an hcp lattice with an open (0001) surface [20]. In realistic reactor conditions, the Be targets will likely contain defects or impurities which are not accounted for in either simulation methods. In addition, the possibility of surface roughness and morphology is not considered in any of the simulations used here, as will be further discussed below.

Wall geometry and materials
The ITER simulations are performed in a 20 • toroidal sector and assume periodic boundary conditions. The simulation geometry is shown in the 3D view of figure 2(a), which also indicates the wall materials, the FW panel numbers counted clockwise from 1 to 18, and the surface cell resolution (∼30-50 mm).
To not over-complicate the already computer memory and CPU-time demanding 3D geometry, fine structures such as castellation of the Be surfaces (introduced to alleviate the effects of thermo-mechanical stress [25]) are neglected. On a local scale, the castellations can influence erosion in a non-trivial way via the particle fluxes and angular distributions [26]. In a similar manner, the surfaces are assumed to be perfectly smooth, ruling out the possibility of a micro-scale roughness which could as well influence erosion [27][28][29].
The poloidal view in figure 2(b) indicates the poloidal projection of the PFCs, together with the two different magnetic configurations used for this work (further discussed below in section 3.4). The figure also marks several locations of interest that are referenced in other figures below.
It should be noted that the wall model is not gas-tight: there are poloidal and toroidal gaps between the tiles and an array of circular protrusions in the panel centres. In addition, the large openings at the outboard equatorial and upper port locations (used for auxiliary heating, diagnostics and test blanket modules [30]) are left empty in the model. In the ERO2.0 simulations, the openings in the FW have the effect that a fraction of the Be test particles can escape the simulation volume in the outward direction. These particles are considered as deposited 'in gaps' and are removed from the simulation. The fraction of particles removed in such way is typically only a few percent for the plasma conditions considered here (see row 8 in table 3), thus the influence on the situation outcome is small. As a simplification for the present study focussed on the main chamber, the modelling does not consider the build-up and re-erosion of Be deposited layers in the divertor, which consists of pure W for the entire simulation. Also W erosion and redeposition in the divertor is not considered.

Plasma and neutral backgrounds
ERO2.0 requires an input plasma and neutral background (PBG), which contains the spatial distribution of plasma parameters in the boundary region: the electron density n e , electron and ion temperatures T e ,T i , the ion fluid flow velocity v and the magnetic field B. Electric fields are calculated automatically in the sheath using the formulas from reference [31], while electric fields outside the sheath were not considered in this study. The input PBG furthermore contains the impinging particle fluxes (neutral and ionised H, D, He) and their mean energies, which are used to obtain the erosion source. The PBG is also required for the Be impurity transport simulations, where the atomic processes and collisions depend strongly on plasma density and temperature.
The PBGs used for this work, provided by the ITER organisation, were calculated using the 2D code package OEDGE [32]. The latter uses plasma solutions from the 2D plasma boundary code SOLPS-4.3 [33] and extends the grid to the wall, using a combination of onion skin modelling and the kinetic neutral transport code EIRENE [34]. In addition, OEDGE uses the so-called ITER heat and nuclear load specifications (HNLS) [1] as a constraint for the far-SOL plasma.
Different plasma species are considered in the respective simulation scenarios (D-T, H, He). In the case of D-T plasmas, the SOLPS-4.3 PBG simulations approximate these as pure D plasmas and He ash from the fusion reaction is ignored, thus the same approximations are also used in the ERO2.0 modelling. Furthermore, Be is considered in ERO2.0 to be the only impurity, though the PBGs include neon seeding for power dissipation in the divertor, which can be neglected for the current study focussing on the main chamber.
Due to the 3D shaping of the FW panels considered in the ERO2.0 modelling, some modifications of the OEDGE data are required. The shaping leads to a gap of a few cm between the outer OEDGE grid and the FW surface cells. In a simple approach, the profiles of n e ,T e ,T i and v are extended to the 3D surface using constant ('flat') extrapolation, as shown in the outer mid-plane (OMP) profiles in figure 3. This kind of extrapolation is the most conservative, since it leads to the highest possible plasma density and temperature at the surface and therefore the highest erosion rate. A discussion of this assumption and effects of less conservative extrapolation methods is provided in reference [14]. Figure 3 also shows that the temperature profiles in the OEDGE extended grid region (marked by the red area) are flat in order to match the highly conservative HNLS, e.g. for case #1: T e = 10 eV and T i = 20 eV.
In the ERO2.0 simulations, the erosion due to impact of background main species particles (H, D or He depending on the case) is treated fundamentally differently from the erosion due to impurities (Be), which are described here solely by the traced test particles and are thus obtained self-consistently in each ERO2.0 simulation. For the background particles in contrast, input from the PBG and some additional assumptions are required, where the model distinguishes between background plasma ions and energetic charge-exchange neutrals (CXNs). In the case of the background neutrals impinging on the FW, the fluxes are toroidally uniform, as shown in figure 4(a). In the case of the background ions however, the magnetic selfshadowing of certain surface areas has to be accounted for. A detailed description of the shadowing model, which is based on the connection length patterns obtained by numerical tracing of the magnetic field lines in 3D, is given in reference [9]. The shadowing model leads to a plasma-wetted area that is only about 1/10th of the total FW area. A comparison of the fluxes of background neutrals and ions is shown in figure 4 on the example of simulation case #1. One can see that the neutral fluxes are generally lower than the ion fluxes by about one order of magnitude. However, neutral fluxes cover a much larger surface area which compensates the lower flux, thus for the total (surface integrated) erosion they are not negligible. This is also shown by the toroidally averaged flux profiles  in figures 5(a) and (b) and the total impact rates provided in table 2. Note that for the He case #8, the ion fluxes are split between He + and He 2+ as discussed below in section 3.5. The erosion by background particles also depends on their energies. The average energies for ions and neutrals are shown as poloidal profiles in figures 5(c) and (d). To calculate the sputtering, Maxwellian distributions around T i are used. Due to the acceleration by the sheath, the ion energies are additionally increased by the sheath potential 3k b T e , thus the average energies shown in figure 5(c) are given by The average neutral energies are provided as a direct output by the EIRENE part of OEDGE. Since their energy spectra were not available, their energy is for simplicity taken here as equal to the average energy at each poloidal location. The average energies of both ion and neutrals are, for most of the poloidal profiles, well above the sputtering threshold o f H, D, He on Be [35].
Finally, some assumptions are made for the impact angular distributions of the background particles, which can have a large effect on the simulation outcome (due to the strong dependence of sputtering yields on the impact angle) as shown by dedicated sensitivity studies [14]. In the current work, the following assumptions are used for calculating the erosion: for ions, the impact angular distributions are obtained for each cell using the so-called sheath tracing module built into ERO2.0. In this module, which is based on the approach described in [31], a number of ions is initialised with a Maxwellian energy distribution at the entrance of the sheath for certain plasma conditions (n e , T e , T i ), magnetic field strength B and inclination angle θ B . The impact angular distributions are then obtained by following the gyro-orbits of the ions between the sheath and the surface taking into account the electric field of the sheath. For example, for typical conditions of n e = 10 18 m −3 , T e = 10 eV, T i = 20 eV, B = 6 eV and θ B = 85 • , the D + ion impact angles are distributed in the range θ = 50 • -90 • with a maximum at θ 65 • -70 • [14].
For neutrals, the angle is set to the magnetic inclination angle. The latter is a strong simplification, which could be significantly improved in the future when numerical impact energy and angular distributions for neutrals become available from EIRENE. Furthermore, for both ions and neutrals the surface was assumed to be perfectly smooth (no roughness considered). A rough surface would, among other effects, influence the angular distributions and thus the erosion. One should note that the present set of assumptions, together with the mostly grazing magnetic field line incidence on the FW (θ B ∼ 80 • -90 • [9]), leads to a high fraction of background particles with grazing impact angles in the range 80 • -90 • where the maximum is found for the sputtering yield curve of D → Be as shown in figure 1(a).

Simulation results
In the following subsections, the simulation results for each case are discussed using different representations of the data. Figure 6 shows the Be net erosion/deposition flux in a 3D view. This detailed representation is useful to see the locations of maximum net erosion or deposition, which are sometimes found at the same poloidal location. Figure 7 shows the same Be net erosion/deposition fluxes as a poloidal profile. This representation is more compact and allows easier comparison of different cases, note however that due to the toroidal averaging some information on the peak fluxes is lost, and in general the averaging shifts the net fluxes closer to zero. Finally , table 3 shows gross erosion and deposition rates obtained after integrating over the entire surface, and also indicates the maximum net erosion fluxes and locations. Note that the net erosion is defined per surface cell, thus the maximum net erosion fluxes and locations are defined for the single surface cell with the highest net erosion (see section 3.1). In general, the maximum net erosion was found to have a statistical uncertainty of up to ∼5% due to Monte-Carlo noise-other quantities given in table 3 have a lower uncertainty of below 1% due to being surface-integrated.

Reference scenario (case #1)
Case #1 is a baseline D-T high-density discharge scenario. However, as mentioned above, the D-T mixture is replaced by pure D in the PBG (as well as for the other D-T cases #2-5). According to the data from reference [23], the sputtering yields for T impact are about factor 1.5 larger than those for D in the relevant impact energy and angle range. Assuming that the same ratio applies to the 'ERO-min' yields, and a homogeneous 50:50 D-T mixture, the total erosion by background species for cases #1-5 is underestimated by factor 1.25.
Significant erosion is found in multiple regions at the inner wall (panels 4-6), outer wall (panels 11, 18) and the apex (panels 8,9). These are also the regions intersected by the secondary separatrix where the main plasma-wall interaction occurs (dashed red line in figure 2(b)). The total Be gross erosion obtained after integrating over the entire 760 m 2 Be FW area amounts to 1.5 × 10 23 Be s −1 . The dominant erosion mechanism is Be self-sputtering, causing 56% of the gross erosion, followed by sputtering via background D CXN (23%) and background D ions (21%).
The  in figure 8. Note that in general, net erosion fluxes are significantly lower than the gross erosion fluxes due to the Be redeposition. In case #1, 90% of the Be eroded from the FW is again redeposited on the FW. The remaining 10% of the eroded Be is deposited in the divertor, of which the major part of 84% is deposited on the outer target and 16% on the inner target. As mentioned above, the build-up and re-erosion of Be layers in the divertor is not considered (making the divertor a sink for the Be test particles). Considering these processes can somewhat change the Be distribution in the divertor.

Influence of the SOL flow (case #2)
Be transport in the ERO2.0 simulations is strongly affected by the presence of background parallel plasma flows via friction and entrainment, described by a Fokker-Planck term. The background flow field v is imported as part of the PBG. Most OEDGE cases in the present study use the 'sourcesink' driven background flows extracted from the SOLPS simulations, which do not include plasma fluid drifts (although drifts are always included for the Be transport in ERO2.0) and may therefore underestimate the background flow velocity. To study the possible effect of a high background flow, the PBG in case #2 uses a high imposed SOL flow v (Mach number M = 0.5 at the inner midplane), with otherwise identical parameters as in case #1 (see references [30,32]). A comparison of the background flows for cases #1 and #2 is shown in figure 9. As one can see, the flow is not zero in case #1, but significantly lower than in case #2 in the SOL region between the inner target and the stagnation point near the OMP.
The most striking difference compared to case #1 is the increase of the fraction of Be deposited in the divertor by about factor 2, now amounting to 21% of the eroded Be. The presence of higher SOL background flows leads to more long-range Be transport due to entrainment. The deposition increases in both inner and outer divertor, but the increase is larger in the inner divertor, since the SOL flow in case #2 is directed predominantly anti-clockwise (see figure 9). Therefore, the distribution between inner (51%) and outer (49%) targets becomes almost equal (thus significantly more deposition occurs at the inner target, compared to case #1). In addition, the lower flux of Be to the FW in case #2 means that the Be self-sputtering of the FW is reduced. Therefore, both gross and net erosion are reduced compared to case #1 (see table 3).

Influence of far-SOL density (case #3)
Cases #1 and #3 can be seen as a SOL high-density case (HDC) and low-density case (LDC), respectively, of the same baseline D-T discharge. With otherwise identical parameters, in case #1 n e and T e /T i at the OMP secondary separatrix are set to 1.5 × 10 19 m −3 and 10 eV/20 eV, respectively, while for case #3 they are set to 0.1 × 10 19 m −3 and 20 eV/40 eV, as seen in table 1 and figure 3. As mentioned in [4], the HDC and LDC represent physics-based limits defined in the ITER HNLS [1], with the LDC representing the case of low convective transport and the HDC representing highly turbulent, filamentary transport resulting in broad 'shoulders' in the SOL density profiles. Note that in a similar way, cases (#4 and 5) and (#6 and 7) are paired as high and low density versions, see table 1. While the higher SOL temperature (and thus particle impact energy) in case #3 increases sputtering yields, this effect is outweighed by the lower fluxes, so that the total gross erosion is reduced by factor 3 to 4.8 × 10 22 Be s −1 compared to case #1.
A significant effect of the reduced density is the stronger long-range transport of Be and thus larger fraction of Be deposited in the divertor, which is now 25% (even higher than in case #2). This can be primarily attributed to the larger penetration depths of the sputtered Be neutrals, before they are ionised and carried towards the divertor. For example, at the OMP location R − R sep = 10 cm in figure 3, the plasma parameters are n e 10 19 m −3 and T e 10 eV for case #1, and n e 10 18 m −3 and T e 20 eV for case #3. Using the ADAS [36] ionisation coefficients for neutral Be and assuming an initial energy of 10 eV for the sputtered Be atoms, we get ionisation lengths of 2.5 cm for case #1 and 21 cm for case #3 (the difference gets larger when closer to the FW, and lower closer to the separatrix where the profiles for cases #1 and #3 converge, as shown in figure 3).

Influence of the magnetic configuration (cases #4 and 5)
Two different magnetic configurations were considered in this study, characterized by the parameter Δr sep , which is the OMP distance between the primary and secondary separatrices ('wall clearance'). Cases #4 (HDC) and #5 (LDC) use Δr sep = 4 cm, which is the minimum allowed value at ITER [1]. The other cases use a magnetic configuration with Δr sep 11 cm, which represents an approximate upper limit designed to reduce the peak stationary heat fluxes experienced on the upper FW in the secondary divertor region (FW panels 8 and 9). However, as shown here (see e.g. figures 6 and 8), this can lead to rather high erosion on the inboard panels 4 and 5. In reality, a compromise in Δr sep will be sought during ITER operation to provide adequate power flux management while minimizing net erosion. Ultimately, these will both depend sensitively on the PBG, which will not be known until ITER operates. The difference between the two configurations with 'narrow' Δr sep = 4 cm and 'wide' Δr sep = 11 cm is shown in figure 2(b).
In the narrow Δr sep configuration, the ion flux is increased on panels 8 and 9 at the machine apex. Also the electron and ion temperature are increased in that region (40 eV and 80 eV, respectively), leading to increased ion impact energies as shown in figure 5. This leads to a substantial increase in the Be erosion and redeposition in that area, as seen in figures 6 and 7(b). Here, regions of strong net erosion and net deposition are found very closely to each other. Both total gross erosion rates and maximum net erosion fluxes for case #4 and #5 are the highest of all cases as shown in table 3. This shows that the narrow Δr sep configuration is not viable for long-term operation from the PSI point of view, if the assumed PBG is a reasonable approximation to what will be found during ITER operation. Interestingly, case #5 has an even higher erosion than case #4 (despite being the LDC), which is attributed to the special role of Be self-sputtering in that case. The higher electron temperature in case #5 increases the sheath potential drop and thus the impact energies of highlycharged Be at the machine apex. Due to the very focussed PSI region, this makes for a strong mechanism of Be selfsputtering (responsible for 96% of the entire erosion in that case) which eventually outweighs the lower D flux (compared to case #4).

Influence of the heating power and plasma species (cases #6-8)
Cases #6 (HDC) and #7 (LDC) correspond to L-mode H plasmas. The low power cases #6 and #7 consequently have lower fluxes and plasma temperatures at the FW compared to cases #1-5. The surface-integrated Be gross erosion is 1.9 × 10 22 Be s −1 and 1.3 × 10 21 Be s −1 for case #6 and #7, respectively, thus about an order of magnitude lower than in the corresponding D plasma HDC and LDC #1 and #3. Note that, as mentioned above in section 2.2, sputtering yields for D impact are used as an approximation of the yields for H impact. According to the data from reference [23], H impact leads to about factor 2 lower sputtering yields than D impact in the relevant impact energy and angle range (the factor is larger for energies below 100 eV, however CXN is the dominant background erosion mechanism for cases #6 and 7, thus only energies above 100 eV according to figure 5(d) are relevant).
The He plasma case #8 uses, as mentioned above, the same PBG as case #7. The approximation is not unreasonable, since dedicated analysis of SOLPS-4.3 solutions for H and He cases shows that the radial n e and T e profiles are in fact very similar for the two cases. Furthermore, the analysis showed that the charge state ratio on the radially outermost SOL rings is approximately constant with He + /He 2+ 0.5. Since in ERO2.0 it is currently not possible to define spatially varying charge state ratios, a constant value He + /He 2+ = 0.5 was used for case #8. The total ion impact rate of 3.4 × 10 21 ions s −1 from table 2 is therefore composed of 1.1 × 10 21 He + s −1 and 2.3 × 10 21 He 2+ s −1 .
As shown in figure 1(a), the sputtering yields for He on Be are significantly higher than for D on Be for most impact angles, except for very shallow impact (θ = 80 • -90 • ) where the D on Be yields are larger than zero due to the 50% fuel concentration assumption. In addition, the impact energies of He 2+ are higher than those of singly charged ions shown in figure 5(c), which further increases the sputtering yields. Nevertheless, the total gross erosion for case #8, being 1.1 × 10 21 Be s −1 , is lower than for case #7. This can be explained by the above-mentioned set of assumptions for the background particle impact angles, which leads to grazing incidence in the central, plasma-wetted parts of the FW panels, resulting in high sputtering by H or D but low sputtering by He. On the other hand, case #8 shows higher erosion at the FW panel edges due to the closer-to-normal magnetic field line incidence in those areas, as seen in figure 6. This does not, however, outweigh the lower erosion in the central parts of the panels. These rather unexpected findings highlight the delicate interplay of the basic sputtering yield data and impact angular distributions, and motivates future, more elaborate studies using angular-resolved CX neutral sources and considering the possible effects of roughness on the particle impact angles.
Finally, one can note that the PFPO cases #6-8 are the ones with the largest fractions of Be deposited in the divertor (the highest being 53% in case #7). As discussed above for case #3, this can be attributed to the longer Be 0 penetration depths due to the lower electron density in the SOL.

Use of the results in synthetic spectroscopy
In order to provide synthetic signals for ITER diagnostics (one of the main goals of the ERO2.0 study described here), the Be 0 and Be + densities in the plasma edge were calculated in eight cases on a regular cylindrical grid with steps ΔR = 1.1 cm, ΔZ = 1.0 cm and Δϕ = 0.1 • . Then, the emissivities of different spectral lines were calculated using the latest available metastable unresolved photon emission coefficients from ADAS [36], which imply a Boltzmann equilibrium. Note, however, that in [37] the initial population of the 3P-metastable level of Be atoms after sputtering turned out to be equal to 1/3rd of the population in the ground state, and not zero, as is assumed in the photon emission coefficients used here. The results obtained for the 20 • toroidal sector are extrapolated to the entire 360 • torus. These emissivities can be used to test the performance of various optical diagnostics aimed at measuring Be fluxes which will form part of the ITER diagnostic set, such as the Vis/IR wide-angle viewing system (WAVS) [38], the main chamber visible spectroscopy [6] or the divertor impurity-flux monitor [39]. Here, the application of ERO2.0 results to synthetic diagnostics is briefly demonstrated through the example case of simulations of the Be I 457.39 nm intensity (transition 2s 1 3d 1 1 D 2.0 → 2s 1 2p 1 1P 1.0 ) for the Vis/IR WAVS filtered camera.
An important characteristic of all optical diagnostics in ITER is their ability to cope with significant levels of stray light due to reflections from metallic surfaces of the FW [40]. Figure 10 shows the ray-traced intensity of the Be I 457 nm emission in the field of view of the ITER Vis/IR WAVS leftview filtered camera in equatorial port 12 in the eight ERO2.0 simulation cases. These synthetic images were obtained using the Raysect ray-tracing framework [41]. The FW light reflection properties are taken equal to those used in [42] for the simulation of D α signals.
In the case of Be emission, the reflected light does not exceed the direct signal anywhere in the areas of high Be erosion in the main chamber. For example, in case #2 on the sight lines that orthogonally observe the plasma-wetted area of FW panel 5, the stray-light-to-signal ratio is about 0.75. This contradicts the previous results in reference [40], where the toroidally symmetric Be density was simulated with the DIVIMP code included in the OEDGE suite of codes [32] applied to simulations for the 2D wall geometry. In case 'd' of reference [40] (see table 1 therein), characterized by low density in the far SOL, the high emissivity of Be in the divertor compared to that in the main chamber resulted in high stray-light-to-signal ratios in the range 3-6 on the same sight lines. In case 'o', characterized by high density in the far SOL, the stray-light-to-signal ratio was about 1.6-2. It should be stressed that in these DIVIMP simulations the divertor plates were assumed to be fully coated with Be. To assess how the Be-coated divertor increases the stray light in the present simulations, an ERO2.0 run for a subcase of case #2 with a Be divertor was carried out. The Be coating of the divertor leads to a factor 6 increase of the power radiated in the divertor in the Be I 457 nm line, but the stray-light-to-signal ratio on the above-mentioned sight lines increased by only 15% from 0.75 to 0.86. The remaining two-fold difference with case 'o' can be attributed to the effects caused by the 3D FW shaping. In the ERO2.0 simulations, the emission is concentrated only near the plasma-wetted side of the protruding part of FW panel 5, thus increasing the useful signal there. Despite the fact that ERO2.0 simulations demonstrate a smaller impact of stray light on the measurements compared to DIVIMP, it is nevertheless still necessary to take stray light into account when interpreting the measured signals to avoid an overestimation of Be erosion. An accuracy analysis for the assessment of global Be erosion via tomographic reconstruction of 3D Be emission profiles in ITER based on the results of these ERO2.0 simulations will be published soon as a separate paper.

Comparison to JET-ILW
Since 2011, JET has been equipped with an ITER-like wall (ILW) with a Be main chamber and W divertor, making it an ideal test bed for Be erosion and migration studies in an ITER-relevant tokamak environment [43]. Net erosion measurements from the ILW-1 campaign are available from long-term samples [44]. As shown in reference [45], by extrapolating the local measurements to the entire 27.5 m 2 JET-ILW Be FW surface area, one can obtain a lower estimate for the total net Be main chamber source in the diverted configuration: 21.2 g, or 1.40 × 10 24 Be atoms. After dividing by the total ILW-1 divertor phase discharge time of 4.51 × 10 4 s, an average Be net erosion rate of Γ net (JET-exp) = 3.1 × 10 19 Be s −1 is obtained. It should be stressed that the long-term samples contain information integrated over an entire experimental campaign with very different discharges, and thus can only serve for a rough comparison with steady-state ITER erosion rates presented here.
In addition, the results from previous ERO2.0 simulations for JET-ILW can be considered [46]. For the H-mode, the predicted net erosion rate is Γ net (JET-ERO2) = 4.9 × 10 19 Be s −1 , which is somewhat higher than the experimental estimate but still in good agreement given the uncertainties.
These numbers can be compared to the total net erosion rates in ITER, which are obtained from table 3 by taking the FW gross erosion rate (row 2) for each case and subtracting the respective fraction of Be redeposited on the FW (row 6). For case #1, the net erosion rate is Γ net (ITER) = 1.5 × 10 22 Be s −1 , thus a factor 17 or 11 higher than the respective experimental and simulated rates for JET-ILW. The LDC #3 is somewhat closer to JET-ILW, with Γ net (ITER) = 1.2 × 10 22 Be s −1 . A closer match is obtained for the low-power PFPO cases #6 and #7, with Γ net (ITER) = 6.4 × 10 21 Be s −1 and 7.3 × 10 20 Be s −1 , respectively. This is because in these cases, the SOL conditions are closer to the JET ILW-1 campaign conditions where relatively low heating powers were used (e.g. P NBI = 6 MW for the JET 'Be monitoring pulse' simulated with ERO2.0) with correspondingly lower plasma temperatures. For instance, the H-mode OMP profile used for the modelling in reference [46] had T e 5 eV at the wall, resulting in lower sputtering yields.
Finally, one can note that both JET-ILW experiments [45,47] and modelling by WallDYN [45,48] and ERO2.0 [46] indicate stronger Be deposition in the inner divertor (predominantly on top of the 'apron') than in the outer divertor. This is contrasted by the present ITER simulations, where the Be deposition rate is higher in the outer divertor for most cases. An exception is case #2, where a high SOL background flow is used and the deposition becomes approximately equal for inner and outer divertor. Therefore, detailed consideration of plasma flows including drift effects could prove important for the Be transport predictions. One can further note that no equivalent mechanism to the strong Be deposition on the inner divertor 'apron' of JET-ILW is found for ITER (including case #2). Instead, the deposition in the inner divertor is focussed in a several cm wide region above the strike line, as seen in figure 6-a finding that seems to be related to the specifics of the ITER geometry and the magnetic configurations considered here.

Comparison to WallDYN predictions for ITER
Comprehensive studies have been also performed in the past with WallDYN to predict the global erosion and redeposition in ITER using the DIVIMP code for trace impurity transport modelling [2,[49][50][51][52][53]. Since a simulation for the same conditions as case #1 was available, a code-code-benchmark between ERO2.0 and WallDYN was performed for that case. To make the two models comparable, dedicated simulations with some adjustments were performed: on the ERO2.0 side, the 3D wall geometry was exchanged for a 2D (toroidally symmetric) wall geometry with the same poloidal bins as used in WallDYN. Furthermore, since the WallDYN model assumed a constant impact angle of 40 • for all particles, including the Be impurities, to account for roughness (based on experience from dedicated roughness studies, see reference [27]), the same was assumed in ERO2.0 (in all other simulations, a perfectly smooth surface was assumed). Finally, the 'ERO-min' sputtering yields were substituted for the ones published in reference [23]. On the WallDYN side, the Be sublimation was switched off, and T i = 2T e was used instead of T i = T e .
After these adjustments, ERO2.0 gives a total Be gross erosion of 2.3 × 10 23 Be s −1 (factor 1.5 higher than in table 3) and WallDYN gives 2.7 × 10 23 Be s −1 . The difference of 15% in the gross erosion can be attributed to differences in the Be transport. It is instructive to compare the fraction of Be eroded in the main chamber that goes into the divertor: 7% for ERO2.0 (three percentage points lower than in table 3) and 3% for WallDYN. Due to the lower Be fluxes on the FW in ERO2.0, the Be self-sputtering is reduced compared to WallDYN and thus the total erosion. Understanding the origin of the differences in Be transport in the two codes requires much more detailed benchmarking, which goes beyond the scope of this work and is suggested for future studies. For now, one can state that the WallDYN and ERO2.0 predictions are in a satisfactory agreement, with the remaining difference in total Be gross erosion being in the frame of the typical modelling uncertainties e.g. due to input sputtering yields or plasma parameters.

Comparison to ERO predictions for ITER
In the past, predictions focussing on Be maximum net erosion and armour lifetime were performed for FW panel 11 using ERO [4,54] (the predecessor of ERO2.0), which in turn was successfully benchmarked against similar predictions by the LIM code [3]. Since each of these studies features a simulation case with the same PBG as case #1, a comparison can be attempted.
The most recent ERO work [54], using the same 'ERO-min' sputtering yields assumptions as the present ERO2.0 work, reports a maximum net erosion rate of 0.04 mm h −1 . The corresponding ERO2.0 maximum net erosion rate found in the present work on panel 11 is 0.02 mm h −1 , thus a factor 2 lower. Detailed analysis has shown that the discrepancy is due to a lower D ion flux to the surface in ERO2.0, which in turn is due to lower electron density and more oblique magnetic inclination angles at the surface. These differences originate from the assumed FW panel 11 surface geometry and its positioning in the machine. In the previous studies this was approximated by analytic formulas [55] whereas an updated geometry from technical drawings was used in the present study. Considering this difference in the input assumptions, the present ERO2.0 net erosion predictions can be seen as consistent with the earlier results from LIM and ERO. However one should stress that the previous studies did not take into account the full FW and were instead local simulations at the level of a single FW panel. According to the present study, the region of overall highest net erosion is found on panels 4 and 5 with a net erosion flux of 2 × 10 21 Be m −2 s −1 corresponding to 0.057 mm h −1 , thus about three times higher than on panel 11. This, however, is due to the use of an equilibrium with an extremely large value of Δr sep , which is a pessimistic use case and is unlikely to be used in real ITER operation.

Summary and conclusions
A comprehensive numerical study using ERO2.0 has been performed to assess the Be FW armour steady-state erosion in ITER for a variety of plasma scenarios. The considered variations of SOL conditions (plasma density, temperature), magnetic configuration, heating power and fuel species lead to large variations over three orders of magnitude in the total gross erosion (1.1 × 10 21 -1.6 × 10 24 Be s −1 ). The variation in the maximum net flux goes even over four orders of magnitude (1.5 × 10 19 -1.7 × 10 23 Be m −2 s −1 ). The lowest erosion is found for the low-power PFPO cases (#6-8), the highest erosion is found for the FPO cases with the magnetic configuration with minimum allowed Δr sep (cases #4 and 5).
For the reference D-T scenario with wide Δr sep (case #1), the predictions are consistent with earlier global modelling by WallDYN, as well as local modelling by LIM and ERO. In the latter case, higher net erosion rates (up to 0.057 mm h −1 ) are found on the inner wall (panels 4 and 5) which was not considered by those studies, which were local simulations concentrating on a single panel in the secondary divertor region and unable to account for the full 3D wall geometry. As discussed in reference [3], such net erosion rates can already be problematic for the ITER FW lifetime. The narrow Δr sep cases #4 and 5, with up to a factor 100 higher maximum net erosion rates (located on the upper FW, panels 8 and 9), would not feasible from the PSI point of view. One should note that in addition to the armour lifetime, the Be erosion will also contribute to T retention due to co-deposition with Be [56] in the main chamber and divertor. However, since this process is highly dependent (among others) on the surface temperature T surf , numerical estimates would require a spatial T surf distribution on the surface which was not considered in the present work.
It is important, however, to note that pessimistic modelling assumptions were used for all cases to generate the background plasmas. As noted above, the two magnetic configurations used here with 'wide' and 'narrow' Δr sep are limiting cases that lead to extreme erosion on the inner and upper FW, respectively. An intermediate value of Δr sep will likely reduce the erosion. In addition, flat temperature profiles in the SOL were assumed (e.g. with T e = 10 eV, T i = 20 eV at the FW surface in the reference D-T scenario) corresponding to highly turbulent, filamentary transport. These assumptions provide useful conservative estimates but are not necessarily realistic-for instance, simulations and experiments from the JET-ILW suggest significantly lower erosion rates.
Significant uncertainties also persist regarding the basic Be erosion yields. These are highly dependent on the surface state, including saturation with fuel, oxidation state and rough-ness. For example, using sputtering yields with Eckstein fit parameters obtained for a clean Be surface (instead of 50% D content as used here) would increase the erosion at normal incidence by about factor 4 [19], but at the same time decrease the erosion at grazing incidence, thus making it difficult to predict whether D surface content increases or decreases the total erosion without further dedicated studies. Molecular dynamics simulations of beryllium oxide sputtering [57] show the occurrence of additional sputtering channels compared to a non-oxidized Be surface [58]-however, since these processes are highly dependent on surface temperature and impact velocity, without additional ERO2.0 simulations it is again difficult to estimate whether the overall erosion would be increased or decreased by the presence of oxygen. In contrast, Be surface roughness can be expected to decrease the erosion by up to a factor of 2, as shown by ion beam experiments and simulations [59] (a factor 2 decrease in erosion was also found by ERO2.0 roughness studies for molybdenum and tungsten [60,61]).
Despite the uncertainties and caveats associated with the modelling reported here, the results can serve as a first guideline to global Be erosion, and also highlight the operational difficulties of working with Be (from a PSI point of view) as well as the need for a better understanding of the SOL plasma. The plasma, currently originating purely from scenario design, can be optimized for PSI e.g. by reducing the triangularity or the value of Δr sep . Additionally, replacement of certain Be FW panels by W at a later stage is under consideration [32].