# SOLEDGE3X full vessel plasma boundary simulations of ITER non-active phase plasmas

N. Rivals, P. Tamain, Y. Marandet, X. Bonnin, J.-S. Park, H. Bufferand, R.A.

Pitts, G. Falchetto, H. Yang, G. Ciraolo

► To cite this version:

N. Rivals, P. Tamain, Y. Marandet, X. Bonnin, J.-S. Park, et al., SOLEDGE3X full vessel plasma boundary simulations of ITER non-active phase plasmas. Nuclear Fusion, 2025, 65 (2), pp.026038. 10.1088/1741-4326/ada2a8. cea-05032697

HAL Id: cea-05032697

https://cea.hal.science/cea-05032697v1

Submitted on 14 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# SOLEDGE3X full vessel plasma boundary simulations of ITER non-active phase plasmas

N. Rivals<fnref n="1" />,<fnref n="2" />,* ID, P. Tamain<fnref n="1" />, Y. Marandet<fnref n="2" /> ID, X. Bonnin<fnref n="3" />, J.-S. Park<fnref n="4" /> ID, H. Bufferand<fnref n="1" /> ID, R.A. Pitts<fnref n="3" /> ID, G. Falchetto<fnref n="1" />, H. Yang<fnref n="1" /> ID and G. Ciraolo<fnref n="1" />

E-mail: nicolas.rivals@cea.fr

Received 12 April 2024, revised 26 November 2024

Accepted for publication 22 December 2024

Published 10 January 2025

<figure><img src="image_2.png" /></figure>

CrossMark

## Abstract

The onset of detachment in the ITER machine is analyzed in this work through the help of 2D-axisymmetric boundary plasma simulations with the SOLEDGE3X-EIRENE code, which features a numerical domain for the plasma solver extending up to the first wall. The plasma boundary is computed in scenarios from the first non-active phase of ITER, in pure H and at 20 MW. This set of simulations is used in two aspects: first, to study the plasma detachment in the divertor, and second, the plasma conditions, fluxes, and beryllium erosion at the first wall. Here, the code results are also compared to those obtained with the well-established SOLPS-ITER code, which includes a plasma numerical domain only covering the main SOL. Results show an increase in the SOL width $\lambda_q$ with increasing density, and a detailed analysis is carried out, for the first time, on each of the different plasma-neutral interactions in the code's physics model in EIRENE. The gross beryllium erosion rates of first wall panels are estimated from 2D simulations, with the aim of assessing their sensitivity to two parameters: the divertor density regime, and the presence of density shoulders in the far-SOL formed by enhanced perpendicular transport at this location. The erosion contributions from neutrals and ions are considered in each case, and the charge-exchange atoms fluxes and energy distributions are provided, highlighting the two atom populations (cold and charge-exchange).

Keywords: edge plasma, simulation, ITER, first-wall, scrape-off layer, plasma-neutral interactions

(Some figures may appear in colour only in the online journal)

<figure><img src="image_3.png" /></figure>

<footnote n="1">¹ CEA, IRFM, F-13108 Saint-Paul-lez-Durance, France</footnote>

<footnote n="2">² Aix-Marseille Univ, CNRS, PIIM, Marseille, France</footnote>

<footnote n="3">³ ITER Organization, CS 90 046 - 13067 St Paul Lez Durance Cedex, France</footnote>

<footnote n="4">⁴ Oak Ridge National Laboratory, Oak Ridge, TN 37831, United States of America</footnote>

<footnote>* Author to whom any correspondence should be addressed.</footnote>
# 1. Introduction

## 1.1. Context & motivation

Understanding the plasma boundary is critical for designing future nuclear fusion power plants, as interactions between the plasma and the machine's wall determine both the plasma performance and the lifetime of in-vessel components. For ITER, many efforts revolve around the numerical simulation of the behavior of its Scrape-off Layer (SOL) in the foreseen operating scenarios. Recently, extensive numerical simulation work has been carried out with the SOLPS4.3 and SOLPS-ITER codes [1, 2]. Notably, The divertor's performance was also assessed with the SOLPS-ITER code for the first, non-active Pre-Fusion Power Operating phase (PFPO-1) of the machine in [3].

In this work, we provide a study based on simulations of the same PFPO-1 operating phase at 20MW of injected power [4] with the SOLEDGE3X-EIRENE [5] code which describes the boundary plasma in the whole vessel, i.e. up to all plasma-facing components. This study builds upon the recently expanded capability of the code to simulate ITER plasmas [6, 7]. It uses a fluid description for the plasma, via equations based on the work of Zhdanov [8, 9], and a kinetic Monte-Carlo approach for plasma-wall and plasma-neutral interactions, through coupling to the EIRENE code [10]. While previous works on simulating plasma conditions at the first wall of ITER [11] relied on extending a fixed main-SOL solution (which includes the divertor targets) obtained with the SOLPS code with a dedicated tool [12] to model plasma in the far-SOL up to the first wall, SOLEDGE3X-EIRENE can describe the plasma self-consistently in both regions of interest: the divertor, and at the first wall in the main chamber<fnref n="5" />.

This study analyzes a set of SOLEDGE3X simulations along three topics:

* a comparison with results obtained with SOLPS-ITER

* a study on the divertor detachment and the plasma-neutral interaction processes leading to it

* a sensitivity analysis of erosion of the first wall to the change in divertor regime and the formation of density shoulders in the far-SOL.

## 1.2. Simulation setup

The simulation set consists of SOLEDGE3X-EIRENE cases at 20 MW of input power in pure H (i.e. without intrinsic or seeded impurities), run in 2D axisymmetric transport mode, without fluid drifts. The main code parameters with which the presented cases were run are summarized in table A1.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** SOLEDGE3X-EIRENE plasma domain mesh for the ITER machine. Left: the plasma solver SOLEDGE3X domain, with a field-aligned mesh, locally aligned with the targets at the active strike points (alignment poloidal distance: 46 cm at the inner target and 35 cm at the outer target), as illustrated in the insert at the outer target. Right: the EIRENE domain. The red segment represents the pumping surface, the green triangles represent gas puffing locations. Cyan dashed lines are semi-transparent surfaces with a reflection probability of neutral particles of 50%, used to represent the non-axisymmetric dome support structures (pillars).</figcaption></figure>

The flux-surface aligned mesh (see figure 1) is based on the magnetic equilibrium of PFPO discharges using the baseline $q_{95} = 3$ ($I_p = 5$ MA, $B_T = 1.76$ T) scenario, in which the geometry of the field lines is identical to the one in the Full Power Operation (FPO) $Q = 10$ plasma scenario. Only the magnitude of the magnetic field is scaled down by a factor of 3 in PFPO-1. This scaling does not affect the considered transport simulations, as fluid drifts are not included. However, it indirectly impacts the choice of transverse transport coefficients which depends on the targeted SOL width, itself expected to scale as $1/I_p$ [14]. In our cases, the values of the prescribed perpendicular diffusion coefficients were selected to match current ITER projections for L-mode SOL width and are also identical to the ones used in the ITER-SOLPS simulation database [3].

The atomic and molecular model used in EIRENE is the reaction set presented by Kotov *et al* in [15] and newly implemented in SOLEDGE3X. It contains additional processes which were not included in the reduced set used until now in SOLEDGE3X: the molecular charge exchange, molecule-ion and neutral-neutral elastic collisions. Those reactions (available in EIRENE since the early/mid-nineties) were selected and included in this newer set specifically to better describe divertor detachment in JET and ITER with the SOLPS code [15]. The following section will highlight that indeed, also in SOLEDGE3X, those processes are found to be important in high-density regimes.

The detachment study presented here (section 3) uses the gas puff rate as scanned parameter, ranging from $1.46 \times 10^{21}$ to $1.76 \times 10^{23}$ e s⁻¹.

<footnote n="5">⁵ The SOLPS-ITER code has been recently improved to include a 'wide-grid' numerical domain now also extending to the first wall [13], so that such numerical simulations up to the first wall should also be achievable in the near future with SOLPS-ITER.</footnote>
Both the first wall material and divertor targets are modeled with Beryllium surfaces, assuming eroded Beryllium particles from the first wall will be transported with the plasma flows and redeposited at the target plates without being significantly re-eroded to expose the tungsten (W) below, which is in line with assumptions in the burning plasma database used in for the divertor physics basis [2, 16]. In fact, the present study was executed before the recent consideration of an W First Wall in ITER, however even in the case of W First Wall, a study carried out with SOLPS-ITER in [3] showed only minor differences when the divertor target surfaces are changed from Be to W, so the results found here with Be targets should also still remain of interest.

### 1.3. Improvements to the code

Even though the SOLEDGE3X code has been widely used to model medium-sized tokamaks [5, 17–20], running ITER simulations proved to be a challenge significantly harder to tackle. Plasma parameters for PFPO-1 cases are not entirely different from what is obtained in high-power experiments for JET, therefore, a probable root cause for the difficulty in simulating the ITER plasmas with the code considered here is the size of the machine. The following sections briefly discuss the rationale for putting forward this assumption. While this should be confirmed by a more in-depth analysis that is outside the scope of this work, the actual observed effect was the considerable increase and sensitivity of collisional processes (both electron-ion exchanges and plasma-neutral interactions), which proved to be challenging to describe. Practically speaking, this translated into very slow and unstable simulations, rendering it impossible to converge cases. Thus, some elements both in terms of physics description and numerical scheme were revised, and are briefly summarized here. The reader may refer to [21] for details on those.

* A ‘multi-grid’-like acceleration scheme: meshes with different resolutions are used, a converged solution is quicker to obtain on the coarse mesh to approximate the solution, which is then interpolated on a finer mesh. The case is then rerun until convergence on this finer mesh starting from this approximate solution as an initial condition. This produced a gain of around a factor of 3 on the computation time.

* An implicit scheme for the computation of collisional terms in the fluid closure: those terms include the collisional energy equipartition (here between the ions and the electrons), and the friction forces between species<fnref n="6" /> which were found to become very large in the lower parts of the divertor, requiring time steps of the order of nanoseconds to describe the plasma parallel fluid velocity adequately. Their expressions in Zhdanov’s work [8] have been linearized and integrated into an implicit scheme in the code.

* An improved gradient diffusion operator with slope limiters following [22] was implemented to prevent anti-diffusive behavior in the case of large gradients diagonal to the mesh. This improves the robustness of the numerical scheme especially near the targets, where the local alignment of the plasma solver mesh to surfaces can generate very slanted cells.

* The number of calls to the EIRENE Monte-Carlo solver was optimized using a scheme that assesses the characteristic time over which the EIRENE sources can be kept constant.

* The transport of the $H_2^+$ molecular ions was included in its description in EIRENE, as their mean free path, in this atomic & molecular model setup, was found to be much larger than the cell size above the reflector plates in the divertor. This became important because of the accumulation of $H_2^+$ ions in these regions, whose density becomes significant compared to the main ion density. This may lead to further investigations in future work regarding the relevance of the inverse molecular charge-transfer process ($H_2^+ + H \rightarrow H_2 + H^+$) that is not included here and which would consume those molecular ions, and secondly, the impact of such high $H_2^+$ density on the electro-neutrality assumption on the plasma solver side, which does not take them into account.

## 2. Comparison with SOLPS-ITER results

To provide some context on the results obtained with SOLEDGE3X, a comparison is carried out with the existing SOLPS-ITER database on similar cases. Six cases with different throughput values (see table 1) are selected for this exercise. The relevant SOLPS-ITER cases were taken from the ITER public IMAS database and include those described in [3] (the last four cases). They are run with the same magnetic equilibrium, which corresponds to the scenario considered prior to the 2021 update, with a distance between the main and secondary separatrices at the outer midplane $dr_{sep} = 11$ cm (i.e. a different equilibrium than the one that is used for the detachment study in the next sections, which uses the $dr_{sep} = 6.5$ cm updated equilibrium). The computing time for those cases was of the order of a few weeks, similar to SOLPS-ITER runs.

### 2.1. Setup matching between both codes

In terms of setup, SOLEDGE3X parameters were chosen to match as closely as possible the SOLPS-ITER runs, both in terms of included physics and parameter values. Most notably, fluid drifts are not activated, and anomalous perpendicular diffusion coefficients for particles, momentum, and heat ($D_\perp, v_\perp, \chi_{\perp e}, \chi_{\perp i}$) are uniform over the computational domain and their values matched between both codes. Sheath heat transfer coefficients and heat flux limiters values for the parallel conductivity are also identical. $H_2$ molecules are injected from the same location (top port), and pumping is modeled with a fixed recycling coefficient (<1) surface below the dome. The list of most important parameters in the runs is summarized in table A1.

<footnote n="6">⁶ In the considered case here there is no friction force accounted for because of the absence of impurities. However, this improvement was necessary for Neon-seeded cases that are a work-in-progress at the time of writing.</footnote>
<table><caption>**Table 1.** Identifiers of the SOLPS-ITER cases considered for the comparison with SOLEDGE3X. Gas puff ($H_2$ molecules) injection rates in electrons per second. For molecules per second, these values need to be divided by two. The corresponding SOLPS case numbers in the IMAS database are provided in the first column. The case at $1.76 \times 10^{23}$ was excluded from the comparison due to the oscillatory nature of the divertor neutral pressure in the SOLPS-ITER run.</caption><thead><tr><td>SOLPS IMAS shot</td><td>Gas puff rate [e s<sup>-1</sup>]</td></tr></thead><tbody><tr><td>103 312</td><td>1.46 &times; 10<sup>21</sup></td></tr><tr><td>103 313</td><td>3.81 &times; 10<sup>21</sup></td></tr><tr><td>103 027</td><td>8.85 &times; 10<sup>21</sup></td></tr><tr><td>103 028</td><td>1.13 &times; 10<sup>22</sup></td></tr><tr><td>103 029</td><td>3.31 &times; 10<sup>22</sup></td></tr><tr><td>103 030</td><td>6.80 &times; 10<sup>22</sup></td></tr></tbody></table>

While the best effort was made to make these runs comparable, it is still worth mentioning that these two codes are widely different in terms of algorithmic implementation, numerical discretization, and treatment of the underlying equations. Consequently, both approaches do not involve the same types of approximations. Therefore, while a perfect match cannot be reasonably expected, the objective is to estimate how much the two codes agree on the selected runs.

The following list identifies the main differences between the two codes that might impact simulation results. Other differences surely exist, but their identification would require meticulous analysis of both codes, an exercise that is outside the scope of this work.

(i) **Full vessel numerical domain:** the considered version of SOLPS-ITER (for the purposes of this study) does not extend radially to all the wall components and therefore applies a radial boundary condition consisting of a set decay length of 3 cm for densities and temperatures into the far-SOL domain. The potential effect of having such boundary conditions is not trivial to assess. It could be estimated by comparing cases with and without the newly developed 'extended grid' feature of SOLPS-ITER [13], but cannot be done with SOLEDGE3X as SOLEDGE3X currently does not include the possibility of restricting the domain before the wall.

(ii) **Parallel currents:** SOLPS-ITER cases solve for the parallel current, while this option is deactivated in the SOLEDGE3X cases, and ambipolarity is assumed. The presence of the parallel current impacts the electron convective heat flux between target plates, and produces limited heat transfer from electrons to ions.

(iii) **Sheath heat transmission:** In SOLEDGE3X the boundary conditions for plasma energy fluxes onto surfaces in the parallel direction use fixed prescribed values for the sheath heat transmission coefficients (4.5 for electrons and 2.5 for ions, see table A1), while the one for electrons in SOLPS-ITER depends on the plasma electrostatic potential, which is not solved for in SOLEDGE3X cases.

(iv) **Core particle outflux:** The SOLEDGE3X runs use a constant particle outflux from the core at $6.13 \times 10^{20}$ ion/s to represent pellet fueling, in addition to the neutral

particle flux reaching the core boundary of the numerical domain which is re-injected as ions. SOLPS-ITER uses lower and different values for the prescribed additional flux from the core, depending on the case. However, this core flux amplitude is only comparable to the puff in the two lowest throughput cases (103 312, 103 313), and therefore those differences should have little impact on the solutions, especially in the SOL. They might, on the other hand, explain the different density levels observed in the closed flux surfaces and at the separatrix.

(v) **Wall temperature:** In SOLEDGE3X, EIRENE is set up so that the wall temperature is uniform at 500 K, while in SOLPS-ITER, the majority of the wall is at 580 K, with a peaked profile up to 1160 K at the strike points. The wall temperature changes the energy of the recycled molecules that thermalize with the wall (but not the velocity of the fast specular reflection of atoms). This temperature difference is very low with respect to the plasma temperatures, therefore the impact is minimal on the plasma response. A few test cases with different wall temperature confirmed that this parameter has very little impact, even on the divertor neutral pressure.

(vi) **Molecular ion transport & $H_2^+$-$H^+$ elastic collisions:** As mentioned in section 1.3, including the transport of $H_2^+$ molecular ions proved to be necessary to stabilize the code, especially during transients. In SOLPS-ITER runs, a static approximation is used. The transport of $H_2^+$ molecular ions also activates $H_2^+$-$H^+$ elastic collisions, but their contribution is found to remain minor (see section 3.4). At the time of writing, it is not possible in EIRENE to activate the tracking of molecular ions without $H_2^+$-$H^+$ elastic collisions, so they are included in SOLEDGE3X-EIRENE simulations even though not present in the SOLPS-ITER runs.

(vii) **Perpendicular diffusion stencil:** The perpendicular diffusion operator used in the SOLEDGE3X code for radial diffusive fluxes is based on a 9-point stencil (which was improved in the context of this work with the addition of slope limiters, as mentioned above) while SOLPS-ITER cases were run using a 5-point stencil scheme (a 9-point stencil option was recently implemented in SOLPS-ITER [23]). The 9-point stencil is expected to reduce numerical errors in the case of very slanted cells, as is the case near the outer target and baffle (in both codes), which could lead to different diffusive behaviors at the outer target.

The next sections compare the results of the two codes across three locations: midplane conditions, conditions and fluxes at the targets, and neutral pressure in the sub-divertor volume. The cases range from attached conditions at low throughput to partially detached conditions at high throughput, where the varying parameter for the scan is the gas puff rate.

## 2.2. Plasma presence in the far-SOL regions

Figure 2 presents 2D plots of the plasma density and the electron temperature for a case after the rollover to provide an
<figure><img src="image_5.png" /><figcaption>**Figure 2.** 2D plots of density, electron temperature and volume radiation for a medium-throughput case (after rollover, $3.31 \times 10^{22} \text{ e}^{-} \cdot \text{s}^{-1}$). Black rectangles indicate the regions of the first wall where plasma-wall interactions are expected to be significant.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Density and temperatures at the outer midplane. Solid lines represent SOLEDGE3X cases, dashed lines represent SOLPS cases. The first row is on a linear scale for the vertical axis, the second row is on a log scale.</figcaption></figure>

illustration of the solutions. This also highlights plasma-dense regions in the far-SOL (black rectangles) which would be outside the plasma solver's numerical domain for non-extended grid SOLPS-ITER simulations. Significant plasma densities are present near the outer bottom first wall plate above the outer baffle and at the machine top near the secondary X-point. Consequently, one can expect increased plasma-wall interaction at these locations, as will be seen in the following sections. In the divertor, the plasma extends radially in the private flux region only up to ~10 cm, but even this small layer made the improvements to the code mentioned in the

previous section necessary. The radial extension of the main chamber plasma towards the first wall is further discussed further on.

## 2.3. Outer midplane conditions

Figure 3 shows that the obtained profiles at the outer midplane are quite similar, and separatrix values are within 15% of difference as shown in figure 4. Separatrix temperatures seem to match, while SOLPS obtains slightly hotter values of $T_i$ in the far-SOL. This may be an effect of the finite radial extension
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Left: Outer midplane separatrix density as a function of throughput (gas injection rate) in SOLEDGE3X and SOLPS simulations. Right: Outer midplane separatrix electron and ion temperatures as a function of separatrix density.</figcaption></figure>

of the SOLPS-ITER grid in the SOL (10 cm at the midplane, vs. 15 cm distance to the FW in SOLEDGE3X), where a prescribed 3 cm decay length of $T_i$ and $T_e$ is applied at this location, and the decay length found in SOLEDGE3X and in both codes at the start of the far-SOL is rather of the order of 1.5 to 2.0 cm.

The left part of figure 4 shows the relation between the throughput and the outer midplane separatrix density. One important feature to note is the strong saturation of the separatrix density towards ~1.9 × 10¹⁹ m⁻³ with increasing throughput obtained in both codes. With a two-point model analysis, one can find that the outer midplane density grows with the square root of the divertor neutral pressure, and experimental data in ASDEX-Upgrade found a relation with an exponent of 0.3 [24], i.e. with an even stronger saturation (vs. 0.5 in the two-point model). Here, since the throughput is roughly linear with the divertor pressure (not plotted here), a similar relation is found, but with an even lower exponent of 0.2. This faster saturation in ITER than in ASDEX-Upgrade may be an effect of the machine size, where the particle source is better confined in the divertor. This observation can also be interpreted as rising SOL opacity, which was observed experimentally in the Alcator C-Mod tokamak and retrieved in SOLPS-ITER simulations in [25].

While the latter graph illustrates that the midplane density response to throughput is not exactly the same in the two codes (which could also be interpreted as a difference in the resulting fueling efficiency), the trend is almost identical.

to SOLPS-ITER may be a product of the absence of parallel current in the SOLEDGE3X cases. However, this case is less relevant to the machine's operational space as target temperatures are too high and would induce too much sputtering. On the inner target, locally at the strike point, SOLEDGE3X shows 30% to 50% lower $T_e$ (e.g. 0.6 eV vs. 1.1 eV for SOLPS-ITER case #103 027) around the rollover, but are more different for the most attached case on the outer target (10 vs 30 eV at the strike point, a factor 2 on the density, with a slightly higher midplane density). While at the strike points, ion temperatures are similar, when moving away along the baffles, SOLEDGE3X produces a factor two or three lower $T_i$ for attached cases. This difference is strongly reduced when increasing the throughput. However, in low-density cases, particle fluxes in these regions are low, leading to a minor impact on the plasma as well as on plasma-wall interaction. Notably, this increased $T_i$ in SOLPS-ITER is in line with the differences found at the mid-plane in the far SOL. Since the values for the sheath heat transfer coefficients have been matched, this difference may be due to the imposed decay length radial boundary condition in the far-SOL in SOLPS-ITER or a difference in the plasma-wall interaction model is unknown at this time. The radial decay length in SOLPS-ITER may lead to lower radial ion and electron heat fluxes. In this case, a stronger discrepancy for $T_i$ can be expected than for $T_e$ (as is the case here), because of ions' much slower parallel thermal conductivity compared to the electrons', leading to lower parallel heat losses at the wall for ions.

## 2.4. Target plasma conditions

At the targets (figure 5), density and electron temperature profiles match very well across the scan except for the lowest density case. For this particular very low-density case (blue), the absence of asymmetry in SOLEDGE3X compared

## 2.5. Target particle and heat fluxes

Figure 6 shows the perpendicular axisymmetric flux density profiles at the targets for plasma particles (top row) and the total heat. The ‘total’ heat flux denoted $q_{\perp,\text{sym,tot}}$ consists of the net sum of all heat deposition channels to the
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Density and temperature at the inner (top row) and outer (bottom row) target plates, for six color-coded throughputs (in e⁻·s⁻¹) of the throughput scan. See figure 3 caption for the color coding.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Target perpendicular toroidally symmetric electron particle flux and total heat flux (including contributions from plasma, neutrals, volumetric radiation, and wall recombination, but excluding radiation from the core plasma that is not in the numerical domain) at inner (left) and outer (right) targets. See figure 3 caption for the color coding.</figcaption></figure>

wall: kinetic from ions, electrons, and neutrals, plus the radiation from volume processes inside the plasma and the latent heat release from recombination of incident ions at the wall. Plasma particle fluxes show very good agreement at the targets, except for the lowest density case. At high throughput, fluxes in SOLEDGE3X are slightly higher. The particle fluxes

for neutrals (not pictured) show more difference: atom and molecule fluxes are larger in SOLEDGE3X, and these differences increase with the throughput, up to +20% on atoms and a factor of 2 on molecules at high throughput.

Total power deposition profiles on the targets figure 6 show good agreement for both codes especially around the rollover.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Two left panels: neutral pressure in the #103 030 (6.80 × 10²² e⁻.s⁻¹, corresponding to the rightmost points for each code on the right panel) case for SOLPS (left) and SOLEDGE3X (middle), red cells are used for computing the average. Right panel: divertor neutral pressure as a function of the separatrix density for SOLPS and SOLEDGE3X, see in the text for the definition of this measure.</figcaption></figure>

At the lowest density, a 50% difference is found coming from the reversed asymmetry between the inner and outer side, which is a consequence of the significantly lower plasma temperatures in SOLEDGE3X at the outer target for this case.

At high throughput, power fluxes are higher by 30%–40% in SOLEDGE3X, but those values matter less because power fluxes are already strongly reduced in these partially detached conditions. The main difference was found to come from the relative weight of the plasma and the neutrals' contributions to the heat deposition at the inner target. While on the outer target, the agreement is good, the inner target plasma power deposition peaks in SOLPS-ITER are a factor 2 below those of SOLEDGE3X, even in attached cases. This difference is compensated by the heat flux from neutrals, so that the total power deposited remains similar. SOLPS-ITER plasma heat fluxes exhibit an inner/outer asymmetry which is less present in SOLEDGE3X. Neither SOLEDGE3X nor SOLPS runs included drifts for the cases considered here. The cause of this difference is currently under investigation, but it is worth noting that similar results have already been reported in the literature during comparisons between edge plasma codes (see [26] for a comparison between UEDGE, SOLEDGE2D, and SOLPS-ITER, figures 3(d) and (h)).

## 2.6. Divertor neutral pressure

In ITER cases run with SOLPS-ITER and mentioned in the literature, the quantity called divertor neutral pressure $p_{\text{div}}$ is computed as the average of the neutral pressure in cells at the edge of the SOLPS plasma domain between the target and the dome umbrella [2]. For comparability, the same average has been computed with SOLEDGE3X across the scan as shown in figure 7. On the right side is shown the relation between the midplane separatrix density and the divertor neutral pressure, showing a very good agreement on this behaviour, especially at low throughput. At high throughput, this pressure as a function of the separatrix density rises faster in SOLEDGE3X than in SOLPS-ITER, but the trend is very similar, even if differences exist between cases, because the response to throughput is not identical (especially at high throughput).

The two left plots show the 2D distribution of the neutral pressure in the divertor for both codes at the highest throughput case. Neutral pressures below the dome are very similar, even if minor differences may exist. Those can exist while still resulting in the same throughput because the pumping is modeled with a fixed particle absorption probability, and not a pumping speed. An increased penetration of atoms along the flux tubes just below the separatrix is found in SOLPS, which is not observed in SOLEDGE3X. Tests were carried out with different wall temperatures in SOLEDGE3X (not pictured here), but did not lead to any significant change in divertor neutral pressure. As a consequence, this parameter difference between the two codes can probably be excluded when attempting to explain the difference in $p_{\text{div}}$. The root cause for this difference is under investigation.

## 2.7. Take-aways from the comparison with SOLPS-ITER

In summary, even though some quantitative differences have been found between the two codes, a good match was obtained for a large part of the scan, and the key results are not changed. The global patterns, trends, and amplitudes are well-matched. Therefore, one can conclude that no contradiction is found with published results with SOLPS-ITER, thus strengthening confidence in the results obtained with both codes. Similar<fnref n="7" /> simulations are the basis for the detachment and plasma-neutral processes study in the next section, and the first wall erosion sensitivity analysis in the last section.

# 3. Divertor detachment and plasma-neutral processes study

## 3.1. Objectives and approach

The transition from the high-recycling regime to the detached regime is a multi-faceted phenomenon with many parameters,

<footnote n="7">⁷ The detachment study simulations are not strictly the same because they use a slightly different equilibrium and wall, even if the rest of the setup is identical (see simulation setup section). A quick comparison between the two geometries did not highlight significant differences.</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Left and center: Peak particle and total heat flux density from plasma, neutrals, and radiation at the outer and inner targets as a function of the separatrix density at the outer midplane across the throughput scan. Right: Evolution of the SOL width $\lambda_{q,e}$ with the separatrix density. $\lambda_{q,e}$ is defined as the exponential decay length of the electron parallel heat flux taken at the poloidal location where it is maximum. See definition in the text.</figcaption></figure>

and some of its driving mechanisms are still under discussion [27–32]. More specifically, while it is certain that the plasma-neutral interactions play a crucial role in the reduction of the particle and energy fluxes to the targets, there is no definite consensus yet on the role and importance of each underlying process in all the different types of plasma scenarios and machines. This work attempts to clarify what happens in the context of the ITER simulations considered here. This section focuses on plasma-neutral processes in the simulations, to build a high-level, schematic picture of the phenomena driving the onset of detachment here. Since it involves many other parameters that are not explored here, either because they depend on the machine (target materials, geometry, etc) or the pulse configuration (input power, magnetic equilibrium, fuel isotope, impurity seeding, etc), the results obtained here cannot claim to provide a picture encompassing all scenarios in all machines.

This exercise is based on a throughput scan, run now with the 2021-updated $dr_{sep} = 6.5$ cm magnetic equilibrium and wall contour, and is structured as follows:

(i) First, peak flux densities at the targets and midplane profiles are reviewed to assess the divertor regimes obtained across the scan.

(ii) Then the plasma-neutral interactions are decomposed into individual contributions from each one of the reactions, for three cases representative of the attached, roll-over, and partially detached conditions, in which the key contributors and trends are identified.

(iii) Then, from this analysis, a schematic view is constructed with the most important processes in each regime.

## 3.2. Divertor regimes

Across the throughput scan, a clear particle flux rollover at the targets is seen in left and center graphs of figure 8. It occurs simultaneously on the inner and outer sides, and the peak heat load strongly decreases.

One can note the strong degree of symmetry obtained across this scan. Here, even if the flux densities of particles and heat are higher on the inner side, both targets roll over at the same outer midplane density. One may have expected the inner side to roll over at lower midplane density as is observed in experiments and in SOLEDGE3X simulations of smaller machines, however determining the reason for this noticeable symmetry would require additional investigations outside the scope of this paper.

Before the roll over, the particle flux already starts to saturate with the separatrix density increase, which shows that the scan's first point is already well into the high-recycling regime. In the early stages of the conduction-limited regime, the particle flux is expected to scale with $n_u^2$ [27], which is not the case here. At the last point in the scan, target profiles still have significant plasma density and particle fluxes, but only in the regions further away from the strike point, as is characteristic of the partially detached regime. All quantities are strongly reduced with very low temperatures below 1 eV. In this case, the heat flux is mainly carried by neutrals and radiation.

Turning now to 1D profiles at the outer midplane, figure 3$^8$ shows the expected behavior for this type of transport simulations, i.e. close to exponential decays after a few centimeters from the separatrix. Notably, no density shoulder appears at all in these profiles, which could have been the case even with uniform perpendicular transport coefficients under certain circumstances (for example, if a strong neutral ionization source had been found in the far-SOL). This result is in line with modeling results obtained with SOLEDGE3X and similar set-up (2D transport with uniform coefficients, no drifts) on medium-sized tokamaks for which mean-field simulations

<footnote>$^8$ The plots show the $dr_{sep} = 11$ cm equilibrium for the sake of brevity, as the results on the $dr_{sep} = 6.5$ cm are almost identical. They can be found in [21].</footnote>
have not been able to exhibit density shoulders without manually increasing perpendicular transport coefficients in the far-SOL w.r.t. near SOL values (which are constrained by profile measurements).

While the electron temperature inside the separatrix decreases with throughput, which could have been reasonably expected, one can observe however that further in the SOL, the $T_e$ profile is almost insensitive to the throughput. This effect mainly comes from the more effective electron-ion energy equipartition with increased density and the large $T_i/T_e$ ratio in the far-SOL. Details of this topic are discussed in [21].

### 3.3. $\lambda_q$ increase with separatrix density

The SOL power width $\lambda_q$ is a central parameter to consider when addressing the heat exhaust challenge, as it is one of the determining factors of the peak heat flux arriving at the targets.

In numerical simulations, it is possible to compute it with a direct approach: $\lambda_q$ is evaluated via $\lambda_{q,e}$ the characteristic decay length of the electron parallel heat flux density at the poloidal location where it is maximum and then remapped to the outer midplane. This method, as it involves quantities upstream above the divertor, applies to any regime, not only to low densities, and is also used in other analyses based on the SOLPS code [33]. The results are shown in the right plot of figure 8.

The first case gives a value at 4 mm and is roughly consistent with the well-known multi-machine inter-ELM H-mode power width scaling law [14], which gives around 2 mm for a 5 MA hydrogen H-mode, and needs to be doubled for L-modes [3]. However, it is then found to strongly increase from 4 mm to 6 mm up to the roll-over, and then continue to increase to 6.5 mm, in a much slower fashion. This behavior was also observed in SOLPS-ITER runs [34]. The last large step to 12 mm is only due to the maximum throughput case, where the density rises so much that profiles tend to be significantly flatter.

Therefore, this shows that in our simulations, even with identical perpendicular transport coefficients, changes in the divertor regime have a noticeable effect on the shape of the SOL upstream. In such cases, what happens upstream is not independent anymore from what happens in the divertor. The origin of the change in shape of the power channel upstream can be attributed to two effects: the increase of the perpendicular heat transport with rising density, and the reduction of the parallel one via, first for the lower part of the scan, convective flow reversal, and then for the last partially detached case, conductivity reduction. A more detailed analysis of this matter is provided in [21]. However, it is recalled that those are low-power simulations, without impurities, and the latter observation may not hold for high-power impurity-seeded burning plasma cases as seems to be indicated by SOLPS simulations in [2].

To conclude on this matter, it is emphasized that running transport simulation density scan at fixed perpendicular transport coefficients does not mean running at fixed $\lambda_q$, because the latter increases with density. Thus, if one wishes to carry out scans at fixed $\lambda_q$, transport coefficients must be tuned accordingly.

While the increase of $\lambda_q$ with density found here may remind of the high-density $\lambda_q$ scaling law in [35], the dependency is not the same here because the suspected underlying physics of self-consistent turbulence transport is not included in the presented simulations.

### 3.4. Analysis of the plasma-neutral interaction processes leading to detachment

Detachment, in the sense of the reduction of the particle flux at the targets, stems from the increased interaction between plasma and neutrals, which generates particle, momentum, and energy sources and sinks. This section aims at identifying the key processes driving those three sources, in an attempt to build a high-level picture of plasma-neutral interactions in the considered ITER simulations. Notably, specific attention will be paid to the nature of the particle sinks, as well as their relative importance in the different regimes.

For this exercise, we investigate the contributions from each individual reaction/collision included in the plasma-neutral interaction model to the plasma source terms. The list of those processes can be found in table A2. In the past, the role of different plasma-neutral interaction process have been investigated with EIRENE (for example in [3, 29, 30, 36–38]), but only aggregated by category, i.e. bundled by type of neutral: atom, molecule, molecular ion. Here, for the first time, all of the individual contributions of all processes in the EIRENE model are decomposed, and analysed.

The present section uses plots in figure 9 of the volume sources along the first flux surface above the separatrix (see top left of the figure) in the outer divertor, for the last 10 cm poloidally next to the target. Processes are color coded, for three cases representative of each regime (coded by the type of line): one at the lowest density ($1.46 \times 10^{21} e^{-} s^{-1}$ in solid lines) for the attached plasma, one at rollover ($8.85 \times 10^{21} e^{-} s^{-1}$ in dashed lines), and one at the maximum throughput (at $1.76 \times 10^{23} e^{-} s^{-1}$ in dotted lines) for the partially-detached plasma.

The following summarizes what can be observed in the four sources plots figure 9. Descriptions of the abbreviations of processes can be found in table A2. For a more in-depth description and analysis, the reader is referred to [21].

* **In all regimes:** The strongest electron energy loss is the atom ionization, due to its high cost per event (~25 eV on average, increasing with $T_e$ = decreasing in cold regions, but decreasing with density, vs. 15 eV cost of ionizing a molecule).

* **Attached regime/high-recycling:** Atom ionization is the largest contributor to the ion particle source, followed by the already significant $H_2^+$ dissociation ($H_2^+$ EDE). Because $T_e$ is large enough, $H_2^+$ ions are created by electron ionization, and thus the $H_2^+$ dissociation is a net ion source. All plasma-neutral reactions happen directly at the targets. The electron energy loss is mainly composed of the expensive atom ionization, the molecule ionization, and the $H_2^+$ dissociation. The ion energy and momentum sources are only negative (sinks), driven by the atom charge exchange.
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Top-left plot: The thick black line is the location of the reaction profiles plotted on the right-side graphs. Right: volume ion particle source (center top), ion momentum source (top right), ion (center bottom) and electron (bottom right) energy source along the outer divertor leg. Plasma-neutral processes are color coded, and cases are coded with line styles (see legend).</figcaption></figure>

*   **At rollover:** At rollover, the atom ionization front detaches, while the $H_2^+$ dissociation increases and stays at the target. However, $T_e$ decreases, and $n_i$ increases enough so that the $H_2^+$ creation mechanism completely switches from electron impact to ion impact (i.e. through molecule charge transfer). Molecule electron ionization disappears, thus the $H_2^+$ dissociation becomes a net zero particle source. Molecules are now also significantly dissociated directly into two atoms further away from the plate. The rest of the $H_2^+$ are consumed by dissociative recombination (MAR), which becomes significant at the target as a particle sink. Electron-H$^+$ recombination remains very weak. Ion-molecule elastic scattering is now a large momentum sink, and the energy transported by the atom charge exchanges starts to heat the cold target plasma after taking energy upstream.

*   **Partially detached:** Finally, at maximum throughput, the atom ionization is largely spread out up to the X-point, and does not occur close to the strike point. The $H_2^+$ dissociation front detaches and spreads. The molecule dissociation into two atoms disappears. Electron-H$^+$ recombination becomes a large particle sink near the target, and the $H_2^+$ dissociative recombination (MAR) detaches from the target where the $H_2^+$ density is high (but still several orders of magnitudes below $n_i$). Electron energy now includes a small source directly at the plate from the three-body recombination. Strong

losses include atom ionization on the whole divertor leg, and the $H_2^+$ dissociation near the target. The atom charge transfer adds now momentum and energy in a larger region between the X-point and the target, while the molecule-ion elastic scattering makes a large and wide momentum sink at this location, which also results in pushing the molecules back to the target. Neutral-neutral collisions become large along the divertor leg at this stage.

In order to highlight two specific points about the investigations of the detachment mechanism, it is observed that:

(i) MAR is the particle sink that appears at the rollover while electron-H$^+$ recombination is not significant, and the latter becomes dominant on the field line only at detachment, and closer to the target.

(ii) The molecule-ion elastic collisions and atom charge exchange are each responsible for both very large momentum and energy sources, and these two types of sources are directly linked and cannot be addressed as independent effects.

As means of summary, a synthetic view of previous description is proposed in figure 10.
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Overview of detachment process with key reactions involved, for the three divertor regimes. The height of the boxes represent the spatial extension of where the processes are important, between the target (bottom) and the X-point (top). The ‘+’ and ‘−’ signs for the charge-exchange orange box indicate whether this process contributes positively to the plasma energy (energy source, plasma is heated by atoms) or negatively (energy sink, plasma is cooled by atoms). Note that the labels in color represent the different regimes, but the behaviors represented are derived from the analysis of a single field line. Thus the ‘partial detachment’ label relates to the regime, while the considered field line is itself simply ‘detached’. Arrows represent the $H_2^+$ creation-consumption reaction chain. See table A2 for the list of abbreviations.</figcaption></figure>

### 3.5. Relation between momentum and energy fractions

From this density regime scan, it is also possible to evaluate the fraction of plasma parallel momentum $f_{mom}$ and energy $f_{power}$ transmitted from the midplane SOL to the target.

Those fractions are defined as such:

$$f_{mom} = \frac{p_{i}^{tot}}{p_{u}^{tot}} \qquad (1)$$

$$p^{tot} = n(T_i + T_e)e + \frac{1}{2} \frac{\Gamma^2}{n} \qquad (2)$$

$$f_{power} = \frac{q_{\parallel,t}f_{x}}{q_{\parallel,u}} \qquad (3)$$

In the latter expression, $f_x$ is the flux expansion between midplane and the target plate.

The results for the throughput scan are plotted in figure 11. As is shown in the left graph, both $f_{mom}$ and $f_{power}$ decrease drastically with increasing midplane density as expected. More notably, when plotting $f_{power}$ against $f_{mom}$ (the right part of the figure), one can see that for high-density regimes near and after the rollover $f_{power}$ is proportional to the square of $f_{mom}$. This observation seems to support the point that, in the considered simulations, at high-density regimes both fractions are closely linked. Thus, they cannot be arbitrarily varied independently from each other in reduced models

to derive conclusions. This squared relation recalls the relation between momentum $p$ and energy $E \propto p^2$, illustrating an energy loss associated with momentum losses. However, more detailed investigations would be needed to better understand this observed link between the two fractions.

## 4. First wall gross erosion simulations

### 4.1. Motivation & approach

In ITER, estimating the First Wall erosion rate in the considered plasma scenarios constitutes therefore a crucial task, and those computations are based on the foreseen properties of the boundary plasma. Studies of this topic have been carried out to assess the erosion/reposition patterns with WALLDYN [16], and identification of potential constraints regarding excessive tungsten (W) erosion in certain magnetic configurations. Recent studies with the ERO2.0 code [11] also showed very large variations in the estimations of the erosion and transport of Beryllium (Be) to the divertor depending on the plasma conditions, magnetic configurations and fuel isotope, in the aim of evaluating dust formation, fuel retention, and identifying the wall tiles which are critical to the wall lifetime.

While previous simulations of the plasma reaching the First Wall in ITER relied on extending main-SOL SOLPS [1] simulations into the far-SOL using the OEDGE code suite [11, 12], SOLEDGE3X-EIRENE has the particularity to use a single numerical domain for the plasma, extending up to all

<footnote n="9">⁹ In the literature, the symbol $f_{power}$ may also refer to the $1 - f_{power}$ used here. The notation here is defined to be consistent between momentum and energy.</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Fraction of momentum and energy at the target with respect to the midplane SOL density (as defined in equation (1)), and energy fraction as a function of the momentum fraction for the throughput scan. The green arrow indicates the direction of increasing separatrix density.</figcaption></figure>

plasma-facing components. Therefore, the solution produced is spatially consistent (with respect to the considered physics model), and can be used to study the impact of the scenario on plasma quantities at the First Wall. This is the objective of this study, which aims at producing two sensitivity analyses for the First Wall gross erosion rate:

(i) In section 4.4, on the impact of the divertor regime from high-recycling to partial detachment.

(ii) In section 4.5, on the impact of enhanced perpendicular transport in the far-SOL, representing the presence of density shoulders.

While one of the goals of carrying the present simulations is to provide the fluxes and plasma backgrounds in the far-SOL and at the first wall to erosion-specialized codes [11, 16, 17] so they can assess the impurity generation, transport, and redeposition, it is possible to set up EIRENE in SOLEDGE3X so that it computes a 'gross' beryllium erosion rate (i.e. without describing the fate of the eroded particles) from physical sputtering (no chemical sputtering is taken into account) on each of the wall segments, using the Roth-Bohdansky-Eckstein model [39, 40]. It is, of course, less precise than what the specialized codes can do because of the absence here of 3D effects (toroidal geometry of panels, tile shaping), the use of a simplified wall geometry contour, the absence of beryllium self-sputtering or more precise Plasma-Wall Interaction database models for sputtering yields (e.g. taking into account Chemically-Assisted Physical Sputtering [41]). The motivation for using these estimations with SOLEDGE3X here is to make a sensitivity study regarding the beryllium sources, rather than precise computations of impurity fluxes and transport.

## 4.2. Erosion modeling setup

The gross erosion rates only contain two source mechanisms, whose contributions are investigated:

(i) Sputtering from $H^+$ plasma ions reaching the First Wall, accelerated in the electrostatic sheath.

(ii) Sputtering from fast $H$ neutrals, in this case high-energy $H$ atoms that have undergone charge exchange in hot plasma regions.

Computing this gross erosion rate is done as a post-treatment, as no feedback on the plasma is included, i.e. the sputtered particles of Beryllium are not tracked, no Beryllium plasma ions or neutral gas is included, and associated radiative losses have not been accounted for in the plasma solution.

## 4.3. On the choice of Beryllium as wall material in the simulations

In the presented simulations, the wall material was set up as Beryllium (Be), for historical reasons. Following the recent consideration for a change to a Tungsten (W) first wall in ITER, the role of Be vs W in the simulations is briefly discussed in this section.

The choice of material here only impacts the way hydrogen ions are recycled into neutrals<fnref n="10" /> at surfaces in the plasma-wall interaction model in EIRENE (share of back-scattered atoms vs. thermally re-emitted molecules, fast-reflection angles and energies), since the Beryllium impurity is not taken into account in the plasma model. A study of similar plasmas with the SOLPS-ITER code investigated changes between

<footnote n="10">¹⁰ Reflection properties of neutrals from surfaces are also affected, but this impact is assumed to be minor.</footnote>
Beryllium and Tungsten for the target materials [3], and showed that, indeed, no major difference was found in the key quantities of interest (plasma parameters, divertor neutral pressure), apart from the fraction of atoms vs. molecules in the neutral content. Therefore, it can also be expected that in the present cases, switching from a Be wall to a W wall would not result in large differences in the plasma solution. In fact, including the sputtered wall material impurity may entail larger changes in the SOL plasma solution with a Be wall than with a W wall, due to the tendency of Be impurities to radiate predominantly in the SOL, while W radiation would rather be located in the core region which is not included in the numerical model here.

The gross erosion rate, as computed here, is expected to differ significantly between a W wall and a Be wall. However, since only a gross rate is considered here, which depends on the same drivers (particle fluxes, energy of particles) whichever the wall material, the trends and variations observed here, which are the focus of this study, should remain true with a W wall, while absolute values would be strongly scaled down.

Lastly, one may also consider Be here as a proxy for low-Z elements such as boron (B) ($Z_{Be} = 4$, $Z_{B} = 5$), which may be used as a part of the ITER wall conditioning strategy [42, 43] in a full-W environment.

Running ITER simulations with a W wall material is currently a work in progress and will be the focus of a subsequent paper.

### 4.4. Impact of density/divertor regime

In this section the impacts of increasing the separatrix density are studied, with radially uniform transport coefficients. The throughput is chosen as the parameter scanned, with seven simulations ranging from $1.46 \times 10^{21}$ at low density, in attached but already high-recycling conditions (in the sense that the target particle flux does not follow the $n_{sep}^2$ relation anymore, but starts to saturate), up to $1.76 \times 10^{23}$ e s$^{-1}$ at partially detached conditions (i.e. the particle flux density and temperatures are strongly reduced, see figure 8), and the majority of the SOL radiation is still located near the targets. The rollover of the particle flux is located at $8.85 \times 10^{21}$.

#### 4.4.1. Plasma conditions at the first wall.
Figure 12 shows the effect of increasing the midplane density on ion particle flux density (perpendicular, toroidally symmetric), electron and ion temperatures, plotted along the first wall using a curvilinear coordinate running from the outer side to the inner side. The effect on plasma quantities is relatively mild: the increase in throughput translates to an increase by a factor of 10 on the ion particle flux at the machine top, but the values remain still very low ($\sim 1 \times 10^{19}$ s$^{-1}$m$^{-2}$ vs. $\sim 1 \times 10^{20}$ s$^{-1}$m$^{-2}$). $T_e$ is almost unchanged at around a few eV, except at the top, where it drops from 10 to 1 eV. This, along with the high plasma

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Impact of separatrix density (increased with throughput at values in legend) on plasma ion particle flux density (perpendicular, toroidally symmetric) electron and ion temperatures and incident atom particle flux density (perpendicular, toroidally symmetric) at the First Wall as a function of the position along the wall.</figcaption></figure>

density observed there, would indicate that the secondary X-point is also undergoing a sort of detachment process, similarly to a divertor. The ion temperature sees a more significant effect and decreases from 10–20 eV to a few eV. The effect of density is much more visible on neutral particles, as seen on the incident atom fluxes (bottom graph of figure 14). They increase drastically with density, gaining up two orders of magnitude from $\sim 1 \times 10^{19}$ s$^{-1}$m$^{-2}$ to $\sim 1 \times 10^{21}$ s$^{-1}$m$^{-2}$ at the top. Incident molecule fluxes (not pictured) follow the same trend but are less relevant for erosion considerations due to their low energy.

#### 4.4.2. Gross erosion estimations.
Figure 13 summarises the results obtained for the gross erosion rate, integrated over all First Wall segments. On the left is shown the contributions from the neutrals (CX atoms) in blue and from the ions in red, with the total plotted in black. The total gross sputtering rate increases with density by a factor of 6 between the scan extremes, with $\sim 90\%$ coming from the CX atoms. After the roll-over at $\sim 1.35 \times 10^{19}$ m$^{-3}$ the CX atoms contribution completely dominates and the ion contribution is almost extinguished. Indeed, even if the number of ions arriving at the first wall increases slightly with density, the ion temperature decreases, especially at the machine top ($T_e$ also decreases
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Impact of throughput on the gross erosion rate computed with EIRENE. Left: Total integrated rate across the first wall, with the ion contribution in red, and atoms contribution in blue. Right: 2D representation of gross erosion rate densities along the first wall contour.</figcaption></figure>

there). Those particle fluxes remain extremely low across the whole scan, with the major part of the First Wall receiving less than $10^{18} \text{s}^{-1} \cdot \text{m}^{-2}$, with local increases at $10^{20} \text{s}^{-1} \cdot \text{m}^{-2}$ above the outer baffle and at the machine top. As a consequence, the ion contribution to the sputtering is almost independent of density.

Figure 13 also shows on the right the areas where sputtering is most important. For the case with the highest obtained gross erosion rate (highest density), those are the puff location at the top, and the tile above the outer baffle. In reality, the puff location consists of a port and not a flat tile, therefore 3D or more precise wall geometry effects may be important there, which are not captured with the simplified wall contour used in our simulations. Also, the particle injection via the puff is assumed axisymmetric here, and not localized toroidally.

### 4.4.3. Fluxes & average energy of charge exchange atoms.

CX atoms being the largest contributors, the focus is turned to them in this section. Figure 14(bottom) shows the atom incident particle fluxes obtained on the first wall, along with their average energy. The atom particle flux increases, as expected with increasing gas injection, especially at the machine top and in the vicinity of the upper outer tiles (wall coordinate ~9 m), where it increases by a factor of 100 across the scan.

The average energy decreases drastically however, as shown on the top plot of figure 14. It changes by a factor of 10 to 100 (depending on the location along the wall) from 200–500 eV down to 20–30 eV and even as low as 3 eV at the machine top (where the particle flux is the largest). A more detailed discussion and analysis of the origin of those CX

atoms reaching the FW with a simple reduced model (shown with stars symbols in figure 14) can be found in [21].

### 4.4.4. Energy distribution of incident atoms.
The energy distribution of the atoms incident to two sections of the first wall, the lower outer wall and the machine top (in columns), are shown in figure 15, for the three cases representative of the throughput scan (in rows): attached, rollover, and partially detached. It is readily seen that two different populations exist, in the same way as was measured in experiments [44, 45]:

(i) A 3–4 eV ‘cold’ population which consists of two almost mono-energetic peaks at 2.7 eV and 4.0 eV. These populations come from the dissociation of molecules and $H_2^+$ molecular ions, which release atoms at fixed Franck-Condon energies.

(ii) A hot ‘CX’ population, ranging from tens of eV to hundreds of eV, and even up to 2 keV in some cases. These are the results of charge exchange reactions in the much hotter regions of the plasma, often behind the separatrix and from ion distribution tails.

Furthermore, the flux-weighted average energies (vertical black dashed lines and values) show that the estimates in the previous 1D graph figure 14 were quite representative of the energy of the charge exchange population, because of its much higher energy. The CX atoms are almost always the majority population (between 50% and 75% of the incident particle flux), except for the machine top at maximum density.

Notably, the obtained CX distribution shape is not exactly a Maxwellian distribution (in dashed green in the graphs).
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Impact of density (with increasing throughput at values in legend) on FW atoms. Top: poloidal distribution of the average energy of incident atoms along the first wall contour (travelled counter-clockwise). Star symbols represent the approximation computed from the reduced model in [21]. Bottom: poloidal distribution of the atom toroidally symmetric particle flux normal to the wall along the first wall contour. The curvilinear coordinate runs along the first wall contour, starting from above the outer baffle, and ending above the inner baffle.</figcaption></figure>

Such a distribution could have been expected because ions are assumed to follow a Maxwellian distribution in the fluid plasma solver, and are sampled as such (then convoluted with a weighting from the cross section $\sigma^{Atom\ HCX}$ which decreases at high energies) in the Monte-Carlo procedure. However, because fast CX atoms can have mean free paths larger than ion temperature gradient lengths, this distribution is generated from interactions with Maxwellian ions at varying temperatures on different flux surfaces. In such cases the resulting distribution is not a single-temperature Maxwellian anymore. Additionally, the obtained CX distribution is missing the lower-end of the Maxwellian distribution, and is closer to a Gaussian law in logarithmic space (red dashed lines). This absence can be interpreted as stemming from the additional collisions in the SOL which are more probable for lower energy atoms and thus would filter them out more effectively than high energy ones. Further, the CX location vary with the penetration depth of the atoms so the sampling of CX neutrals actually stems from a superposition of (cross-section filtered) Maxwellians at various temperatures. From the Gaussian fit, the high bound of the energy to capture 95% (i.e. within $2\sigma$) of the flux is evaluated (red arrows) to assess the most energetic atoms in the distribution. Very high values above 1 keV are found for the attached case because the plasma density is much lower, and as a consequence, atoms can penetrate further into the main plasma. However, these fluxes are quite low.

Apart from some exceptions at maximum density, a general trend seems to appear: an increased density does not seem to change much the shape of the distributions (i.e. the ratio of CX to cold population, and the width of the CX distribution), whereas it does shift the center of the CX distribution toward lower energies. Such decrease in energies could indicate that

sputtering in a case with a W FW, which has a higher threshold energy, may strongly reduced. Such investigations are currently in progress, and the first elements are reported in [46].

### 4.5. Impact of density shoulders

#### 4.5.1. Motivation.
This section considers the impact of increasing the anomalous perpendicular diffusion coefficients in the far SOL. This exercise follows the many reports made on the formation of density shoulders in the SOL since the 1990's [47], and under high divertor density regimes [48, 49]. In such conditions, plasma is strongly transported towards the main wall in the far SOL, increasing plasma-wall interactions there. Although the exact conditions in which SOL density shoulders appear and their quantitative impact on plasma transport properties remain open questions for ITER, the current understanding suggests that they might be present in ITER high-density regimes. This calls for a sensitivity study of erosion rates as a function of far SOL transverse transport coefficients, leveraging again the advantage of the ability to self-consistently simulate the plasma (modeling turbulent transport using a diffusive approach) at this location with SOLEDGE3X. It is the aim of the present section, which builds upon the initial plasma conditions estimations carried out in [7].

#### 4.5.2. Modeling density shoulders.
The understanding of the formation and characteristics of density shoulders is the topic of active research. The flattening of the density profile is currently understood as being related to convective LFS filamentary transport, and several proposed mechanisms are invoked in recent literature [48, 49]. Those works study the
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Energy distribution of atoms incident on two of the first wall sections, for three cases (in rows) representing the attached, rollover, and partial detachment state. The black vertical line is the total average energy. Two fits are carried out: for a Maxwellian distribution ($f_E(E) = 2\sqrt{\frac{E}{\pi}}(\frac{1}{T})^{3/2}\exp(-\frac{E}{T})$ with T and E in eV) in green with moderate success, and a more successful fit with a Gaussian law in logarithmic space in red. In each panel, the vertical red dashed line represents the center of the fit, and the red arrow shows the energy upper bound required to capture 95% of the distribution (i.e. $2\sigma$).</figcaption></figure>

correlations between the amplitude of the density shoulders and other parameters: divertor collisionnality, recycling in the far SOL, divertor regime, radial electric fields, particle source in the main chamber, plasma current, safety factor, and midplane neutral pressure. Experimentally, the appearance of density shoulders is correlated with the transition to high-density regimes [48, 49]. However, all across the throughput scan in our 2D transport simulations with uniform and constant prescribed transport coefficients, which extends after the rollover into partial detachment, no sign of shoulder is observed (see the midplane density profiles in figure 3).

This implies two possibilities: either a density shoulder would not form in these situations in ITER (this is not known), or the simulation lacks a necessary ingredient to form such shoulders. The latter is consistent with the interpretation of density shoulders being driven by a local change in the amplitude of transverse turbulent transport as previous simulations were obtained at fixed transverse transport coefficients. A fully consistent modelling would require a self-consistent description of turbulent transport in the SOL, which is to date not achievable for ITER scale simulations in realistic parameters with the relevant divertor physics. We can nevertheless
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Outer midplane radial profile of far SOL perpendicular diffusivity coefficients for the far SOL diffusivity scan: coefficients are increased to higher values (from 0.3 m² s⁻¹ to either 2.0 m² s⁻¹ or 10.0 m² s⁻¹) from 1 or 3 cm after the separatrix at the outer midplane. These coefficient radial profiles are then propagated along flux surfaces.</figcaption></figure>

approach the problem by scanning artificially transverse transport coefficients in our mean-field simulations to assess the sensitivity of the first wall fluxes to the possible formation of a shoulder.

Hereafter, the enhanced transport giving rise to density shoulders is modeled by increasing the perpendicular diffusivity coefficients in the far SOL from the base value of $D_\perp = 0.3 \text{ m}^2 \text{s}^{-1}$ (which was used uniformly in the throughput scan). Two parameters are varied:

* the value in the far SOL, either $D_\perp^{\text{far-SOL}} = 2.0 \text{ m}^2 \text{s}^{-1}$ or $10.0 \text{ m}^2 \text{s}^{-1}$

* the distance from the separatrix at the outer midplane where this increase starts: either 1 cm or 3 cm (translating respectively into $1.8\lambda_q$ vs. $5.4\lambda_q$, with a $\lambda_q$ of 5.5 mm for the reference case).

The increase from the base value to the far SOL value is linear over a 1 cm ramp in all cases. These midplane profiles are then propagated poloidally along the flux surfaces across the numerical domain<fnref n="11" />, including in the divertor below the X-point. Other diffusion coefficients ($\chi_e, \chi_i, \nu_\perp$) are scaled proportionally to the $D_\perp$ increase from its base value of $0.3 \text{ m}^2 \text{s}^{-1}$, as shown in figure 16. The rationale for choosing such values comes from density profiles observed in machines such as JET, where the density decay length at the outer midplane $\lambda_n$ increases by a factor of 3 to 4 [49], corresponding to factors of 9 and 16 on $D_\perp$. Similar factors for the cases considered here would give 3 to 5 $\text{m}^2 \text{s}^{-1}$, however to include additional uncertainties, this range is extended to 2 to 10 $\text{m}^2 \text{s}^{-1}$. The simulations here could be therefore understood as

'best' and 'worst' cases in terms of diffusion where shoulders form. In total, 4 cases are considered and compared to the reference case #103027 with constant uniform coefficients ($D_\perp = 0.3 \text{ m}^2 \text{s}^{-1}$) at a low throughput of $8.85 \times 10^{21} \text{ e}^- \text{ s}^{-1}$.

Of course, these simulations still rely on a diffusive model in the shoulder even though this is not necessarily a faithful description of the perpendicular transport there, which may be more ballistic ($\Gamma_\perp \sim n\bar{v}_{\text{eff},\perp}$). However, the choice of closure for transverse fluxes matters little for the present study, as long as the chosen transport coefficients induce the desired gradient lengths.

A supplementary caveat is to be noted here: the choice was made here to model the density shoulders by increasing also the perpendicular heat transport coefficients, which translates to a shoulder also in the temperature profiles. This may not be a standard assumption as the presence of density shoulders in experiments is not always associated with a temperature shoulder. Applying also a temperature shoulder was rather done in the spirit of a sensitivity study in mind, to model 'pessimistic' cases.

4.5.3. *Impact on plasma profiles.* Figure 17 shows profiles obtained at the outer midplane, where the temperature and density profiles flatten as expected in the far SOL. This changes dramatically the values of the density and temperatures at the outer wall, by a factor of up to 1000 for the density, and a factor of up to 10 for temperatures (1 to 10 eV for $T_e$ and 3 to 30 eV for $T_i$). Inside the separatrix, the stronger the far SOL transport (closer shoulder, higher values of coefficients), the more the density rises (up to +50%) the lower temperatures get (up to -25%). A potential contributor to this is the main wall recycling (most significantly from the lower outer wall), which is found to drastically increase the particle source inside the separatrix, by up to a factor of 60 with

<footnote n="11">¹¹ This enhanced transport is therefore also propagated to the high field side, which may lead to overestimating perpendicular fluxes there, as anomalous transport is usually found to balloon toward the low field side in experiments.</footnote>
<figure><img src="image_20.png" /><figcaption>**Figure 17.** Impact of increasing far SOL diffusion coefficients on outer midplane conditions (plasma density, electron temperature, and ion temperature).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 18.** Impact of enhanced transport on first wall erosion. Left: gross erosion rate split between ions and atoms contributions for all the considered cases with increased far SOL transport. Right: 2D plot of the gross erosion rate for the most diffusive case. The arrow illustrates the factor of increase w.r.t. the reference case.</figcaption></figure>

respect to the prescribed core influx. This effect is discussed in [21].

A notable fact is that even though profiles change inside and outside the separatrix, values at the separatrix for $n_e$, $T_e$ and $T_i$ remain almost unchanged (less than ~10% difference) over the range of diffusivities considered.

In the divertor, the shoulder start location seems to be the most determining factor, as results using enhanced transport from 3 cm after the separatrix (red and green curves) are almost indistinguishable from the reference case in black. In all cases, little change is observed on particles (both density and particle flux). The energy (total heat flux and temperature) is more sensitive and peak values decrease by ~35% for the higher diffusive case. Notably, this decrease is seen near the strike points, where the diffusion coefficients are yet not increased, but only for cases where the diffusion enhancement starts at the midplane within less than $2\lambda_q$ ( = 1.2 cm, see figure 8).

The hierarchy of plasma-neutral interactions in the outer divertor leg remains unchanged with increased diffusivity, though with slightly scaled down amplitudes.

#### 4.5.4. Impact on gross erosion.
As was done in the throughput section 4.4.2, the sensitivity of the computed gross erosion rate to the enhanced transport is studied.

Figure 18(left panel) shows the total integrated gross erosion rate from the first wall, with the ion contribution in red and the atom contribution in blue. On the right is illustrated the location in 2D of the highest gross erosion rate density in each part of the wall. Contrary to the throughput scan with uniform diffusion, sputtering from ions is the dominant contribution (50%–70%) in cases with enhanced far SOL transport. This increase is much larger than the one found in the throughput
<figure><img src="image_22.png" /><figcaption>**Figure 19.** Impact of enhanced far SOL transport (shoulders) on plasma ion flux (top), electron (middle) and ion (bottom) temperatures at the first wall as a function of the position along the wall contour.</figcaption></figure>

scan (factor 6 increase), with a factor of 40 between the scan's extremes.

Looking at the spatial distribution, one finds like in the throughput scan that the top and outer lower tiles (#8,#9, & #18) show the highest sputtering rate.

### 4.5.5. Impact of the presence of shoulders on plasma fluxes.
To interpret the strong increase in the ion erosion contribution, figure 19 shows the influence of increased diffusion coefficients on the ion plasma flux and electron and ion temperatures at the first wall. From this increased transport, the ion particle fluxes are dramatically increased by factors of 100 or 1000, to reach $2 \times 10^{21} \text{m}^{-2} \text{s}^{-1}$ at the machine top and above the outer baffle. Temperatures are also further raised from the increased heat perpendicular transport coefficients to a few eV for electrons and to 10–20 eV for ions.

### 4.5.6. Impact of the presence of shoulders on fluxes & average energy of atoms.
Figure 20 shows, as was done for the throughput scan, the atom particle fluxes across the first wall, and the associated average energy. A significant increase in the incident atom flux is seen with increased far SOL transport, all along the first wall contour, by factors of 10 to 20.

Nevertheless, the average energy of these atoms at the first wall is almost completely unchanged. Investigations of the CX rate profiles at the midplane in [21] indicates that most of their collisions occur within 5 cm of the separatrix. In the previous section, it was seen that the presence of shoulders almost

does not affect the temperatures around the separatrix, and at 5 cm in the SOL, profiles have not yet changed dramatically. This way, neither the origin location of high energy charge-exchange atoms nor the plasma conditions there have significantly changed, which likely explains why the average energy of atoms at the first wall barely changes.

However, neutral densities vastly increase with diffusion, as is shown in the left of figure 21. Molecule density at the wall is raised by a factor of 100 between the scan's extremes, while atom density gains a factor of 10 at the separatrix. The reaction plots, in the middle and right plots of the figure, also show an increase up to a factor of 100 of the interactions as a consequence. This illustrates the increase in main chamber recycling, and the appearance of molecular processes in the far SOL, such as molecular ionization and dissociation, shows that the hierarchy of the plasma-neutral reactions near the first wall also changes with enhanced transport. On the right plot, in the far SOL, many charge exchanges occur, which was not the case in the reference simulation, indicating that atoms at lower energies (i.e. with a shorter mean free path) are getting increasingly screened by the colder far SOL plasma.

### 4.5.7. Energy distribution of incident atoms.
As for the throughput scan, the energy distribution of the atoms incident on two sections of the First Wall are shown in figure 22, for the two extreme cases in rows: the reference case, and the highest diffusivity case.

The following effect is observed: on the lower outer wall, the shape of the distribution does not vary and is simply scaled
<figure><img src="image_23.png" /><figcaption>**Figure 20.** Impact of enhanced far SOL transport on FW incident atoms. Top: poloidal distribution of the average energy of incident atoms along the first wall. Star symbols represent the approximation computed from the reduced model in [21]. Bottom: poloidal distribution of the atom toroidally symmetric normal particle flux along the first wall contour. The curvilinear coordinate spans across the first wall, starting from above the outer baffle, and ending above the inner baffle.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Figure 21.** Left: Atom, molecule, and total neutral density profiles at the outer midplane for the enhanced far SOL transport scan. Middle and right: reactions rates at the outer midplane for the reference case (flat diffusion coefficient profiles) and the most diffusive case (shoulder at 1 cm of the separatrix, $D_\perp = 10\,\text{m}^2\,\text{s}^{-1}$ in the far SOL). Note: the spikes in the far SOL are due to statistical noise, reactions are rare in this location without shoulders.</figcaption></figure>

up in amplitude, while at the machine top, two effects are opposed: the share of the CX population increases with diffusivity, but its width narrows, and its center energy is lowered. The global effect is a reduction by a factor two of the average energy (vertical dashed lines), as was observed in figure 20(at wall coordinate ~ 11 m).

To conclude this gross erosion section, it was seen that when increasing throughput with identical transport coefficients, two effects balanced each other: neutral particle fluxes were increased, but the energy of atoms which were inducing erosion was decreased (the contribution from ions was low). Now, with enhanced plasma far SOL transport, both ion particle and energy transports increase, and the two effects

accumulate. Ion-related sputtering increases drastically and now becomes dominant. The neutral-induced sputtering is raised as well, but not as strongly as the ions.

### 4.6. $T_i/T_e$ ratio

In this last part, the focus is turned to the ratio of ion to electron temperature in the SOL, and especially in the far SOL at the wall, as it is a key assumption for the erosion-specialised codes. The following paragraphs can be used to compare the assumptions which were used in erosion simulations versus what is obtained here from the up-to-the-wall plasma simulations.
<figure><img src="image_25.png" /><figcaption>**Figure 22.** Energy distribution of atoms incident on two of the first wall sections, for two cases (in rows) representing the reference case at rollover with uniform transport coefficients, and the most diffusive case with a shoulder at 1 cm of the separatrix, $D_{\perp} = 10 \text{ m}^2 \text{ s}^{-1}$ in the far SOL. The black vertical line is the total average energy. Two fits are carried out: from a Maxwellian distribution in green with moderate success, and a more successful fit with a Gaussian law in logarithmic space in red. The vertical red dashed line represents the center of this fit, and the red arrow shows the energy upper bound required to capture 95% of the distribution (i.e. $2\sigma$).</figcaption></figure>

Figure 23 shows the outer midplane profile for this quantity, for 3 cases to try to cover extremes for the parameters explored in this study. A common assumption is to take $T_i/T_e = 2$ by default in the SOL. The plot shows that for models with uniform diffusion coefficients (high and low throughput cases on the graph), this assumption is higher than what is found in the simulation in the near SOL (within 5 cm of the separatrix) which is around 1.2 to 1.6, decreasing with throughput.

In the far SOL however, the ratio largely increases going towards the wall, and the value at the first wall is actually between 2 for maximum throughput and 5 at the lowest throughput (the first wall is at 15 cm from the separatrix).

Because those field lines are sheath-limited and thus driven by the sheath (and not by the conductivity), this behavior comes from the difference in the sheath heat transfer coefficients between the electrons and the ions ($\gamma_e = 4.5$, $\gamma_i = 2.5$), which makes the electrons lose more energy at the wall than the ions under the same conditions. Because of this, the difference between the electron and ion temperatures will tend to increase radially, reaching a maximum at the wall. This is amplified by the fact that the collisional energy equipartition between the two species tends to decrease as the density decreases in the far SOL, so that the energy transfer from ions to electrons becomes less effective far from the separatrix.
<figure><img src="image_26.png" /><figcaption>**Figure 23.** Left: $T_i/T_e$ ratios at the outer midplane for three cases: low throughput (which corresponds to rollover at $8.85 \times 10^{21} \text{e}^-\text{s}^{-1}$), high throughput at $1.76 \times 10^{23} \text{e}^-\text{s}^{-1}$, and high diffusivity with strongly enhanced far-SOL transport (shoulder @ 1 cm, $D_\perp = 10 \text{ m}^2 \text{s}^{-1}$ in far-SOL with low throughput). Right side: the associated plasma density profiles.</figcaption></figure>

In terms of shoulder formation, not much changes in the near SOL, however in the far SOL and at the first wall, the ratio tends to reduce towards the values observed for large throughput. This is due to the increase of the density which enhances the collisional energy exchange between the two species and tends to bring the ratio towards unity.

In conclusion, estimates for the $T_i/T_e$ ratio could be provided as follows, depending on a macro ‘density’ indicator of the SOL with respect to the machines’ window of operation:

(i) In the near SOL: 1.6 to 1.2, decreasing with density.

(ii) In the far SOL: 2 to 5, decreasing with density.

## 5. Conclusion

In this work, 2D transport simulations of the ITER first non-active phase at 20 MW were studied. A throughput scan was carried out with the SOLEDGE3X-EIRENE code to describe plasmas scenarios in the vessel up to the first wall, ranging from high-recycling to partial detachment divertor conditions. Due to the specificity of the ITER cases when compared to current smaller machines (mainly related to the size of the machine), improvements were made to the code, both in terms of physics and numerical scheme. Most notably, the plasma-neutral interaction model setup in EIRENE was improved to better describe detachment, and the treatment of collisional terms in the plasma-solver was rendered more stable. A comparison of results with those obtained with the well-established SOLPS-ITER code (without grid extension) was carried out,

and a good agreement was found between the two codes. The evolution of the SOL plasma was analyzed. Key observations at the midplane included the expected separatrix density saturation with throughput, and an increase of the SOL power decay length $\lambda_q$ with throughput. Then, an analysis of the breakdown of the plasma-neutral interaction processes in the divertor led to the construction of a synthetic picture describing which reactions in the plasma-neutral interaction model are important in each regime and in each location in the divertor. This highlighted the role of the $H_2^+$ ion as an intermediary, and the strong link between the momentum loss fraction and the power loss fraction, showing that these parameters cannot be varied as independent variables in reduced models for high density regimes.

Two sensitivity studies have been carried out regarding the gross erosion rates for the beryllium first wall: the impact of throughput, and the presence of density shoulders in the far SOL. With increasing throughput, the gross erosion increases moderately, mainly driven by CX neutrals. Their incident particle flux increases, but is counteracted by their decreasing energy. With enhanced far SOL transport, the gross erosion increases drastically, now mainly from ions, due to the combined increased particle flux and ion energy.

The atom energy distribution was also investigated, in which it was found that two populations exist, one ‘cold’ (Frank-Condon) at 3–4 eV, and one hot ‘CX’ population. This hot population can have a very wide energy distribution, especially for low-density plasmas. With increasing throughput, this hot population energy is very strongly cooled, and can
even disappear collisionnally at the machine top. As the atom charge exchange reaction rate is almost independent of $T_i$, the plasma density is the main parameter for their localization. They are predominantly born from regions of plasma density at around ~ $1 \times 10^{19} \text{m}^{-3}$, which is around the separatrix. Therefore, since increased far SOL transport does not change separatrix conditions, the charge exchange atom energy at the wall does not change with the presence of density shoulders.

## Acknowledgment

The authors would like to thank the referees for their feedback and suggestions to improve this paper.

This work was performed under the auspices of the ITER Scientist Fellow Network and received funding from an ITER Organization Implementing Agreement (No. IO/IA/20/4300002306). The views and opinions expressed herein do not necessarily reflect those of the ITER organization. This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No. 101052200 - EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.
# Appendix

<table><caption>**Table A1.** Summary of simulation parameters.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>Particle perpendicular diffusivity D<sub>⊥</sub></td><td>0.3 m<sup>2</sup>.s<sup>-1</sup></td></tr><tr><td>Perpendicular viscosity ν<sub>⊥</sub></td><td>0.2 m<sup>2</sup>.s<sup>-1</sup></td></tr><tr><td>Perpendicular electron heat conductivity χ<sub>⊥,e</sub></td><td>1.0 m<sup>2</sup>.s<sup>-1</sup></td></tr><tr><td>Perpendicular ion heat conductivity χ<sub>⊥,i</sub></td><td>1.0 m<sup>2</sup>.s<sup>-1</sup></td></tr><tr><td>Parallel viscosity</td><td>ON</td></tr><tr><td>Conduction heat flux limiter coeff. for electrons</td><td>0.2</td></tr><tr><td>Conduction heat flux limiter coeff. for ions</td><td>OFF</td></tr><tr><td>Viscosity heat flux limiter coeff. for ions</td><td>0.375</td></tr><tr><td>Neutral model</td><td>Kinetic (EIRENE)</td></tr><tr><td>Fluid drifts</td><td>OFF</td></tr><tr><td>Collisional closure</td><td>Zhdanov</td></tr><tr><td>Turbulence</td><td>OFF</td></tr><tr><td>3D</td><td>OFF</td></tr><tr><td>Content</td><td>Pure H</td></tr><tr><td>Input power</td><td>20 MW (10MW e<sup>-</sup> + 10MW H<sup>+</sup>)</td></tr><tr><td>Particle input from the core</td><td>6.13 × 10<sup>20</sup> s<sup>-1</sup> + neutrals reaching the core</td></tr><tr><td>Sheath heat transmission factors</td><td>e<sup>-</sup> : 4.5 ; H<sup>+</sup> : 2.5</td></tr><tr><td>Target, first wall, dome material</td><td>Beryllium</td></tr><tr><td>Target, first wall, dome material temperature</td><td>500 K (uniform)</td></tr><tr><td>Target, first wall, dome recycling coefficient</td><td>1.0</td></tr><tr><td>Dome pillars material</td><td>Iron</td></tr><tr><td>Dome pillars reflection probability</td><td>0.5</td></tr><tr><td>Pump material</td><td>Iron</td></tr><tr><td>Pump recycling coefficient</td><td>0.9928</td></tr><tr><td>Gas injection location</td><td>Top outer port</td></tr></tbody></table>

<table><caption>**Table A2.** Summary of reactions included in the EIRENE setup (also referred to as the ‘Kotov model’).</caption><thead><tr><th>Reaction<sup>a</sup></th><th>IAEA reaction classification</th><th>Name</th><th>EIRENE identifiers<sup>b</sup></th></tr></thead><tbody><tr><td>H + e<sup>-</sup> → H<sup>+</sup> + 2e<sup>-</sup></td><td>Ionization (EIN)</td><td>Atom EIN</td><td>AMJUEL H.4,10 2.1.5</td></tr><tr><td>H + e<sup>-</sup> → H<sup>+</sup> + e<sup>-</sup> → H + e<sup>-</sup> + hν</td><td>Excitation (EEX)</td><td>included in above</td><td>included in above</td></tr><tr><td>H<sup>+</sup> + e<sup>-</sup> → H + hν</td><td>Radiative Recombination (ERR)</td><td>Ion ERR/T</td><td>AMJUEL H.4,10 2.1.8</td></tr><tr><td>H<sup>+</sup> + 2e<sup>-</sup> → H + e<sup>-</sup></td><td>Three-Body Recombination (ERT)</td><td>included in above</td><td>included in above</td></tr><tr><td>H + H<sup>+</sup> → H<sup>+</sup> + H</td><td>Charge Transfer (HCX)</td><td>Atom HCX</td><td>HYDHEL H.1,3 3.1.8</td></tr><tr><td>H<sub>2</sub> + e<sup>-</sup> → H<sub>2</sub><sup>+</sup> + 2e<sup>-</sup></td><td>Ionization (EIN)</td><td>Molecule EIN</td><td>AMJUEL H.4 2.2.9</td></tr><tr><td>H<sub>2</sub> + e<sup>-</sup> → 2H + e<sup>-</sup></td><td>Dissociation (EDS)</td><td>Molecule EDS</td><td>AMJUEL H.4 2.2.5 g</td></tr><tr><td>H<sub>2</sub> + e<sup>-</sup> → H + H<sup>+</sup> + 2e<sup>-</sup></td><td>Dissociative Ionization (EDI)</td><td>Molecule EDI</td><td>AMJUEL H.4 2.2.10</td></tr><tr><td>H<sub>2</sub> + H<sup>+</sup> → H<sub>2</sub> + H<sup>+</sup></td><td>Elastic Scattering (HES)</td><td>Molecule-Ion HES</td><td>AMJUEL H.0,1,3 0.3 T</td></tr><tr><td>H<sub>2</sub> + H<sup>+</sup> → H<sub>2</sub><sup>+</sup> + H</td><td>Charge Transfer (HCX)</td><td>Molecule HCX</td><td>AMJUEL H.2 3.2.3</td></tr><tr><td>H<sub>2</sub><sup>+</sup> + e<sup>-</sup> → 2H + hν</td><td>Dissociative Recombination (EDR)</td><td>H<sub>2</sub><sup>+</sup> EDR</td><td>AMJUEL H.4,8 2.2.14</td></tr><tr><td>H<sub>2</sub><sup>+</sup> + e<sup>-</sup> → H + H<sup>+</sup> + e<sup>-</sup></td><td>Dissociative Excitation<sup>c</sup> (EDE)</td><td>H<sub>2</sub><sup>+</sup> EDE<sup>c</sup></td><td>AMJUEL H.4 2.2.12</td></tr><tr><td>H<sub>2</sub><sup>+</sup> + e<sup>-</sup> → 2H<sup>+</sup> + 2e<sup>-</sup></td><td>Dissociative Ionization (EDI)</td><td>H<sub>2</sub><sup>+</sup> EDI</td><td>AMJUEL H.4 2.2.11</td></tr><tr><td>H + H → H + H</td><td>Elastic Scattering (HES)</td><td>Atom-Atom HES</td><td>AMMONX R-H-H</td></tr><tr><td>H + H<sub>2</sub> → H + H<sub>2</sub></td><td>Elastic Scattering (HES)</td><td>Atom-Molecule HES</td><td>AMMONX R-H-H2</td></tr><tr><td>H<sub>2</sub> + H<sub>2</sub> → H<sub>2</sub> + H<sub>2</sub></td><td>Elastic Scattering (HES)</td><td>Molecule-Molecule HES</td><td>AMMONX R-H2-H2</td></tr></tbody></table>

<footnote><sup>a</sup> The standard naming system maintained by the IAEA is available at https://amdis.iaea.org/databases/processes/.</footnote>

<footnote><sup>b</sup> Notation in EIRENE databases (see www.eirene.de).</footnote>

<footnote><sup>c</sup> The runs are not resolved in terms of excitation levels, however this rate represents the reaction involving excited atom product.</footnote>
ORCID iDs

N. Rivals https://orcid.org/0000-0002-4467-5857

Y. Marandet https://orcid.org/0000-0002-9970-538X

J.-S. Park https://orcid.org/0000-0003-0871-7527

H. Bufferand https://orcid.org/0000-0002-9731-5855

R.A. Pitts https://orcid.org/0000-0001-9455-2698

H. Yang https://orcid.org/0000-0002-4983-1249

References

[1] Bonnin X., Dekeyser W., Pitts R., Coster D., Voskoboynikov S. and Wiesen S. 2016 Presentation of the new SOLPS-ITER code package for tokamak plasma edge modelling *Plasma Fusion Res.* **11** 1403102

[2] Pitts R.A. et al 2019 Physics basis for the first ITER tungsten divertor *Nucl. Mater. Energy* **20** 100696

[3] Park J.-S., Bonnin X. and Pitts R. 2021 Assessment of ITER divertor performance during early operation phases *Nucl. Fusion* **61** 016021

[4] I Organization 2018 ITER research plan within the staged approach (level III - provisional version) (available at: www.iter.org/doc/www/content/com/Lists/ITER%20Technical%20Reports/Attachments/9/ITER_Research_Plan_within_the_Staged_Approach levIII PROVersion.pdf)

[5] Bufferand H. et al 2021 Progress in edge plasma turbulence modelling-hierarchy of models from 2D transport application to 3D fluid simulations in realistic tokamak geometry *Nucl. Fusion* **61** 116052

[6] Rivals N. et al 2022 SOLEDGE3X full vessel plasma simulations for computation of ITER first-wall fluxes *Contrib. Plasma Phys.* **62** e202100182

[7] Rivals N. et al 2022 Impact of enhanced far-SOL transport on first wall fluxes in ITER from full vessel edge-plasma simulations *Nucl. Mater. Energy* **33** 101293

[8] Zhdanov V. 2002 *Transport Processes in Multicomponent Plasma* (CRC Press)

[9] Bufferand H. et al 2022 Implementation of multi-component Zhdanov closure in SOLEDGE3X *Plasma Phys. Control. Fusion* **64** 055001

[10] Reiter D., Baelmans M. and Börner P. 2005 The EIRENE and B2-EIRENE codes *Fusion Sci. Technol.* **47** 172–86

[11] Romazanov J. et al 2022 Beryllium erosion and redeposition in ITER H, He and D–T discharges *Nucl. Fusion* **62** 036011

[12] Stangeby P. et al 2003 Interpretive modeling of as-simple-as-possible plasma discharges on DIII-D using the OEDGE code *J. Nucl. Mater.* **313-316** 883–7

[13] Dekeyser W. et al 2021 Plasma edge simulations including realistic wall geometry with SOLPS-ITER *Nucl. Mater. Energy* **27** 100999

[14] Eich T. et al 2013 Scaling of the tokamak near the scrape-off layer H-mode power width and implications for ITER *Nucl. Fusion* **53** 093031

[15] Kotov V., Reiter D., Pitts R.A., Jachmich S., Huber A. and Coster D.P. JET-EFDA Contributors 2008 Numerical modelling of high density JET divertor plasma with the SOLPS4.2 (B2-EIRENE) code *Plasma Phys. Control. Fusion* **50** 105012

[16] Schmid K., Krieger K., Lisgo S., Meisl G. and Brezinsek S. JET Contributors 2015 WALLDYN simulations of global impurity migration in JET and extrapolations to ITER *Nucl. Fusion* **55** 053015

[17] Di Genova S. et al 2021 Modelling of tungsten contamination and screening in WEST plasma discharges *Nucl. Fusion* **61** 106019

[18] Yang H. et al 2022 Numerical modelling of the impact of leakage under divertor baffle in WEST *Nucl. Mater. Energy* **33** 101302

[19] Ciraolo G., Baschetti S., Bufferand H., Falchetto G., Ghendrih P., Yang H., Rivals N., Tamain P. and Serre E. 2022 A self-consistent cross-field transport model for edge tokamak plasma simulations: SOLEDGE modelling and comparison with experiments (available at: https://hal-cea.archives-ouvertes.fr/cea-03766468)

[20] Bufferand H. et al 2019 Description of the model application to 2D transport modelling application to 3D turbulence modelling (available at: https://hal.science/hal-02903014)

[21] Rivals N. 2023 Understanding density regimes and divertor detachment in ITER relevant conditions (available at: https://cea.hal.science/tel-04461920)

[22] Sharma P. and Hammett G.W. 2007 Preserving monotonicity in anisotropic diffusion *J. Comput. Phys.* **227** 123–42

[23] Dekeyser W., Bonnin X., Lisgo S.W., Pitts R.A. and LaBombard B. 2019 Implementation of a 9-point stencil in SOLPS-ITER and implications for Alcator C-Mod divertor plasma simulations *Nucl. Mater. Energy* **18** 125–30

[24] Kallenbach A. et al 2018 Parameter dependences of the separatrix density in nitrogen seeded ASDEX Upgrade H-mode discharges *Plasma Phys. Control. Fusion* **60** 045006

[25] Reksoatmodjo R., Mordijick S., Hughes J., Lore J. and Bonnin X. 2021 The role of edge fueling in determining the pedestal density in high neutral opacity Alcator C-Mod experiments *Nucl. Mater. Energy* **27** 100971

[26] Moscheni M. et al 2022 Cross-code comparison of the edge codes SOLPS-ITER, SOLEDGE2D and UEDGE in modelling a low-power scenario in the DTT *Nucl. Fusion* **62** 056009

[27] Stangeby P.C. 2000 *The Plasma Boundary of Magnetic Fusion Devices* (Plasma Physics Series) (Institute of Physics Publishing)

[28] Krasheninnikov S.I., Kukushkin A.S. and Pshenov A.A. 2016 Divertor plasma detachment *Phys. Plasmas* **23** 055602

[29] Kotov V. and Reiter D. 2009 Two-point analysis of the numerical modelling of detached divertor plasmas *Plasma Phys. Control. Fusion* **51** 115002

[30] Krasheninnikov S.I., Kukushkin A.S. and Pshenov A.A. 2022 Divertor plasma detachment: roles of plasma momentum, energy and particle balances *Plasma Phys. Control. Fusion* **64** 125011

[31] Kukushkin A.S. and Krasheninnikov S.I. 2019 Bifurcations and oscillations in divertor plasma *Plasma Phys. Control. Fusion* **61** 074001

[32] Kukushkin A., Pacher H. and Pitts R. 2015 Characteristics of divertor detachment for ITER conditions *J. Nucl. Mater.* **463** 586–90

[33] Kaveeva E. et al 2020 SOLPS-ITER modelling of ITER edge plasma with drifts and currents *Nucl. Fusion* **60** 046019

[34] Park J.-S., Bonnin X. and Rivals N. 2021 private communication

[35] Eich T. et al 2020 Turbulence driven widening of the near-SOL power width in ASDEX Upgrade H-Mode discharges *Nucl. Fusion* **60** 056016

[36] Krasheninnikov S.I. and Kukushkin A.S. 2022 1D model of tokamak scrape-off layer and divertor plasmas: Dimensionless parameters and the results of numerical simulations *Phys. Plasmas* **29** 122502

[37] Stangeby P. 2019 The roles of power loss and momentum-pressure loss in causing particle-detachment in tokamak divertors: I. A heuristic model analysis *Plasma Phys. Control. Fusion* **62** 025012
[38] Guillemaut C. et al 2014 Influence of atomic physics on EDGE2D-EIRENE simulations of JET divertor detachment with carbon and beryllium/tungsten plasma-facing components *Nucl. Fusion* **54** 093012

[39] Reiter D. 2023 The EIRENE code user manual (available at: http://eirene.de/Documentation/eirene.pdf)

[40] Eckstein W., Garcia-Rosales C., Roth J. and Ottenberger W. 1993 *Sputtering Data* (Max-Planck-Institut für Plasmaphysik)

[41] Brezinsek S. et al 2014 Study of physical and chemical assisted physical sputtering of beryllium in the JET ITER-like wall *Nucl. Fusion* **54** 103001

[42] Bodner G. et al 2022 Initial results from boron powder injection experiments in WEST lower single null L-mode plasmas *Nucl. Fusion* **62** 086020

[43] Krieger K., Balden M., Bortolon A., Dux R., Griener M., Hegele K., Laggner F., Rohde V. and Wampler W. 2023 Wall conditioning effects and boron migration during boron powder injection in ASDEX upgrade *Nucl. Mater. Energy* **34** 101374

[44] Hey J., Chu C. and Hintz E. 2000 Spectroscopic studies of cold atomic hydrogen and deuterium produced in a tokamak edge plasma *Contrib. Plasma Phys.* **40** 9–22

[45] Marandet Y., Genesio P., Koubiti M., Godbert-Mouret L., Felts B., Stamm R., Capes H. and Guirlet R. 2004 Characterization of tokamak edge plasmas using spectroscopic line profiles *Nucl. Fusion* **44** S118–22

[46] Sureshkumar S. et al 2024 First SOLEDGE3X-EIRENE simulations of the ITER Neon seeded burning plasma boundary up to the first wall *Nucl. Mater. Energy* **41** 101780

[47] Rowan W.L., Bengtson R., Bonnin X., Edmonds P., Hurwitz P., Solano E., Tsui H., Uglum J. and Wootton A. 1995 Particle balance in diverted plasmas in TEXT-U *J. Nucl. Mater.* **220-222** 668–71

[48] Vianello N. et al 2020 Scrape-off layer transport and filament characteristics in high-density tokamak regimes *Nucl. Fusion* **60** 016001

[49] Wynn A. et al 2018 Investigation into the formation of the scrape-off layer density shoulder in JET ITER-like wall l-mode and h-mode plasmas *Nucl. Fusion* **58** 056001