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

HAL
open science

# Deuterium uptake, desorption and sputtering from W(110) surface covered with oxygen

E.A. Hodille, B. Pavec, J. Denis, A. Dunand, Y. Ferro, Marco Minissale, Thierry Angot, Christian Grisolia, Régis Bisson

► To cite this version:

E.A. Hodille, B. Pavec, J. Denis, A. Dunand, Y. Ferro, et al. Deuterium uptake, desorption and sputtering from W(110) surface covered with oxygen. Nuclear Fusion, 2024, 64 (4), pp.046022. 10.1088/1741-4326/ad2a29. hal-04521576

HAL Id: hal-04521576

https://hal.science/hal-04521576

Submitted on 26 Mar 2024

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# Deuterium uptake, desorption and sputtering from W(110) surface covered with oxygen

To cite this article: E.A. Hodille et al 2024 *Nucl. Fusion* **64** 046022

View the [article online](#) for updates and enhancements.

You may also like

- Water activation products generation and transport in DEMO divertor

P. Chiovaro, A. Quartararo, P. Avona et al.

- Parallel flows as a key component to interpret Super-X divertor experiments

M. Carpita, O. Février, H. Reimerdes et al.

- Core ion measurements with collective Thomson scattering for DEMO burn control

J. Rasmussen and S.B. Korsholm
PAPER • OPEN ACCESS

# Deuterium uptake, desorption and sputtering from W(110) surface covered with oxygen

To cite this article: E.A. Hodille et al 2024 *Nucl. Fusion* **64** 046022

View the [article online](#) for updates and enhancements.

You may also like

- Water activation products generation and transport in DEMO divertor

P. Chiovaro, A. Quartararo, P. Avona et al.

- Parallel flows as a key component to interpret Super-X divertor experiments

M. Carpita, O. Février, H. Reimerdes et al.

- Core ion measurements with collective Thomson scattering for DEMO burn control

J. Rasmussen and S.B. Korsholm
# Deuterium uptake, desorption and sputtering from W(110) surface covered with oxygen

E.A. Hodille<fnref n="1" />,*, B. Pavec<fnref n="1" />,<fnref n="2" />, J. Denis<fnref n="3" />, A. Dunand<fnref n="3" />, Y. Ferro<fnref n="3" />, M. Minissale<fnref n="3" />, T. Angot<fnref n="3" />, C. Grisolia<fnref n="1" /> and R. Bisson<fnref n="3" />

E-mail: etienne.hodille@cea.fr

Received 11 December 2023, revised 30 January 2024

Accepted for publication 16 February 2024

Published 5 March 2024

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

CrossMark

## Abstract

Rate equation modelling is performed to simulate $D_2$ and $D_2 + D_2^+$ exposure of the W(110) surface with varying coverage of oxygen atoms (O) from the clean surface up to 0.75 monolayer of O. Density Functional Theory (DFT) calculated energetics are used as inputs for the surface processes and desorption energies are optimized to best reproduce the Thermal Desorption Spectrometry (TDS) experiments obtained for $D_2$ exposure. For the clean surface, the optimized desorption energies (1.10 eV–1.40 eV) are below the DFT ones (1.30 eV–1.50 eV). For the O covered surface, the main desorption peak is reproduced with desorption energies of 1.10 eV and 1.00 eV for 0.50 and 0.75 monolayer of O respectively. This is slightly higher than the DFT predicted desorption energies. In order to simulate satisfactorily the total retention obtained experimentally for $D_2 + D_2^+$ exposure, a sputtering process needs to be added to the model, describing the sputtering of adsorbed species (D atoms) by the incident D ions. The impact of the sputtering process on the shape of the TDS spectra, on the total retention and on the recycling of D from the wall is discussed. In order to better characterize the sputtering process, especially its products and yields, atomistic calculations such as molecular dynamics are suggested as a next step for this study.

Keywords: plasma-wall interactions, hydrogen, multi-scale modelling

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

## 1. Introduction

In a fusion reactor, Plasma Facing Components (PFCs) are exposed to high heat flux ($\approx 10 \text{ MWm}^{-2}$) carried in part by an intense fluxes of hydrogen (H) isotopes [1, 2]. The most

intense part of the plasma wall interactions occurs in the diver- tor for which the PFCs material is foreseen to be polycrystal- line tungsten (W) [3]. This intense particle flux leads to the implantation and trapping of hydrogen isotopes in PFCs caus- ing two main issues. First, the retention of tritium in PFCs (the so-called tritium inventory) creates a safety issue as tritium is radioactive. Then, the fuel retained in the wall can desorb thermally to the plasma (the so-called fuel recycling), which may act as an uncontrollable source of fuel affecting the con- trol of the core plasma density (and thus of the overall reactiv- ity of the plasma) [4]. Thus, controlling the inventory of hydro- gen isotopes in the PFCs and the release processes from the

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

Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ CEA, IRFM, F13108 Saint Paul Lez Durance, France</footnote>

<footnote n="2">² Ecole Polytechnique, Route de Saclay, 91128 Palaiseau Cedex, France</footnote>

<footnote n="3">³ Aix-Marseille Univ, CNRS, PIIM, Marseille, France</footnote>

<footnote>* Author to whom any correspondence should be addressed.</footnote>
# Deuterium uptake, desorption and sputtering from W(110) surface covered with oxygen

E.A. Hodille<fnref n="1" />,*, B. Pavec<fnref n="1" />,<fnref n="2" />, J. Denis<fnref n="3" />, A. Dunand<fnref n="3" />, Y. Ferro<fnref n="3" />, M. Minissale<fnref n="3" />, T. Angot<fnref n="3" />, C. Grisolia<fnref n="1" /> and R. Bisson<fnref n="3" />

E-mail: etienne.hodille@cea.fr

Received 11 December 2023, revised 30 January 2024

Accepted for publication 16 February 2024

Published 5 March 2024

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

CrossMark

## Abstract

Rate equation modelling is performed to simulate $D_2$ and $D_2 + D_2^+$ exposure of the W(110) surface with varying coverage of oxygen atoms (O) from the clean surface up to 0.75 monolayer of O. Density Functional Theory (DFT) calculated energetics are used as inputs for the surface processes and desorption energies are optimized to best reproduce the Thermal Desorption Spectrometry (TDS) experiments obtained for $D_2$ exposure. For the clean surface, the optimized desorption energies (1.10 eV–1.40 eV) are below the DFT ones (1.30 eV–1.50 eV). For the O covered surface, the main desorption peak is reproduced with desorption energies of 1.10 eV and 1.00 eV for 0.50 and 0.75 monolayer of O respectively. This is slightly higher than the DFT predicted desorption energies. In order to simulate satisfactorily the total retention obtained experimentally for $D_2 + D_2^+$ exposure, a sputtering process needs to be added to the model, describing the sputtering of adsorbed species (D atoms) by the incident D ions. The impact of the sputtering process on the shape of the TDS spectra, on the total retention and on the recycling of D from the wall is discussed. In order to better characterize the sputtering process, especially its products and yields, atomistic calculations such as molecular dynamics are suggested as a next step for this study.

Keywords: plasma-wall interactions, hydrogen, multi-scale modelling

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

## 1. Introduction

In a fusion reactor, Plasma Facing Components (PFCs) are exposed to high heat flux ($\approx 10 \text{ MWm}^{-2}$) carried in part by an intense fluxes of hydrogen (H) isotopes [1, 2]. The most

intense part of the plasma wall interactions occurs in the diver- tor for which the PFCs material is foreseen to be polycrystal- line tungsten (W) [3]. This intense particle flux leads to the implantation and trapping of hydrogen isotopes in PFCs caus- ing two main issues. First, the retention of tritium in PFCs (the so-called tritium inventory) creates a safety issue as tritium is radioactive. Then, the fuel retained in the wall can desorb thermally to the plasma (the so-called fuel recycling), which may act as an uncontrollable source of fuel affecting the con- trol of the core plasma density (and thus of the overall reactiv- ity of the plasma) [4]. Thus, controlling the inventory of hydro- gen isotopes in the PFCs and the release processes from the

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

Original Content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ CEA, IRFM, F13108 Saint Paul Lez Durance, France</footnote>

<footnote n="2">² Ecole Polytechnique, Route de Saclay, 91128 Palaiseau Cedex, France</footnote>

<footnote n="3">³ Aix-Marseille Univ, CNRS, PIIM, Marseille, France</footnote>

<footnote>* Author to whom any correspondence should be addressed.</footnote>
material are crucial for the safe operation of a fusion reactor. For this purpose, macroscopic models have been developed to describe the kinetics of H retention and outgassing from PFCs [2, 5, 6] under tokamak conditions. These models are usually parametrized by simulating laboratory experiments like Thermal Desorption Spectrometry (TDS) [7–9] and can take inputs from atomistic scale calculations like Density Functional Theory (DFT) or Molecular Dynamics (MD) for the energetics of the H/material interactions.

The interaction of H with the W surface can potentially affect the release of fuel from the W PFCs. It is thus important to accurately model this interaction of H with the surface to have a good estimation of the overall fuel retention. In order to describe more accurately the different processes at the surface, recent rate equation models [10–15] started to use the kinetic surface model based on Pick and Sonnenberg model [16] which describes the different transition rates from the sub-surface to the surface and the recombination process. This model can straightforwardly introduce the energetics of the surface processes calculated by DFT [17–20] which allows a comparison of the DFT energetics with experimental TDS data.

In this work, we simulate the experiment presented by Dunand et al [21] where single crystal tungsten oriented along the (110) surface has been exposed to both D₂ and D₂⁺ + D₂ at room temperature. The W(110) surface was prepared with different coverages of adsorbed oxygen (O), from clean surface up to 0.75 monolayer (ML) of O. The D₂ exposure suggests that the presence of adsorbed O blocks the sticking of D and reduces the desorption energies, as already captured by DFT calculations [19], ab-initio MD calculations [22] and previous experiments [23]. The D₂⁺ + D₂ exposure shows a reduction of the total retention which suggests the existence of another process which removes particles during D₂⁺ interaction with the surface. The objectives of this paper is to simulate these experimental results to (i) compare the DFT results from [17, 19] with the experimental TDS, (ii) determine the trapping parameters in well prepared W(110) single crystal and (iii) add a new surface process (involving sputtering of adsorbed D by D₂⁺) to reproduce the difference of retention between D₂ and D₂⁺ + D₂ exposure. To do so, the MHIMS (Migration of Hydrogen Isotopes in MaterialS) code [7] is used with the already implemented kinetic surface model [11, 24] which has been modified to allow D₂ exposure [25] and to which the sputtering process is added.

## 2. Model description

The MHIMS code [7, 11] intends to simulate the transport of any H isotopes in materials. The model implemented in the code has already been presented and discussed before [7, 11, 24–27]. In this study, different features of the code are improved, notably with the addition of different terms.

The model implemented in MHIMS can be divided in two parts:

(i) The transport of D in the bulk of the material, which couples the diffusion of D and the interaction with bulk defects.

(ii) The interaction of D with the surface of the material. This part acts as the boundary condition on the boundary $\partial\Omega$.

In this work, we simulate a single crystal of W with W(110) surface exposed to D₂⁺ ions and D₂ molecules.

D₂⁺ ions are implanted in the bulk through a volume source S (m⁻³ s⁻¹) characterized by an incident flux $\phi_{inc}$ (m⁻²s⁻¹), an incident energy $E_{inc}$ (eV), and an incident angle $\alpha_{inc}$:

$$S = (1 - r(E_{inc}, \alpha_{inc})) \phi_{inc} f(E_{inc}, \alpha_{inc}) \quad (1)$$

where $r(E_{inc}, \alpha_{inc})$ (dimensionless) is the reflection of ions at the surface and $f(E_{inc}, \alpha_{inc})$ (m⁻¹) the spatial distribution of the implanted D in the material. The implanted flux $\phi_{imp} = (1 - r(E_{inc}, \alpha_{inc})) \phi_{inc}$ is the part of the incident flux which stays in the material, diffuses and can be trapped. Both $r(E_{inc}, \alpha_{inc})$ and $f(E_{inc}, \alpha_{inc})$ are evaluated with SRIM [28] while $\phi_{inc}$ is defined by the experimental parameters of [21].

D₂ molecules interact with the surface as a flux of molecules $\Gamma_{D₂}$ (D₂m⁻²s⁻¹) calculated from the D₂ pressure $p_{D₂}$ (Pa) and the temperature of the gas $T_{gas}$ (K):

$$\Gamma_{D_2} = \frac{p_{D_2}}{\sqrt{2\pi m_{D_2} k_B T_{gas}}} \quad (2)$$

with $m_{D₂}$ (kg) the mass of the molecule impinging the wall and $k_B = 8.6175 \times 10^{-5}$ eV K⁻¹ the Boltzmann constant.

The model detailed below is implemented in the MHIMS code in which the equations are solved numerically using the DLSODE package [29, 30].

### 2.1. Bulk model

The model describing H isotopes migration in the materials is based on the McNabb and Foster equations [31]. It considers two types of D atoms: the D atoms in interstitial sites (IS), which are mobile with the concentration $c_m$ (m⁻³) and the trapped D in the trapping site of type i; with the concentration $c_{t,i}$ (m⁻³). The equations governing the spatial and time evolutions of both quantities are given here:

$$\frac{\partial c_m}{\partial t} = \nabla \cdot (D(T) \nabla c_m) + S - \sum_i \frac{\partial c_{t,i}}{\partial t} \quad (3)$$

$$\frac{\partial c_{t,i}}{\partial t} = \nu_{t,i}(T) c_m (n_i - c_{t,i}) - \nu_{dt,i}(T) c_{t,i} \quad (4)$$

with the parameters described in table 1. To account for the difference of mass between H and D in the diffusivity, the diffusivity calculated for H by Fernandez et al [32] is scaled by $1/\sqrt{2}$ [34]. Since DFT calculates that H sits in tetrahedral position in bcc W [32, 33], $n_{IS} = 6\rho_W$. The trapping and detrapping pre-exponential frequencies are taken equal to $10^{13}$ s⁻¹ based on phonon calculations and harmonic transition state theory [32, 33].
material are crucial for the safe operation of a fusion reactor. For this purpose, macroscopic models have been developed to describe the kinetics of H retention and outgassing from PFCs [2, 5, 6] under tokamak conditions. These models are usually parametrized by simulating laboratory experiments like Thermal Desorption Spectrometry (TDS) [7–9] and can take inputs from atomistic scale calculations like Density Functional Theory (DFT) or Molecular Dynamics (MD) for the energetics of the H/material interactions.

The interaction of H with the W surface can potentially affect the release of fuel from the W PFCs. It is thus important to accurately model this interaction of H with the surface to have a good estimation of the overall fuel retention. In order to describe more accurately the different processes at the surface, recent rate equation models [10–15] started to use the kinetic surface model based on Pick and Sonnenberg model [16] which describes the different transition rates from the sub-surface to the surface and the recombination process. This model can straightforwardly introduce the energetics of the surface processes calculated by DFT [17–20] which allows a comparison of the DFT energetics with experimental TDS data.

In this work, we simulate the experiment presented by Dunand et al [21] where single crystal tungsten oriented along the (110) surface has been exposed to both D₂ and D₂⁺ + D₂ at room temperature. The W(110) surface was prepared with different coverages of adsorbed oxygen (O), from clean surface up to 0.75 monolayer (ML) of O. The D₂ exposure suggests that the presence of adsorbed O blocks the sticking of D and reduces the desorption energies, as already captured by DFT calculations [19], ab-initio MD calculations [22] and previous experiments [23]. The D₂⁺ + D₂ exposure shows a reduction of the total retention which suggests the existence of another process which removes particles during D₂⁺ interaction with the surface. The objectives of this paper is to simulate these experimental results to (i) compare the DFT results from [17, 19] with the experimental TDS, (ii) determine the trapping parameters in well prepared W(110) single crystal and (iii) add a new surface process (involving sputtering of adsorbed D by D₂⁺) to reproduce the difference of retention between D₂ and D₂⁺ + D₂ exposure. To do so, the MHIMS (Migration of Hydrogen Isotopes in MaterialS) code [7] is used with the already implemented kinetic surface model [11, 24] which has been modified to allow D₂ exposure [25] and to which the sputtering process is added.

## 2. Model description

The MHIMS code [7, 11] intends to simulate the transport of any H isotopes in materials. The model implemented in the code has already been presented and discussed before [7, 11, 24–27]. In this study, different features of the code are improved, notably with the addition of different terms.

The model implemented in MHIMS can be divided in two parts:

(i) The transport of D in the bulk of the material, which couples the diffusion of D and the interaction with bulk defects.

(ii) The interaction of D with the surface of the material. This part acts as the boundary condition on the boundary $\partial\Omega$.

In this work, we simulate a single crystal of W with W(110) surface exposed to D₂⁺ ions and D₂ molecules.

D₂⁺ ions are implanted in the bulk through a volume source S (m⁻³ s⁻¹) characterized by an incident flux $\phi_{inc}$ (m⁻²s⁻¹), an incident energy $E_{inc}$ (eV), and an incident angle $\alpha_{inc}$:

$$S = (1 - r(E_{inc}, \alpha_{inc})) \phi_{inc} f(E_{inc}, \alpha_{inc}) \quad (1)$$

where $r(E_{inc}, \alpha_{inc})$ (dimensionless) is the reflection of ions at the surface and $f(E_{inc}, \alpha_{inc})$ (m⁻¹) the spatial distribution of the implanted D in the material. The implanted flux $\phi_{imp} = (1 - r(E_{inc}, \alpha_{inc})) \phi_{inc}$ is the part of the incident flux which stays in the material, diffuses and can be trapped. Both $r(E_{inc}, \alpha_{inc})$ and $f(E_{inc}, \alpha_{inc})$ are evaluated with SRIM [28] while $\phi_{inc}$ is defined by the experimental parameters of [21].

D₂ molecules interact with the surface as a flux of molecules $\Gamma_{D₂}$ (D₂m⁻²s⁻¹) calculated from the D₂ pressure $p_{D₂}$ (Pa) and the temperature of the gas $T_{gas}$ (K):

$$\Gamma_{D_2} = \frac{p_{D_2}}{\sqrt{2\pi m_{D_2} k_B T_{gas}}} \quad (2)$$

with $m_{D₂}$ (kg) the mass of the molecule impinging the wall and $k_B = 8.6175 \times 10^{-5}$ eV K⁻¹ the Boltzmann constant.

The model detailed below is implemented in the MHIMS code in which the equations are solved numerically using the DLSODE package [29, 30].

### 2.1. Bulk model

The model describing H isotopes migration in the materials is based on the McNabb and Foster equations [31]. It considers two types of D atoms: the D atoms in interstitial sites (IS), which are mobile with the concentration $c_m$ (m⁻³) and the trapped D in the trapping site of type i; with the concentration $c_{t,i}$ (m⁻³). The equations governing the spatial and time evolutions of both quantities are given here:

$$\frac{\partial c_m}{\partial t} = \nabla \cdot (D(T) \nabla c_m) + S - \sum_i \frac{\partial c_{t,i}}{\partial t} \quad (3)$$

$$\frac{\partial c_{t,i}}{\partial t} = \nu_{t,i}(T) c_m (n_i - c_{t,i}) - \nu_{dt,i}(T) c_{t,i} \quad (4)$$

with the parameters described in table 1. To account for the difference of mass between H and D in the diffusivity, the diffusivity calculated for H by Fernandez et al [32] is scaled by $1/\sqrt{2}$ [34]. Since DFT calculates that H sits in tetrahedral position in bcc W [32, 33], $n_{IS} = 6\rho_W$. The trapping and detrapping pre-exponential frequencies are taken equal to $10^{13}$ s⁻¹ based on phonon calculations and harmonic transition state theory [32, 33].
<table><caption>**Table 1.** Parameters of the bulk model. The W atomic concentration is calculated considering a lattice constant of 3.16 Å.</caption><tr><td><i>D</i>(<i>T</i>)</td><td>Interstitial D diffusivity (m<sup>2</sup>s<sup>-1</sup>)</td><td><sup>1.9×10<sup>-7</sup></sup>&frasl;<sub>√2</sub> exp(-<sup>0.2</sup>&frasl;<sub>k<sub>B</sub>T</sub>)</td><td>[32]</td></tr><tr><td><i>ν</i><sub>t,i</sub></td><td>Trapping rate constant (m<sup>3</sup>s<sup>-1</sup>)</td><td><i>ν</i><sub>0</sub>/<i>n</i><sub>IS</sub> exp(-<sup><i>E</i><sub>t,i</sub></sup>&frasl;<sub>k<sub>B</sub>T</sub>)</td><td>[25]</td></tr><tr><td><i>ν</i><sub>dt,i</sub></td><td>Detrapping rate constant (s<sup>-1</sup>)</td><td><i>ν</i><sub>0</sub> exp(-<sup><i>E</i><sub>dt,i</sub></sup>&frasl;<sub>k<sub>B</sub>T</sub>)</td><td>[7]</td></tr><tr><td><i>n</i><sub>i</sub></td><td>Trap concentration (m<sup>-3</sup>)</td><td>To be determined</td><td></td></tr><tr><td><i>ν</i><sub>0</sub></td><td>Pre-exponential (de)trapping frequency (s<sup>-1</sup>)</td><td>10<sup>13</sup> s<sup>-1</sup></td><td>[32, 33]</td></tr><tr><td><i>E</i><sub>dt,i</sub></td><td>Detrapping energy from trap <i>i</i> (eV)</td><td>To be determined</td><td></td></tr><tr><td><i>E</i><sub>t,i</sub></td><td>Trapping energy to trap <i>i</i> (eV)</td><td>0.2 eV</td><td>[11]</td></tr><tr><td><i>n</i><sub>IS</sub></td><td>Concentration of interstitial site (m<sup>-3</sup>)</td><td>6<i>ρ</i><sub>W</sub></td><td>[32, 33]</td></tr><tr><td><i>ρ</i><sub>W</sub></td><td>W atomic concentration (m<sup>-3</sup>)</td><td>6.3382 × 10<sup>28</sup> m<sup>-3</sup></td><td></td></tr></table>

<table><caption>**Table 2.** Expression of the different flux $\phi_P$ in equations (5) and (6) which are reported in [11, 24, 25].</caption><thead><tr><th colspan="3">φ<sub>P</sub> (m<sup>-2</sup>s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>φ<sub>gas</sub></td><td>2Γ<sub>D<sub>2</sub></sub>s(T)(1 − θ<sub>D</sub>)<sup>2</sup></td><td>D atoms sticking from the gas phase</td></tr><tr><td>φ<sub>des</sub></td><td>2ν<sub>des</sub>(T)c<sub>surf</sub><sup>2</sup></td><td>Desorbing D<sub>2</sub> to the gas phase</td></tr><tr><td>φ<sub>s→b</sub></td><td>ν<sub>s→b</sub>(T)c<sub>surf</sub></td><td>Absorption of adsorbed D to the bulk</td></tr><tr><td>φ<sub>b→s</sub></td><td>ν<sub>b→s</sub>(T)c<sub>m</sub>(∂Ω)(1 − θ<sub>D</sub>)</td><td>(Re-)surfacing: Bulk to surface transition</td></tr></tbody></table>

## 2.2. Surface model

### 2.2.1. Description.
The interaction of D with the W surface is based on the description given by Pick and Sonnenberg [16] and first introduced in MHIMS to simulate atomic D exposure in W [11]. Similar models can be found in other rate equation codes analogous to MHIMS in the literature such as TESSIM(-X) [12, 14] or HIIPC [35]. Since then, the model implemented in MHIMS has been upgraded to take into account the exposure to D<sub>2</sub> gas [25]. It considers two types of D atoms: the adsorbed D atoms with the concentration $c_{surf}$ (m<sup>-2</sup>) and the interstitial D atoms at the boundary $\partial\Omega$, i.e. just beneath the surface, with the concentration $c_m(\partial\Omega)$ (m<sup>-3</sup>). The equations governing the spatial and time evolutions of both quantities are given here:

$$ \frac{dc_{\text{surf}}}{dt} = \varphi_{\text{gas}} - \varphi_{\text{des}} - \varphi_{\text{sput}} - \varphi_{\text{s→b}} + \varphi_{\text{b→s}} \quad (5) $$

$$ \lambda \frac{dc_m(\partial\Omega)}{dt} = \varphi_{s\to b} - \varphi_{b\to s} + D(T)(\nabla c_m)_{\partial\Omega} \quad (6) $$

with $\lambda = 1.117$ Å the distance between two IS in W. The expression of the various fluxes in the right-hand-side of equations (5) and (6) can be, for the most part, found in [11, 24, 25] and are recalled in table 2. The last term of the right-hand-side of equation (6) represents the diffusive flux of interstitial D away from the surface. It connects the surface model to the bulk model described by equations (3) and (4).

The parameters involved in the expression of these fluxes are reported in table 3. The pre-exponential frequencies for desorption, bulk to surface and surface to bulk processes are taken equal to $10^{13}$ s<sup>-1</sup> based on phonon calculations and harmonic transition state theory [36].

In this work, we are investigating D desorption from W(110) with pre-adsorbed O atoms at 3 different coverages $\theta_O$: clean ($\theta_O = 0$), 0.50 ML of O ($\theta_O = 0.5$) and 0.75 ML

of O ($\theta_O = 0.75$). According to DFT calculations and thermodynamics models [18, 37, 38], the saturation occurs for 1 H/W on the W(110) surface. Thus, $n_{\text{surf}(\text{clean})} = 1.416 \times 10^{19}$ m<sup>-2</sup> considering a lattice constant for tungsten of 3.16 Å. According to experimental results [39] and recent DFT calculations [19], when O is pre-adsorbed on the surface, the H saturation occurs for $\theta_D + \theta_O = 1$. Thus, $n_{\text{surf}(\theta_O)} = n_{\text{surf}(\text{clean})}(1 - \theta_O)$ and the values for the three surfaces are reported in table 3.

### 2.2.2. Sputtering in MHIMS.
In this work, based on the suggestion made by Dunand et al [21], a sputtering process has been introduced in the model via the flux $\varphi_{\text{sput}}$ in equation (5). It represents the sputtering of adsorbed D (D<sub>ads</sub>) by incident D ions (D<sub>inc</sub>) from the plasma which is schematically described in figure 1.

Assuming the sputtering reaction D<sub>inc</sub> + D<sub>ads</sub> is an elementary process, its kinetic is proportional to the concentration of its reactants: $c_{\text{surf}}$ is the concentration of D<sub>ads</sub> and the incident flux $\phi_{\text{inc}}$ sets the amount of D<sub>inc</sub> reacting with the surface. For a fully covered surface, the frequency of the sputtering events is $\sigma_{\text{sput}}\phi_{\text{inc}}$ with $\sigma_{\text{sput}}$ (m<sup>2</sup>) the sputtering cross section. It yields:

$$ \varphi_{\text{sput}} = \sigma_{\text{sput}} \phi_{\text{inc}} c_{\text{surf}}. \quad (7) $$

This has the same expression as the abstraction model implemented in MHIMS [11]. Indeed, part of the ensemble of sputtering processes can be compared to the hot-atom and the Eley-Rideal recombinations [41] simulated by the abstraction model, in which the products of the reactions are D<sub>2</sub> molecules. However, the sputtering process can also release single D atoms going to the plasma and dedicated MD calculations are necessary to distinguish the products of the sputtering. The distinction is important as the kinetic energy and the
<table><caption>**Table 1.** Parameters of the bulk model. The W atomic concentration is calculated considering a lattice constant of 3.16 Å.</caption><tr><td>D(T)</td><td>Interstitial D diffusivity (m<sup>2</sup>s<sup>-1</sup>)</td><td><sup>1.9×10<sup>-7</sup></sup>&frasl;<sub>√2</sub> exp(-<sup>0.2</sup>&frasl;<sub>k<sub>B</sub>T</sub>)</td><td>[32]</td></tr><tr><td>ν<sub>t,i</sub></td><td>Trapping rate constant (m<sup>3</sup>s<sup>-1</sup>)</td><td>ν<sub>0</sub>/n<sub>IS</sub> exp(-<sup>E<sub>t,i</sub></sup>&frasl;<sub>k<sub>B</sub>T</sub>)</td><td>[25]</td></tr><tr><td>ν<sub>dt,i</sub></td><td>Detrapping rate constant (s<sup>-1</sup>)</td><td>ν<sub>0</sub> exp(-<sup>E<sub>dt,i</sub></sup>&frasl;<sub>k<sub>B</sub>T</sub>)</td><td>[7]</td></tr><tr><td>n<sub>i</sub></td><td>Trap concentration (m<sup>-3</sup>)</td><td>To be determined</td><td></td></tr><tr><td>ν<sub>0</sub></td><td>Pre-exponential (de)trapping frequency (s<sup>-1</sup>)</td><td>10<sup>13</sup> s<sup>-1</sup></td><td>[32, 33]</td></tr><tr><td>E<sub>dt,i</sub></td><td>Detrapping energy from trap i (eV)</td><td>To be determined</td><td></td></tr><tr><td>E<sub>t,i</sub></td><td>Trapping energy to trap i (eV)</td><td>0.2 eV</td><td>[11]</td></tr><tr><td>n<sub>IS</sub></td><td>Concentration of interstitial site (m<sup>-3</sup>)</td><td>6ρ<sub>W</sub></td><td>[32, 33]</td></tr><tr><td>ρ<sub>W</sub></td><td>W atomic concentration (m<sup>-3</sup>)</td><td>6.3382 × 10<sup>28</sup> m<sup>-3</sup></td><td></td></tr></table>

<table><caption>**Table 2.** Expression of the different flux $\phi_P$ in equations (5) and (6) which are reported in [11, 24, 25].</caption><thead><tr><th colspan="3">φ<sub>P</sub> (m<sup>-2</sup>s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>φ<sub>gas</sub></td><td>2Γ<sub>D<sub>2</sub></sub>s(T)(1 − θ<sub>D</sub>)<sup>2</sup></td><td>D atoms sticking from the gas phase</td></tr><tr><td>φ<sub>des</sub></td><td>2ν<sub>des</sub>(T)c<sub>surf</sub><sup>2</sup></td><td>Desorbing D<sub>2</sub> to the gas phase</td></tr><tr><td>φ<sub>s→b</sub></td><td>ν<sub>s→b</sub>(T)c<sub>surf</sub></td><td>Absorption of adsorbed D to the bulk</td></tr><tr><td>φ<sub>b→s</sub></td><td>ν<sub>b→s</sub>(T)c<sub>m</sub>(∂Ω)(1 − θ<sub>D</sub>)</td><td>(Re-)surfacing: Bulk to surface transition</td></tr></tbody></table>

## 2.2. Surface model

### 2.2.1. Description.
The interaction of D with the W surface is based on the description given by Pick and Sonnenberg [16] and first introduced in MHIMS to simulate atomic D exposure in W [11]. Similar models can be found in other rate equation codes analogous to MHIMS in the literature such as TESSIM(-X) [12, 14] or HIIPC [35]. Since then, the model implemented in MHIMS has been upgraded to take into account the exposure to D<sub>2</sub> gas [25]. It considers two types of D atoms: the adsorbed D atoms with the concentration $c_{surf}$ (m<sup>-2</sup>) and the interstitial D atoms at the boundary $\partial\Omega$, i.e. just beneath the surface, with the concentration $c_m(\partial\Omega)$ (m<sup>-3</sup>). The equations governing the spatial and time evolutions of both quantities are given here:

$$ \frac{dc_{\text{surf}}}{dt} = \varphi_{\text{gas}} - \varphi_{\text{des}} - \varphi_{\text{sput}} - \varphi_{\text{s→b}} + \varphi_{\text{b→s}} \quad (5) $$

$$ \lambda \frac{dc_{\text{m}}(\partial\Omega)}{dt} = \varphi_{\text{s→b}} - \varphi_{\text{b→s}} + D(T)(\nabla c_{\text{m}})_{\partial\Omega} \quad (6) $$

with $\lambda = 1.117$ Å the distance between two IS in W. The expression of the various fluxes in the right-hand-side of equations (5) and (6) can be, for the most part, found in [11, 24, 25] and are recalled in table 2. The last term of the right-hand- side of equation (6) represents the diffusive flux of interstitial D away from the surface. It connects the surface model to the bulk model described by equations (3) and (4).

The parameters involved in the expression of these fluxes are reported in table 3. The pre-exponential frequencies for desorption, bulk to surface and surface to bulk processes are taken equal to $10^{13}$ s<sup>-1</sup> based on phonon calculations and harmonic transition state theory [36].

In this work, we are investigating D desorption from W(110) with pre-adsorbed O atoms at 3 different coverages $\theta_O$: clean ($\theta_O = 0$), 0.50 ML of O ($\theta_O = 0.5$) and 0.75 ML

of O ($\theta_O = 0.75$). According to DFT calculations and thermodynamics models [18, 37, 38], the saturation occurs for 1 H/W on the W(110) surface. Thus, $n_{\text{surf}(\text{clean})} = 1.416 \times 10^{19}$ m<sup>-2</sup> considering a lattice constant for tungsten of 3.16 Å. According to experimental results [39] and recent DFT calculations [19], when O is pre-adsorbed on the surface, the H saturation occurs for $\theta_D + \theta_O = 1$. Thus, $n_{\text{surf}(\theta_O)} = n_{\text{surf}(\text{clean})}(1 - \theta_O)$ and the values for the three surfaces are reported in table 3.

### 2.2.2. Sputtering in MHIMS.
In this work, based on the suggestion made by Dunand et al [21], a sputtering process has been introduced in the model via the flux $\phi_{\text{sput}}$ in equation (5). It represents the sputtering of adsorbed D (D<sub>ads</sub>) by incident D ions (D<sub>inc</sub>) from the plasma which is schematically described in figure 1.

Assuming the sputtering reaction D<sub>inc</sub> + D<sub>ads</sub> is an elementary process, its kinetic is proportional to the concentration of its reactants: $c_{\text{surf}}$ is the concentration of D<sub>ads</sub> and the incident flux $\phi_{\text{inc}}$ sets the amount of D<sub>inc</sub> reacting with the surface. For a fully covered surface, the frequency of the sputtering events is $\sigma_{\text{sput}}\phi_{\text{inc}}$ with $\sigma_{\text{sput}}$ (m<sup>2</sup>) the sputtering cross section. It yields:

$$ \varphi_{\text{sput}} = \sigma_{\text{sput}} \phi_{\text{inc}} c_{\text{surf}}. \quad (7) $$

This has the same expression as the abstraction model implemented in MHIMS [11]. Indeed, part of the ensemble of sputtering processes can be compared to the hot-atom and the Eley-Rideal recombinations [41] simulated by the abstraction model, in which the products of the reactions are D<sub>2</sub> molecules. However, the sputtering process can also release single D atoms going to the plasma and dedicated MD calculations are necessary to distinguish the products of the sputtering. The distinction is important as the kinetic energy and the
<table><caption>**Table 3.** Parameters to express the D fluxes at the surface.</caption><thead><tr><th>θ<sub>D</sub></th><th>D coverage of the surface (dimensionless)</th><th>c<sub>surf</sub>/n<sub>surf</sub></th></tr></thead><tbody><tr><td>Γ<sub>D<sub>2</sub></sub></td><td>D<sub>2</sub> flux on the surface (m<sup>-2</sup>s<sup>-1</sup>)</td><td>Equation (2)</td></tr><tr><td>s(T)</td><td>Sticking coefficient of D<sub>2</sub> (dimensionless)</td><td>s<sub>0</sub> exp(-E<sub>diff</sub>/k<sub>B</sub>T)</td></tr><tr><td>ν<sub>des</sub>(T)</td><td>Desorption rate constant (m<sup>2</sup>s<sup>-1</sup>)</td><td>ν<sub>0</sub>λ<sub>des</sub><sup>2</sup> exp(-E<sub>des</sub>/k<sub>B</sub>T)</td></tr><tr><td>ν<sub>s→b</sub>(T)</td><td>Surface to bulk rate constant (s<sup>-1</sup>)</td><td>ν<sub>0</sub> exp(-E<sub>s→b</sub>/k<sub>B</sub>T)</td></tr><tr><td>ν<sub>b→s</sub>(T)</td><td>Bulk to surface rate constant (ms<sup>-1</sup>)</td><td>ν<sub>0</sub>λ<sub>abs</sub> exp(-E<sub>b→s</sub>/k<sub>B</sub>T)</td></tr><tr><td>s<sub>0</sub></td><td>Pre-exponential sticking coefficient (dimensionless)</td><td>1.0</td></tr><tr><td>E<sub>des</sub>(θ<sub>D</sub>)</td><td>Energy barrier for D<sub>2</sub> desorption (eV)</td><td>Figure 2</td></tr><tr><td>E<sub>b→s</sub></td><td>Energy barrier from bulk to surface (eV)</td><td>E<sub>diff</sub></td></tr><tr><td>E<sub>diss</sub></td><td>Energy barrier for D<sub>2</sub> dissociation (eV)</td><td>0.0 eV</td></tr><tr><td>E<sub>s→b</sub></td><td>Energy barrier from surface to bulk (eV)</td><td>(E<sub>des</sub> - E<sub>diss</sub>) / 2 + E<sub>b→s</sub> - Q<sub>sol</sub></td></tr><tr><td>Q<sub>sol</sub></td><td>Heat of solution of H in IS in W (eV)</td><td>1.0 eV</td></tr><tr><td>λ<sub>des</sub></td><td>Average distance between adsorption sites (m)</td><td>1/√n<sub>surf</sub></td></tr><tr><td>λ<sub>abs</sub></td><td>Distance between subsurface IS and adsorption sites (m)</td><td>n<sub>surf</sub>/n<sub>IS</sub></td></tr><tr><td>n<sub>surf</sub></td><td>Concentration of adsorption sites (×10<sup>19</sup> m<sup>-2</sup>)</td><td>1.416 for clean W(110)<br>0.708 for θ<sub>O</sub> = 0.50<br>0.354 for θ<sub>O</sub> = 0.75</td></tr></tbody></table>

Incident D⁺

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Schematic representation of the sputtering process of adsorbed D (green) by incident D ions (red) on metallic (grey) surface. Two possible products are presented: a single sputtered D atom with an implanted D atom and a sputtered D atom recombined with the incident ion (whether directly or after its resurfacing). Other possible products may be obtained such as the recombination an adsorbed atom with a sputtered atom hoping on the surface.</figcaption></figure>

behavior of D atoms and D₂ molecules in the edge plasma are different in a fusion reactor.

The sputtering process is also quantified by the sputtering yield $Y_{\text{sput}}$ (dimensionless) which takes the following form: $Y_{\text{sput}} = \sigma_{\text{sput}}n_{\text{surf}}$. In that case, the sputtering flux is given as:

$$ \varphi_{\text{sput}} = Y_{\text{sput}} \phi_{\text{inc}} \theta_{\text{D}}. \qquad (8) $$

Thus, $Y_{\text{sput}}$ is the amount of sputtered D adsorbed atoms per incident D ions in the case of a saturated surface ($\theta_{\text{D}} = 1$). In the next, $Y_{\text{sput}}$ is a free parameter to reproduce the experimental results of Dunand et al [21].

### 2.2.3. Coverage dependent $E_{\text{des}}$.
DFT calculations [17, 20] and experiments [10, 23, 42] suggest that the desorption energy of H from the W(110) surface depends on the H surface coverage. This dependency has already been implemented in rate equation models by using a continuous function $E_{\text{des}}(\theta_{\text{D}})$ [15, 24, 43–45]. In MHIMS, this function had initially the form of a Fermi–Dirac distribution [24]. However, recent DFT calculations [19] shown the possibility for $E_{\text{des}}$ to drop further when the local coverage increases above the saturation limits. This can be rendered by the addition of an exponential decrease as implemented in CRDS for Be [43]. Thus, the updated expression of $E_{\text{des}}(\theta_{\text{D}})$ is:

$$ E_{\text{des}}(\theta_{\text{D}}) = E_{\text{FD}}(\theta_{\text{D}})(1 - \alpha \exp(\beta(1 - \theta_{\text{D}}))) \quad (9) $$

$$ E_{\text{FD}}(\theta_{\text{D}}) = E_{0} + \frac{\Delta E}{1 + \exp\left(\frac{\theta_{\text{D}} - \theta_{\text{D0}}}{\delta\theta_{\text{D}}}\right)}. \qquad (10) $$

In equation (9), $E_{\text{FD}}(\theta_{\text{D}})$ is the part already shown in [15, 24] and the second part is the exponential decrease where $\alpha$ (dimensionless) indicates the amplitude of the drop of $E_{\text{des}}$ and $\beta$ (dimensionless) indicate how fast it drops.

This oversaturation of H revealed at the W(110) surface by DFT can be created by H diffusion on the surface or from the bulk to the surface. The desorption at oversaturation calculated by DFT, $E_{\text{oversat}}^{\text{DFT}}$, is in the range of 0.7–0.8 eV for clean surface [17] and 0.6 eV for W(110) with pre-adsorbed O [19]. These values are used to calculate $\alpha$ as:
<table><caption>**Table 3.** Parameters to express the D fluxes at the surface.</caption><thead><tr><th>θ<sub>D</sub></th><th>D coverage of the surface (dimensionless)</th><th>c<sub>surf</sub>/n<sub>surf</sub></th></tr></thead><tbody><tr><td>Γ<sub>D<sub>2</sub></sub></td><td>D<sub>2</sub> flux on the surface (m<sup>-2</sup>s<sup>-1</sup>)</td><td>Equation (2)</td></tr><tr><td>s(T)</td><td>Sticking coefficient of D<sub>2</sub> (dimensionless)</td><td>s<sub>0</sub> exp(-E<sub>diff</sub>/k<sub>B</sub>T)</td><td>[15]</td></tr><tr><td>ν<sub>des</sub>(T)</td><td>Desorption rate constant (m<sup>2</sup>s<sup>-1</sup>)</td><td>ν<sub>0</sub>λ<sub>des</sub><sup>2</sup> exp(-E<sub>des</sub>/k<sub>B</sub>T)</td><td>[11]</td></tr><tr><td>ν<sub>s→b</sub>(T)</td><td>Surface to bulk rate constant (s<sup>-1</sup>)</td><td>ν<sub>0</sub> exp(-E<sub>s→b</sub>/k<sub>B</sub>T)</td><td>[11]</td></tr><tr><td>ν<sub>b→s</sub>(T)</td><td>Bulk to surface rate constant (ms<sup>-1</sup>)</td><td>ν<sub>0</sub>λ<sub>abs</sub> exp(-E<sub>b→s</sub>/k<sub>B</sub>T)</td><td>[11]</td></tr><tr><td>s<sub>0</sub></td><td>Pre-exponential sticking coefficient (dimensionless)</td><td>1.0</td><td>[15]</td></tr><tr><td>E<sub>des</sub>(θ<sub>D</sub>)</td><td>Energy barrier for D<sub>2</sub> desorption (eV)</td><td>Figure 2</td><td></td></tr><tr><td>E<sub>b→s</sub></td><td>Energy barrier from bulk to surface (eV)</td><td>E<sub>diff</sub></td><td>[11]</td></tr><tr><td>E<sub>diss</sub></td><td>Energy barrier for D<sub>2</sub> dissociation (eV)</td><td>0.0 eV</td><td>[19]</td></tr><tr><td>E<sub>s→b</sub></td><td>Energy barrier from surface to bulk (eV)</td><td>(E<sub>des</sub> - E<sub>diss</sub>) / 2 + E<sub>b→s</sub> - Q<sub>sol</sub></td><td>[16]</td></tr><tr><td>Q<sub>sol</sub></td><td>Heat of solution of H in IS in W (eV)</td><td>1.0 eV</td><td>[32, 40]</td></tr><tr><td>λ<sub>des</sub></td><td>Average distance between adsorption sites (m)</td><td>1/√n<sub>surf</sub></td><td>[11]</td></tr><tr><td>λ<sub>abs</sub></td><td>Distance between subsurface IS and adsorption sites (m)</td><td>n<sub>surf</sub>/n<sub>IS</sub></td><td>[11]</td></tr><tr><td>n<sub>surf</sub></td><td>Concentration of adsorption sites (×10<sup>19</sup> m<sup>-2</sup>)</td><td>1.416 for clean W(110)<br>0.708 for θ<sub>O</sub> = 0.50<br>0.354 for θ<sub>O</sub> = 0.75</td><td></td></tr></tbody></table>

Incident D⁺

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Schematic representation of the sputtering process of adsorbed D (green) by incident D ions (red) on metallic (grey) surface. Two possible products are presented: a single sputtered D atom with an implanted D atom and a sputtered D atom recombined with the incident ion (whether directly or after its resurfacing). Other possible products may be obtained such as the recombination an adsorbed atom with a sputtered atom hoping on the surface.</figcaption></figure>

behavior of D atoms and D₂ molecules in the edge plasma are different in a fusion reactor.

The sputtering process is also quantified by the sputtering yield $Y_{\text{sput}}$ (dimensionless) which takes the following form: $Y_{\text{sput}} = \sigma_{\text{sput}}n_{\text{surf}}$. In that case, the sputtering flux is given as:

$$ \varphi_{\text{sput}} = Y_{\text{sput}} \phi_{\text{inc}} \theta_{\text{D}}. \qquad (8) $$

Thus, $Y_{\text{sput}}$ is the amount of sputtered D adsorbed atoms per incident D ions in the case of a saturated surface ($\theta_{\text{D}} = 1$). In the next, $Y_{\text{sput}}$ is a free parameter to reproduce the experimental results of Dunand et al [21].

### 2.2.3. Coverage dependent $E_{\text{des}}$.
DFT calculations [17, 20] and experiments [10, 23, 42] suggest that the desorption energy of H from the W(110) surface depends on the H surface coverage. This dependency has already been implemented in rate equation models by using a continuous function $E_{\text{des}}(\theta_{\text{D}})$ [15, 24, 43–45]. In MHIMS, this function had initially the form of a Fermi–Dirac distribution [24]. However, recent DFT calculations [19] shown the possibility for $E_{\text{des}}$ to drop further when the local coverage increases above the saturation limits. This can be rendered by the addition of an exponential decrease as implemented in CRDS for Be [43]. Thus, the updated expression of $E_{\text{des}}(\theta_{\text{D}})$ is:

$$ E_{\text{des}}(\theta_{\text{D}}) = E_{\text{FD}}(\theta_{\text{D}}) (1 - \alpha \exp(\beta (1 - \theta_{\text{D}}))) \quad (9) $$

$$ E_{\text{FD}}(\theta_{\text{D}}) = E_{0} + \frac{\Delta E}{1 + \exp\left(\frac{\theta_{\text{D}} - \theta_{\text{D0}}}{\delta \theta_{\text{D}}}\right)}. \qquad (10) $$

In equation (9), $E_{\text{FD}}(\theta_{\text{D}})$ is the part already shown in [15, 24] and the second part is the exponential decrease where $\alpha$ (dimensionless) indicates the amplitude of the drop of $E_{\text{des}}$ and $\beta$ (dimensionless) indicate how fast it drops.

This oversaturation of H revealed at the W(110) surface by DFT can be created by H diffusion on the surface or from the bulk to the surface. The desorption at oversaturation calculated by DFT, $E_{\text{oversat}}^{\text{DFT}}$, is in the range of 0.7–0.8 eV for clean surface [17] and 0.6 eV for W(110) with pre-adsorbed O [19]. These values are used to calculate $\alpha$ as:
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Evolution of $E_{des}$ with $c_{surf}$ for the three coverage of pre-adsorbed O. The symbols are DFT data from a [17] and b [19], the dashed line are $E_{des}^{DFT}$ which fit the DFT data with equation (9) and the solid lines are the optimized $E_{des}(\theta_D)$ evolution to reproduce the experimental data (see section 4). For the DFT data from [17], two data points are given per $c_{surf}$ as they calculated the desorption with two methods: nudge elastic band method or the difference of energy between the adsorbate state and the gas state.</figcaption></figure>

$$ \alpha = 1 - \frac{E_{\text{oversat}}^{\text{DFT}}}{E_{\text{FD}}(\theta_{\text{D}} = 1)}. \quad (11) $$

The remaining parameters of $E_{des}$ are obtained by fitting the DFT data from [17, 19]. The obtained evolutions of $E_{des}(\theta_D)$ are reported in figure 2 as a function of $c_{surf}$ for the three surfaces investigated. The effect of the pre-adsorbed O is seen by the maximum achievable D coverage which decreases. These parameters are then optimized to reproduce the experimental data which will be described further in the next sections.

# 3. Method

## 3.1. Description of the simulated scenario

In [21], a W(110) sample with three different coverages of pre-adsorbed O is exposed to $D_2$ molecules or $D_2$ molecules and $D_2^+$ ions. The parameters of the exposure to both species are reported in table 4. For the three different O coverages, the simulation scenario contains three steps:

(i) The exposure to external sources of D (either $D_2$ only or both $D_2$ and $D_2^+$) lasts $t_{imp} = 3000$ s at the surface temperature $T_{imp}$, the pressure is $p_{D_2}$ and the incidence flux (if $D_2^+$ are present) is $\phi_{inc}$ (see table 4).

(ii) The storage: the sample resides in vacuum for $t_{storage} = 1$ h at $T_{storage}$.

(iii) TDS is realized up to 800 K with a ramp of 5 K s⁻¹.

<table><caption>**Table 4.** First row: experimental parameters for the exposure to $D_2$ (only $p_{D_2}$), implantation ($p_{D_2}$ and $\phi_{inc}$) and storage (between the exposure and the TDS) in Dunand et al [21]. Second row: parameters calculated for the MHIMS simulations.</caption><tr><td></td><td></td></tr><tr><td>p<sub>D<sub>2</sub></sub></td><td>2 &times; 10<sup>-5</sup> Pa</td></tr><tr><td>&phi;<sub>inc</sub></td><td>0.5 &times; 10<sup>18</sup> D<sub>2</sub><sup>+</sup> m<sup>-2</sup>s<sup>-1</sup></td></tr><tr><td>E<sub>inc</sub></td><td>250 eV D<sup>-1</sup></td></tr><tr><td>&alpha;<sub>inc</sub></td><td>normal incidence</td></tr><tr><td>fluence</td><td>3 &times; 10<sup>21</sup> Dm<sup>-2</sup></td></tr><tr><td>T<sub>imp</sub></td><td>300 K</td></tr><tr><td>t<sub>storage</sub></td><td>1 h</td></tr><tr><td>T<sub>storage</sub></td><td>300 K</td></tr><tr><td></td><td></td></tr><tr><td>r(E<sub>inc</sub>, &alpha;<sub>inc</sub>)</td><td>0.539</td></tr><tr><td>&phi;<sub>imp</sub></td><td>4.61 &times; 10<sup>17</sup> Dm<sup>-2</sup>s<sup>-1</sup></td></tr><tr><td>&Gamma;<sub>D<sub>2</sub></sub></td><td>1.52 &times; 10<sup>18</sup> D<sub>2</sub>m<sup>-2</sup>s<sup>-1</sup></td></tr><tr><td>t<sub>imp</sub></td><td>3000 s</td></tr><tr><td></td><td></td></tr></table>

## 3.2. Modelling steps

To identify the relevant model parameters, simulation of the experimental data from Dunand et al [21] are performed. These data being presented as raw spectrometer signal in [21], a quantitative evaluation is made which is presented in appendix A. The parametrization is performed through an optimization routine given in appendix B and inspired from [9, 46].

The parameters that one need to optimize can be divided in three groups:

1. The surface parameters $P_{sim}^{surf} = \{E_0, \Delta E, \theta_{D0}, \delta\theta_D, \beta\}$ which set the evolution of the surface energies $E_{des}(\theta_D)$ with equation (9). One set of parameters has to be determined for each surface (clean, $O_{0.50ML}$ and $O_{0.75ML}$).

2. The trap parameters $P_{sim}^{trap,i} = \{E_{dt,i}, n_i\}$ and the number of traps. We assume that the oxygen coverage does not impact the bulk trap parameters, hence only one set of trap parameters has to be determined.

3. The sputtering yields of adsorbed deuterium $P_{sim}^{sput} = \{Y_{sput}\}$ for each surface (clean, $O_{0.50ML}$ and $O_{0.75ML}$).

First, the surface parameters are determined using the results of exposure to $D_2$ for which D atoms will not enter the bulk. Then, the trap parameters are determined using the results of exposure to $D_2$ and $D_2^+$ since this is the only way to insert D atoms in the bulk traps at 300 K. According to Dunand et al [21], the sputtering process seems to have less impact during the $D_2 + D_2^+$ exposure of the $O_{0.75ML}$ surface. Thus, the trap parameters will be determined for the $O_{0.75ML}$ surface exposed to $D_2$ and $D_2^+$. For this case, the simulations revealed that $Y_{sput}^{O_{0.75ML}}$ cannot be determined unambiguously and it will be varied from 0.000 to 3.000 and an optimized set of trap parameters will be obtained for each values of $Y_{sput}^{O_{0.75ML}}$. Finally, the sputtering yields for the clean surface and the $O_{0.50ML}$ surface will be determined for each set of trap parameters (hence for different values of $Y_{sput}^{O_{0.75ML}}$).
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Evolution of $E_{des}$ with $c_{surf}$ for the three coverage of pre-adsorbed O. The symbols are DFT data from a [17] and b [19], the dashed line are $E_{des}^{DFT}$ which fit the DFT data with equation (9) and the solid lines are the optimized $E_{des}(\theta_D)$ evolution to reproduce the experimental data (see section 4). For the DFT data from [17], two data points are given per $c_{surf}$ as they calculated the desorption with two methods: nudge elastic band method or the difference of energy between the adsorbate state and the gas state.</figcaption></figure>

$$ \alpha = 1 - \frac{E_{\text{oversat}}^{\text{DFT}}}{E_{\text{FD}}(\theta_{\text{D}} = 1)}. \quad (11) $$

The remaining parameters of $E_{des}$ are obtained by fitting the DFT data from [17, 19]. The obtained evolutions of $E_{des}(\theta_D)$ are reported in figure 2 as a function of $c_{surf}$ for the three surfaces investigated. The effect of the pre-adsorbed O is seen by the maximum achievable D coverage which decreases. These parameters are then optimized to reproduce the experimental data which will be described further in the next sections.

# 3. Method

## 3.1. Description of the simulated scenario

In [21], a W(110) sample with three different coverages of pre-adsorbed O is exposed to $D_2$ molecules or $D_2$ molecules and $D_2^+$ ions. The parameters of the exposure to both species are reported in table 4. For the three different O coverages, the simulation scenario contains three steps:

(i) The exposure to external sources of D (either $D_2$ only or both $D_2$ and $D_2^+$) lasts $t_{imp} = 3000$ s at the surface temperature $T_{imp}$, the pressure is $p_{D_2}$ and the incidence flux (if $D_2^+$ are present) is $\phi_{inc}$ (see table 4).

(ii) The storage: the sample resides in vacuum for $t_{storage} = 1$ h at $T_{storage}$.

(iii) TDS is realized up to 800 K with a ramp of 5 K s⁻¹.

<table><caption>**Table 4.** First row: experimental parameters for the exposure to $D_2$ (only $p_{D_2}$), implantation ($p_{D_2}$ and $\phi_{inc}$) and storage (between the exposure and the TDS) in Dunand et al [21]. Second row: parameters calculated for the MHIMS simulations.</caption><tbody><tr><td rowspan="10">[21]</td><td>p<sub>D<sub>2</sub></sub></td><td>2 &times; 10<sup>-5</sup> Pa</td></tr><tr><td>&phi;<sub>inc</sub></td><td>0.5 &times; 10<sup>18</sup> D<sub>2</sub><sup>+</sup> m<sup>-2</sup>s<sup>-1</sup></td></tr><tr><td>E<sub>inc</sub></td><td>250 eV D<sup>-1</sup></td></tr><tr><td>&alpha;<sub>inc</sub></td><td>normal incidence</td></tr><tr><td>fluence</td><td>3 &times; 10<sup>21</sup> Dm<sup>-2</sup></td></tr><tr><td>T<sub>imp</sub></td><td>300 K</td></tr><tr><td>t<sub>storage</sub></td><td>1 h</td></tr><tr><td>T<sub>storage</sub></td><td>300 K</td></tr><tr><td>r(E<sub>inc</sub>, &alpha;<sub>inc</sub>)</td><td>0.539</td></tr><tr><td>&phi;<sub>imp</sub></td><td>4.61 &times; 10<sup>17</sup> Dm<sup>-2</sup>s<sup>-1</sup></td></tr><tr><td>&Gamma;<sub>D<sub>2</sub></sub></td><td>1.52 &times; 10<sup>18</sup> D<sub>2</sub>m<sup>-2</sup>s<sup>-1</sup></td></tr><tr><td>t<sub>imp</sub></td><td>3000 s</td></tr></tbody></table>

## 3.2. Modelling steps

To identify the relevant model parameters, simulation of the experimental data from Dunand et al [21] are performed. These data being presented as raw spectrometer signal in [21], a quantitative evaluation is made which is presented in appendix A. The parametrization is performed through an optimization routine given in appendix B and inspired from [9, 46].

The parameters that one need to optimize can be divided in three groups:

1. The surface parameters $P_{sim}^{surf} = \{E_0, \Delta E, \theta_{D0}, \delta\theta_D, \beta\}$ which set the evolution of the surface energies $E_{des}(\theta_D)$ with equation (9). One set of parameters has to be determined for each surface (clean, $O_{0.50ML}$ and $O_{0.75ML}$).

2. The trap parameters $P_{sim}^{trap,i} = \{E_{dt,i}, n_i\}$ and the number of traps. We assume that the oxygen coverage does not impact the bulk trap parameters, hence only one set of trap parameters has to be determined.

3. The sputtering yields of adsorbed deuterium $P_{sim}^{sput} = \{Y_{sput}\}$ for each surface (clean, $O_{0.50ML}$ and $O_{0.75ML}$).

First, the surface parameters are determined using the results of exposure to $D_2$ for which D atoms will not enter the bulk. Then, the trap parameters are determined using the results of exposure to $D_2$ and $D_2^+$ since this is the only way to insert D atoms in the bulk traps at 300 K. According to Dunand et al [21], the sputtering process seems to have less impact during the $D_2 + D_2^+$ exposure of the $O_{0.75ML}$ surface. Thus, the trap parameters will be determined for the $O_{0.75ML}$ surface exposed to $D_2$ and $D_2^+$. For this case, the simulations revealed that $Y_{sput}^{O_{0.75ML}}$ cannot be determined unambiguously and it will be varied from 0.000 to 3.000 and an optimized set of trap parameters will be obtained for each values of $Y_{sput}^{O_{0.75ML}}$. Finally, the sputtering yields for the clean surface and the $O_{0.50ML}$ surface will be determined for each set of trap parameters (hence for different values of $Y_{sput}^{O_{0.75ML}}$).
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Comparison between the experimental and simulated TDS spectra after exposure to D₂ only for (a) the clean surface, (b) the O₀.₅₀ML surface and (c) the O₀.₇₅ML surface. The dashed line are the simulated spectra obtained with Eₐ₇₅(θₐ) derived from the DFT calculations and the solid lines are the results of the optimization on the surface parameters (solid line in figure 2).</figcaption></figure>

# 4. Results

## 4.1. Surface parameters

Figure 3 shows the results of the optimization of the surface parameters, done with the TDS spectra after exposure to D₂ only, for (a) the clean surface, (b) the O₀.₅₀ML surface and (c) the O₀.₇₅ML surface. The cost function $\epsilon$ is evaluated for temperature between 300 K and 600 K to focus on the main part of the TDS spectra. For the clean surface and the O₀.₅₀ML surface, the relative error with the optimized surface parameters is below 9%. For the O₀.₇₅ML surface, the relative error is higher at 22% presumably because the actual desorption measured by the mass spectrometer is only slightly above the fluctuation of its noise signal. The evolution Eₐ₇₅(θₐ) obtained after the optimization of the surface parameters is given in solid lines in figure 2 and the optimized parameter are reported in table 5. For the clean surface, the optimized desorption energy at the peak of the TDS spectrum is 1.20 eV and the high temperature shoulder corresponds to desorption energies in the 1.3–1.4 eV range. For the oxygen covered surfaces, the optimized desorption energies at the peak of the TDS spectra

<table><caption>**Table 5.** Optimized surface parameters P<sub>sim</sub><sup>surf</sup> = {E₀, ΔE, θ<sub>D₀</sub>, δθ<sub>D</sub>, β} to reproduce the experimental TDS spectra (see figure 3). α is calculated using equation (11) using E<sub>oversat</sub><sup>DFT</sup> equal to 0.8 eV for the clean surface [17] and 0.6 eV for the surface with adsorbed O [19].</caption><thead><tr><th></th><th>clean</th><th>O<sub>0.50ML</sub></th><th>O<sub>0.75ML</sub></th></tr></thead><tbody><tr><td>E<sub>0</sub> (eV)</td><td>1.142</td><td>1.111</td><td>1.066</td></tr><tr><td>ΔE (eV)</td><td>0.346</td><td>0.289</td><td>0.234</td></tr><tr><td>θ<sub>D<sub>0</sub></sub></td><td>0.253</td><td>0.113</td><td>0.161</td></tr><tr><td>δθ<sub>D</sub></td><td>0.180</td><td>0.082</td><td>0.057</td></tr><tr><td>β</td><td>8.902</td><td>7.240</td><td>4.144</td></tr><tr><td>α (equation (11))</td><td>0.303</td><td>0.460</td><td>0.437</td></tr></tbody></table>

<table><caption>**Table 6.** Optimized value of the trapping parameters and the sputtering yield of the clean surface and the O₀.₅₀ML surface assuming Y<sub>sput</sub><sup>O₀.₇₅ML</sup> = 1.500.</caption><tr><td>E<sub>dt,1</sub></td><td>0.912 eV</td></tr><tr><td>n<sub>1</sub></td><td>7.70 × 10<sup>-6</sup> at.fr.</td></tr><tr><td>E<sub>dt,2</sub></td><td>1.176 eV</td></tr><tr><td>n<sub>2</sub></td><td>7.90 × 10<sup>-7</sup> at.fr.</td></tr><tr><td>Y<sub>sput</sub><sup>clean</sup></td><td>3.453</td></tr><tr><td>Y<sub>sput</sub><sup>O<sub>0.50ML</sub></sup></td><td>5.820</td></tr><tr><td>Y<sub>sput</sub><sup>O<sub>0.75ML</sub></sup></td><td>1.500</td></tr></table>

are 1.10 eV and 1.05 eV for the O₀.₅₀ML and O₀.₇₅ML surfaces respectively.

## 4.2. Bulk parameters

The bulk parameters are determined with the data for the O₀.₇₅ML surface exposed to D₂ and D₂⁺. Two traps are required to reproduce the experimental TDS spectrum for this surface and exposure conditions. For each trap the detrapping energy E<sub>dt,i</sub> and the trap concentration n<sub>i</sub> are optimized for a given value of Y<sub>sput</sub> between 0.000 and 3.000. For each values of Y<sub>sput</sub>, the trapping parameters have been successfully converged to reach a relative error below 10%. As a set of satisfactory trapping parameters can be obtained for each value Y<sub>sput</sub><sup>O₀.₇₅ML</sup>, we arbitrarily chose the ones for Y<sub>sput</sub><sup>O₀.₇₅ML</sup> = 1.500 as reference values which are reported in table 6. The relative evolution of each trapping parameters compared to these reference values are shown in figure 4.

Regarding the sensitivity of the optimized detrapping energies with respect to the value of Y<sub>sput</sub><sup>O₀.₇₅ML</sup>, the variations are similar for both traps and are small: there is a relative difference from -2% to +0.5% which represents a maximum of 0.025 eV between the extrema. Thus, one can consider the detrapping energies of the two traps to be 0.91 eV and 1.18 eV. Regarding the sensitivity of the optimized trap concentration, the variations are also the same for both traps but the amplitude of the variations is much different compared to the detrapping energies: it goes from -30% to +10%. Indeed, increasing the sputtering yield reduced significantly the amount of mobile
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Comparison between the experimental and simulated TDS spectra after exposure to D₂ only for (a) the clean surface, (b) the O₀.₅₀ML surface and (c) the O₀.₇₅ML surface. The dashed line are the simulated spectra obtained with *E*<sub>des</sub>(*θ*<sub>D</sub>) derived from the DFT calculations and the solid lines are the results of the optimization on the surface parameters (solid line in figure 2).</figcaption></figure>

# 4. Results

## 4.1. Surface parameters

Figure 3 shows the results of the optimization of the surface parameters, done with the TDS spectra after exposure to D₂ only, for (a) the clean surface, (b) the O₀.₅₀ML surface and (c) the O₀.₇₅ML surface. The cost function $\epsilon$ is evaluated for temperature between 300 K and 600 K to focus on the main part of the TDS spectra. For the clean surface and the O₀.₅₀ML surface, the relative error with the optimized surface parameters is below 9%. For the O₀.₇₅ML surface, the relative error is higher at 22% presumably because the actual desorption measured by the mass spectrometer is only slightly above the fluctuation of its noise signal. The evolution $E_{des}(\theta_D)$ obtained after the optimization of the surface parameters is given in solid lines in figure 2 and the optimized parameter are reported in table 5. For the clean surface, the optimized desorption energy at the peak of the TDS spectrum is 1.20 eV and the high temperature shoulder corresponds to desorption energies in the 1.3–1.4 eV range. For the oxygen covered surfaces, the optimized desorption energies at the peak of the TDS spectra

<table><caption>**Table 5.** Optimized surface parameters $P_{sim}^{surf} = \{E_0, \Delta E, \theta_{D0}, \delta\theta_D, \beta\}$ to reproduce the experimental TDS spectra (see figure 3). $\alpha$ is calculated using equation (11) using $E_{oversat}^{DFT}$ equal to 0.8 eV for the clean surface [17] and 0.6 eV for the surface with adsorbed O [19].</caption><thead><tr><th></th><th>clean</th><th>O<sub>0.50ML</sub></th><th>O<sub>0.75ML</sub></th></tr></thead><tbody><tr><td>E<sub>0</sub> (eV)</td><td>1.142</td><td>1.111</td><td>1.066</td></tr><tr><td>ΔE (eV)</td><td>0.346</td><td>0.289</td><td>0.234</td></tr><tr><td>θ<sub>D0</sub></td><td>0.253</td><td>0.113</td><td>0.161</td></tr><tr><td>δθ<sub>D</sub></td><td>0.180</td><td>0.082</td><td>0.057</td></tr><tr><td>β</td><td>8.902</td><td>7.240</td><td>4.144</td></tr><tr><td>α (equation (11))</td><td>0.303</td><td>0.460</td><td>0.437</td></tr></tbody></table>

<table><caption>**Table 6.** Optimized value of the trapping parameters and the sputtering yield of the clean surface and the O₀.₅₀ML surface assuming $Y_{sput}^{O_{0.75ML}} = 1.500$.</caption><tr><td>E<sub>dt,1</sub></td><td>0.912 eV</td></tr><tr><td>n<sub>1</sub></td><td>7.70 × 10<sup>-6</sup> at.fr.</td></tr><tr><td>E<sub>dt,2</sub></td><td>1.176 eV</td></tr><tr><td>n<sub>2</sub></td><td>7.90 × 10<sup>-7</sup> at.fr.</td></tr><tr><td>Y<sup>clean</sup><sub>sput</sub></td><td>3.453</td></tr><tr><td>Y<sup>O<sub>0.50ML</sub></sup><sub>sput</sub></td><td>5.820</td></tr><tr><td>Y<sup>O<sub>0.75ML</sub></sup><sub>sput</sub></td><td>1.500</td></tr></table>

are 1.10 eV and 1.05 eV for the O₀.₅₀ML and O₀.₇₅ML surfaces respectively.

## 4.2. Bulk parameters

The bulk parameters are determined with the data for the O₀.₇₅ML surface exposed to D₂ and D₂⁺. Two traps are required to reproduce the experimental TDS spectrum for this surface and exposure conditions. For each trap the detrapping energy $E_{dt,i}$ and the trap concentration $n_i$ are optimized for a given value of $Y_{sput}$ between 0.000 and 3.000. For each values of $Y_{sput}$, the trapping parameters have been successfully converged to reach a relative error below 10%. As a set of satisfactory trapping parameters can be obtained for each value $Y_{sput}^{O_{0.75ML}}$, we arbitrarily chose the ones for $Y_{sput}^{O_{0.75ML}} = 1.500$ as reference values which are reported in table 6. The relative evolution of each trapping parameters compared to these reference values are shown in figure 4.

Regarding the sensitivity of the optimized detrapping energies with respect to the value of $Y_{sput}^{O_{0.75ML}}$, the variations are similar for both traps and are small: there is a relative difference from -2% to +0.5% which represents a maximum of 0.025 eV between the extrema. Thus, one can consider the detrapping energies of the two traps to be 0.91 eV and 1.18 eV. Regarding the sensitivity of the optimized trap concentration, the variations are also the same for both traps but the amplitude of the variations is much different compared to the detrapping energies: it goes from -30% to +10%. Indeed, increasing the sputtering yield reduced significantly the amount of mobile
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Sensitivity study of the optimized value of the detrapping energies $E_{dt,i}$ and the trap concentration $n_i$ with the sputtering yield $Y_{\text{sput}}^{\text{O}_{0.75\text{ML}}}$ for the O$_{0.75\text{ML}}$ surface. The reference values are taken for $Y_{\text{sput}}^{\text{O}_{0.75\text{ML}}} = 1.500$. The identification of the defect is done by comparing with atomistic data from [47]$^a$, [48]^b and [49]$^c$, see section 5.2.</figcaption></figure>

particle in the bulk, hence reducing the trapping rate: to compensate this lost of trapped particles, the trap concentration has to be increased.

### 4.3. Sputtering yield

To determine the sputtering yield for the clean surface and the O$_{0.50\text{ML}}$ surface, the trapping parameters for each value of $Y_{\text{sput}}^{\text{O}_{0.75\text{ML}}}$ are used and an optimization is done. Thus, it is not possible to perfectly reproduce the experimental TDS spectra with the variation of just one parameter. Instead, we optimized $Y_{\text{sput}}$ to match the experimental integrated retention $I_{\text{exp}} = \int_t f_{\text{exp}}(T(t))dt$. The cost function is then: $\epsilon(Y_{\text{sput}}) = \frac{|I_{\text{exp}} - I_{\text{sim}}|}{I_{\text{exp}}}$. The optimized value of $Y_{\text{sput}}^{\text{clean}}$ and $Y_{\text{sput}}^{\text{O}_{0.50\text{ML}}}$ obtained with the trapping parameters for $Y_{\text{sput}}^{\text{O}_{0.75\text{ML}}} = 1.500$ are reported in table 6 and the relative evolutions of both sputtering yields are reported in figure 5. The comparison between the experimental and optimized simulated TDS spectra are given in figure 6(b).

$Y_{\text{sput}}^{\text{clean}}$ and $Y_{\text{sput}}^{\text{O}_{0.50\text{ML}}}$ increases when $Y_{\text{sput}}^{\text{O}_{0.75\text{ML}}}$ does. Indeed, increasing $Y_{\text{sput}}^{\text{O}_{0.75\text{ML}}}$ means increasing $n_i$ (and $E_{dt,i}$). If $Y_{\text{sput}}$ stays

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Sensitivity study of the optimized value of $Y_{\text{sput}}^{\text{clean}}$ and $Y_{\text{sput}}^{\text{O}_{0.50\text{ML}}}$ with the trapping parameters obtained for value of $Y_{\text{sput}}^{\text{O}_{0.75\text{ML}}}$ between 0.000 and 3.000. The reference is set for $Y_{\text{sput}}^{\text{O}_{0.75\text{ML}}} = 1.500$.</figcaption></figure>

constant, the content of adsorbed D and thus of mobile D in the bulk would stay constant. Thus, the trapping/detrapping equilibrium ratio and consequently the proportion of filled trapped would stay constant. Thus, it would increase the retention as $n_i$ increases.

It is noticeable that the difference is wider for the O$_{0.50\text{ML}}$ surface (from -30% to +20%) than for the clean surface (from -16% to +7%). Indeed, the removal of D from the surface via the sputtering is proportional to $\theta_D$ (equation (8)). The desorption energies decreases with the coverage of O (figure 2) which means $\theta_D$ on the O$_{0.50\text{ML}}$ surface will be smaller than on the clean surface: to remove the same amount of particles with the sputtering, $Y_{\text{sput}}$ has to be higher.

## 5. Discussion

### 5.1. Surface energies

During the optimization of the surface parameters, differences between $E_{\text{des}}$ given by DFT calculations from [17, 19] and the optimized values are observed (figure 2). For the clean surface, the desorption spectrum obtained with the initial guess from DFT data of [17] has the same shape as the experimental one: a first peak at low temperature (420 K for the optimized and 470 K for the DFT desorption energies) and a high temperature shoulder. This shape is due to the increase of $E_{\text{des}}$ as the hydrogen surface coverage decreases during the TDS. However, the desorption energy distribution $E_{\text{des}}(\theta_D)$ has to be shifted downward by about 0.15 eV to match the position of the experimental temperature peak (figure 3(a)): the optimized desorption energy at the 420 K peak of the TDS spectrum is 1.20 eV and the high temperature shoulder corresponds to desorption energies in the 1.3–1.4 eV range. For the O$_{0.50\text{ML}}$ surface, the desorption energy given by DFT calculations [19]
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Sensitivity study of the optimized value of the detrapping energies $E_{dt,i}$ and the trap concentration $n_i$ with the sputtering yield $Y_{\text{sput}}^{\text{O}_{0.75\text{ML}}}$ for the O$_{0.75\text{ML}}$ surface. The reference values are taken for $Y_{\text{sput}}^{\text{O}_{0.75\text{ML}}} = 1.500$. The identification of the defect is done by comparing with atomistic data from [47]$^a$, [48]^b and [49]$^c$, see section 5.2.</figcaption></figure>

particle in the bulk, hence reducing the trapping rate: to compensate this lost of trapped particles, the trap concentration has to be increased.

### 4.3. Sputtering yield

To determine the sputtering yield for the clean surface and the O$_{0.50\text{ML}}$ surface, the trapping parameters for each value of $Y_{\text{sput}}^{\text{O}_{0.75\text{ML}}}$ are used and an optimization is done. Thus, it is not possible to perfectly reproduce the experimental TDS spectra with the variation of just one parameter. Instead, we optimized $Y_{\text{sput}}$ to match the experimental integrated retention $I_{\text{exp}} = \int_t f_{\text{exp}}(T(t))dt$. The cost function is then: $\epsilon(Y_{\text{sput}}) = \frac{|I_{\text{exp}} - I_{\text{sim}}|}{I_{\text{exp}}}$. The optimized value of $Y_{\text{sput}}^{\text{clean}}$ and $Y_{\text{sput}}^{\text{O}_{0.50\text{ML}}}$ obtained with the trapping parameters for $Y_{\text{sput}}^{\text{O}_{0.75\text{ML}}} = 1.500$ are reported in table 6 and the relative evolutions of both sputtering yields are reported in figure 5. The comparison between the experimental and optimized simulated TDS spectra are given in figure 6(b).

$Y_{\text{sput}}^{\text{clean}}$ and $Y_{\text{sput}}^{\text{O}_{0.50\text{ML}}}$ increases when $Y_{\text{sput}}^{\text{O}_{0.75\text{ML}}}$ does. Indeed, increasing $Y_{\text{sput}}^{\text{O}_{0.75\text{ML}}}$ means increasing $n_i$ (and $E_{dt,i}$). If $Y_{\text{sput}}$ stays

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Sensitivity study of the optimized value of $Y_{\text{sput}}^{\text{clean}}$ and $Y_{\text{sput}}^{\text{O}_{0.50\text{ML}}}$ with the trapping parameters obtained for value of $Y_{\text{sput}}^{\text{O}_{0.75\text{ML}}}$ between 0.000 and 3.000. The reference is set for $Y_{\text{sput}}^{\text{O}_{0.75\text{ML}}} = 1.500$.</figcaption></figure>

constant, the content of adsorbed D and thus of mobile D in the bulk would stay constant. Thus, the trapping/detrapping equilibrium ratio and consequently the proportion of filled trapped would stay constant. Thus, it would increase the retention as $n_i$ increases.

It is noticeable that the difference is wider for the O$_{0.50\text{ML}}$ surface (from -30% to +20%) than for the clean surface (from -16% to +7%). Indeed, the removal of D from the surface via the sputtering is proportional to $\theta_D$ (equation (8)). The desorption energies decreases with the coverage of O (figure 2) which means $\theta_D$ on the O$_{0.50\text{ML}}$ surface will be smaller than on the clean surface: to remove the same amount of particles with the sputtering, $Y_{\text{sput}}$ has to be higher.

## 5. Discussion

### 5.1. Surface energies

During the optimization of the surface parameters, differences between $E_{\text{des}}$ given by DFT calculations from [17, 19] and the optimized values are observed (figure 2). For the clean surface, the desorption spectrum obtained with the initial guess from DFT data of [17] has the same shape as the experimental one: a first peak at low temperature (420 K for the optimized and 470 K for the DFT desorption energies) and a high temperature shoulder. This shape is due to the increase of $E_{\text{des}}$ as the hydrogen surface coverage decreases during the TDS. However, the desorption energy distribution $E_{\text{des}}(\theta_D)$ has to be shifted downward by about 0.15 eV to match the position of the experimental temperature peak (figure 3(a)): the optimized desorption energy at the 420 K peak of the TDS spectrum is 1.20 eV and the high temperature shoulder corresponds to desorption energies in the 1.3–1.4 eV range. For the O$_{0.50\text{ML}}$ surface, the desorption energy given by DFT calculations [19]
has to be increased by about 0.1 eV at high D coverage and reaches 1.10 eV to simulate the position of the TDS peak at 400 K. At low D coverage, the desorption energy has to be increased by about 0.3 eV at low coverage up to 1.34 eV to reproduce the high temperature tail which is not present for the initial guess from DFT data. A similar observation is made for the O<sub>0.75ML</sub> surface: the desorption energy given by DFT calculations [19] has to be increased at least by about 0.2 eV and reaches 1.05 eV to reproduce the experimental data with a peak at 380 K. The high temperature desorption tail can only be reproduced by further increasing E<sub>des</sub> at low coverage up to 1.28 eV. It should be noted that with the energy calculated by DFT, the O<sub>0.75ML</sub> surface retained almost no D.

Overall, the trend observed in the optimized value is the same as the trend of the DFT data. As soon as O is adsorbed on W(110), the desorption energy of D decreases. Thus, one can imagine the co-adsorption of O with tritium could ease the tritium recovery as well as slow down the tritium uptake. It could be done by exposing tritiated W surfaces to O after the campaign to help removing tritium from surfaces. However, an important part of tritium can be trapped in the bulk that would be harder to recover with this technique. In addition, experimental observation by Whitten and Gomer [39] suggest the sticking of O on hydrogenated W(110) is drastically reduced at high H coverages. Thus, further study would be required to assess the actual kinetics of the co-adsorption and evaluate the feasibility of this technique.

What is remarkable, is that for both O coverage, E<sub>des</sub> rise up to about 1.3 eV for θ<sub>D</sub> < 0.2. This behavior possibly comes from the structure of O on the O<sub>0.50ML</sub> and O<sub>0.75ML</sub> surfaces observed in calculations [19] and experiments [21]. For the O<sub>0.50ML</sub> surface, a O p(2×1) pattern is observed for the lowest energy configuration. It creates channel of free adsorption site for H isotopes: at high coverages, these channels are filled with H and the desorption is easy but as these channels empty, the recombination is less probable as H would need to jump over these channels or diffuse multiple jumps along them. For the O<sub>0.75ML</sub> surface, the effect is even more pronounced as the channel are periodically closed by the addition of an additional O atom: thus, at low coverage, the recombination of 2 D atoms requires multiple jumps above O atoms. This could be helped by the motion of O atom on the surface creating easier configuration for the recombination of D. However, the energy barrier for this motion is above 1 eV [19, 50, 51] which can only be triggered at high temperature on the TDS spectra.

## 5.2. Nature of the identified defects

The contribution of the two traps can be identified on the TDS spectrum especially for the O<sub>0.75ML</sub> surface as the sputtering impact less the retention and TDS 5.3. As already discussed, the surface contribution of the TDS spectrum for this surface is a main peak at 380 K. Figure 6(b) shows that implanting particles contributed to the shift and growth of the main peak to 420 K which comes from the detrapping from the bulk trap 1 (E<sub>dt</sub> ≈ 0.9 eV). The high temperature shoulder from 500 K

to 600 K corresponds to the detrapping from the bulk trap 2 (E<sub>dt</sub> ≈ 1.2 eV).

These detrapping energies compare relatively well with the intrinsic trap obtained with rate equation simulations done on different W grade [7, 8, 52, 53]. According to atomistic scale calculations, detrapping energies in the range of 0.9–1.2 eV could correspond to impurities [47] or edge dislocation lines [48, 49], assuming E<sub>dt</sub> = E<sub>bind</sub> + E<sub>diff</sub>, E<sub>bind</sub> the binding energy of H with the defect.

In principle, these detrapping energies could also correspond to trapping at grain boundaries [54, 55] or mono-vacancies [32, 56]. However, the experiments have been performed on a single crystal so no grain boundaries are expected and the formation energy of mono-vacancy is to high (>3 eV [32, 56]) to have a significant amount of thermal vacancies. It could be that the ion implantation creates ion-induced vacancies in super-saturated layers [57] but the exposure condition, especially the low flux, is not prone to the creation of these ion-induced defects according to thermodynamic [58]. Secondary defect creation during D<sub>2</sub><sup>+</sup> implantation may also be considered for the creation of bulk mono vacancy near the surface of W(110) when O atoms are adsorbed on the surface [59]. A binary elastic collision between a 250 eV D<sup>-1</sup> ion and an adsorbed O atom can create an energetic 100 eV O that could further collide with a bulk W atom that would receive up to 30 eV of kinetic energy thereby creating a Frenkel pair. Some of the created vacancies could remain if the created W interstitials are accommodated as adatoms on the W(110) surface. One can note that TDS shapes for D<sub>2</sub> + D<sub>2</sub><sup>+</sup> exposed W samples are not influenced by the presence of O atoms prior to D<sub>2</sub><sup>+</sup> implantation. However, in 2023 Dunand et al have shown that multiple 2173 K annealing of the W(110) sample between experiments is needed to remove single vacancy defects [60], which was not realized in the 2022 study. Therefore, we cannot exclude that residual secondary defects could be at the origin of the observed TDS peaks for D<sub>2</sub> + D<sub>2</sub><sup>+</sup> exposure.

Previous TDS simulations [7, 52, 53] suggests that the concentration of intrinsic traps in the bulk is around 10<sup>-5</sup> – 10<sup>-3</sup> at.fr. for polycrystalline tungsten which is 2–3 orders of magnitude higher than what is obtained here. This difference is attributed to (i) the absence of grain boundaries and (ii) the thermal preparation of the SCW sample used in [21] which involves repeated heating up to 2200 K to clean the surface which presumably removed a significant amount of defects.

## 5.3. Impact of the sputtering

### 5.3.1. On the TDS spectra.
Figure 6 shows the comparison between the experimental and simulated TDS spectra after exposure to D<sub>2</sub> (a) and D<sub>2</sub> + D<sub>2</sub><sup>+</sup> (b) for the three different surfaces. To highlight the effect of the sputtering process, simulations without the sputtering process are also shown in figure 6(b). As mentioned earlier, only one parameter is changed to simulate the TDS spectra, Y<sub>sput</sub>. For the clean surface and the O<sub>0.50ML</sub> surface, the position of the TDS peak is not reproduced quantitatively, however the evolution of the
has to be increased by about 0.1 eV at high D coverage and reaches 1.10 eV to simulate the position of the TDS peak at 400 K. At low D coverage, the desorption energy has to be increased by about 0.3 eV at low coverage up to 1.34 eV to reproduce the high temperature tail which is not present for the initial guess from DFT data. A similar observation is made for the O<sub>0.75ML</sub> surface: the desorption energy given by DFT calculations [19] has to be increased at least by about 0.2 eV and reaches 1.05 eV to reproduce the experimental data with a peak at 380 K. The high temperature desorption tail can only be reproduced by further increasing E<sub>des</sub> at low coverage up to 1.28 eV. It should be noted that with the energy calculated by DFT, the O<sub>0.75ML</sub> surface retained almost no D.

Overall, the trend observed in the optimized value is the same as the trend of the DFT data. As soon as O is adsorbed on W(110), the desorption energy of D decreases. Thus, one can imagine the co-adsorption of O with tritium could ease the tritium recovery as well as slow down the tritium uptake. It could be done by exposing tritiated W surfaces to O after the campaign to help removing tritium from surfaces. However, an important part of tritium can be trapped in the bulk that would be harder to recover with this technique. In addition, experimental observation by Whitten and Gomer [39] suggest the sticking of O on hydrogenated W(110) is drastically reduced at high H coverages. Thus, further study would be required to assess the actual kinetics of the co-adsorption and evaluate the feasibility of this technique.

What is remarkable, is that for both O coverage, E<sub>des</sub> rise up to about 1.3 eV for θ<sub>D</sub> < 0.2. This behavior possibly comes from the structure of O on the O<sub>0.50ML</sub> and O<sub>0.75ML</sub> surfaces observed in calculations [19] and experiments [21]. For the O<sub>0.50ML</sub> surface, a O p(2×1) pattern is observed for the lowest energy configuration. It creates channel of free adsorption site for H isotopes: at high coverages, these channels are filled with H and the desorption is easy but as these channels empty, the recombination is less probable as H would need to jump over these channels or diffuse multiple jumps along them. For the O<sub>0.75ML</sub> surface, the effect is even more pronounced as the channel are periodically closed by the addition of an additional O atom: thus, at low coverage, the recombination of 2 D atoms requires multiple jumps above O atoms. This could be helped by the motion of O atom on the surface creating easier configuration for the recombination of D. However, the energy barrier for this motion is above 1 eV [19, 50, 51] which can only be triggered at high temperature on the TDS spectra.

## 5.2. Nature of the identified defects

The contribution of the two traps can be identified on the TDS spectrum especially for the O<sub>0.75ML</sub> surface as the sputtering impact less the retention and TDS 5.3. As already discussed, the surface contribution of the TDS spectrum for this surface is a main peak at 380 K. Figure 6(b) shows that implanting particles contributed to the shift and growth of the main peak to 420 K which comes from the detrapping from the bulk trap 1 (E<sub>dt</sub> ≈ 0.9 eV). The high temperature shoulder from 500 K

to 600 K corresponds to the detrapping from the bulk trap 2 (E<sub>dt</sub> ≈ 1.2 eV).

These detrapping energies compare relatively well with the intrinsic trap obtained with rate equation simulations done on different W grade [7, 8, 52, 53]. According to atomistic scale calculations, detrapping energies in the range of 0.9–1.2 eV could correspond to impurities [47] or edge dislocation lines [48, 49], assuming E<sub>dt</sub> = E<sub>bind</sub> + E<sub>diff</sub>, E<sub>bind</sub> the binding energy of H with the defect.

In principle, these detrapping energies could also correspond to trapping at grain boundaries [54, 55] or mono-vacancies [32, 56]. However, the experiments have been performed on a single crystal so no grain boundaries are expected and the formation energy of mono-vacancy is to high (>3 eV [32, 56]) to have a significant amount of thermal vacancies. It could be that the ion implantation creates ion-induced vacancies in super-saturated layers [57] but the exposure condition, especially the low flux, is not prone to the creation of these ion-induced defects according to thermodynamic [58]. Secondary defect creation during D<sub>2</sub><sup>+</sup> implantation may also be considered for the creation of bulk mono vacancy near the surface of W(110) when O atoms are adsorbed on the surface [59]. A binary elastic collision between a 250 eV D<sup>-1</sup> ion and an adsorbed O atom can create an energetic 100 eV O that could further collide with a bulk W atom that would receive up to 30 eV of kinetic energy thereby creating a Frenkel pair. Some of the created vacancies could remain if the created W interstitials are accommodated as adatoms on the W(110) surface. One can note that TDS shapes for D<sub>2</sub> + D<sub>2</sub><sup>+</sup> exposed W samples are not influenced by the presence of O atoms prior to D<sub>2</sub><sup>+</sup> implantation. However, in 2023 Dunand et al have shown that multiple 2173 K annealing of the W(110) sample between experiments is needed to remove single vacancy defects [60], which was not realized in the 2022 study. Therefore, we cannot exclude that residual secondary defects could be at the origin of the observed TDS peaks for D<sub>2</sub> + D<sub>2</sub><sup>+</sup> exposure.

Previous TDS simulations [7, 52, 53] suggests that the concentration of intrinsic traps in the bulk is around 10<sup>-5</sup> – 10<sup>-3</sup> at.fr. for polycrystalline tungsten which is 2–3 orders of magnitude higher than what is obtained here. This difference is attributed to (i) the absence of grain boundaries and (ii) the thermal preparation of the SCW sample used in [21] which involves repeated heating up to 2200 K to clean the surface which presumably removed a significant amount of defects.

## 5.3. Impact of the sputtering

5.3.1. On the TDS spectra. Figure 6 shows the comparison between the experimental and simulated TDS spectra after exposure to D<sub>2</sub> (a) and D<sub>2</sub> + D<sub>2</sub><sup>+</sup> (b) for the three different surfaces. To highlight the effect of the sputtering process, simulations without the sputtering process are also shown in figure 6(b). As mentioned earlier, only one parameter is changed to simulate the TDS spectra, Y<sub>sput</sub>. For the clean surface and the O<sub>0.50ML</sub> surface, the position of the TDS peak is not reproduced quantitatively, however the evolution of the
<figure><img src="image_8.png" /><figcaption>**Figure 6.** comparison between the experimental and simulated TDS spectra after exposure to $D_2$ (a) and after exposure to $D_2 + D_2^+$ (b). The solid line are simulated spectra with the sputtering taken into account, the dashed line are simulated spectra without sputtering taken into account and the symbols are the experimental data. The trapping and sputtering parameters of these simulations are reported in table 6 ($Y_{\text{sput}}^{\text{O}_{0.75\text{ML}}} = 1.500$).</figcaption></figure>

total amount of desorbed D with increasing O coverage is reproduced almost quantitatively.

When the sputtering process is not activated in the simulation, the TDS spectra are more intense as the surface is not partially depleted from its content. When the sputtering process is accounted for, the lowest temperature part of the TDS spectra for the clean and $O_{0.50\text{ML}}$ surfaces is removed. For the $O_{0.75\text{ML}}$ surface, the effect of the sputtering is much less pronounced indicating that the sputtering is not very efficient for this surface. Indeed, $Y_{\text{sput}}$ is lower for this surface and the desorption energies being low, the coverage is low diminishing the sputtered flux (equation (8)).

In order to differentiate the D retention in the bulk from the D retention on the surface, Dunand et al [21] subtracted the TDS signal obtained in both exposure conditions leading to a so-called differential retention. This differential retention is reported in figure 7 and the integrated differential retention is reported in table 7. In the simulation, when the sputtering is turned off, the (integrated) differential retention is very similar for all surfaces and the integrated differential retention corresponds to the bulk retention: thus, if there was no sputtering, the differential retention would give indeed an idea of

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Differential retention ($(D_2 + D_2^+)$ – (D$_2$)) of TDS spectra shown in figure 6.</figcaption></figure>

the bulk retention. However, the experimental results do not show this trend and the differential retention actually decreases and when diminishing the O coverage, even becomes negative for the clean surface: there is less retention when exposing to both $D_2$ and $D_2^+$. The same behavior is observed in the simulation with the sputtering process turned on. Note that, the integrated differential retention is in the same order of magnitude for the simulation and the experiments, the differences being attributed to the cumulative error resulting from the different optimization steps (surface, bulk and sputtering parameters). The simulated differentiated TDS spectra are reproducing only qualitatively the experimental ones for the clean and $O_{0.50\text{ML}}$ surfaces. Indeed, the negative peaks of the differential TDS appear at a lower temperature in the simulation. These negative peaks of the differential TDS appear because the sputtering prevents the surface coverage to build up and reach high $\theta_D$ and thus low $E_{\text{des}}$. It shows that adding the sputtering process, or any similar process removing D from the surface by the incident flux of ions, allows to reproduce the trend observed in the experiments.

### 5.3.2. On the retention dynamics.
In the previous section, we showed that the sputtering process could explain the differential retention behavior observed in the experiment. To investigate more the effect of the sputtering on the D retention dynamics, we look at the evolution of the D retention during the exposure and the storage phase. Figures 8 and 9 show the evolution of the amount of D on the surface (a), in the bulk (b) and the total amount of D (c) for the clean and $O_{0.75\text{ML}}$ surfaces respectively. For each surface, simulations with and without the sputtering process are shown.

When the sputtering process is turned off, the surface is close to saturation (shown by the shadowed area) for both surfaces during the exposure leading to low value of $E_{\text{des}}$. Thus, during the storage phase at room temperature, desorption takes place on top of the detrapping from the trap in the bulk: the total amount of D decreases. It is worth noting that the retention in the bulk is very similar in both case (about $4 \times 10^{18} \text{Dm}^{-2}$) as shown in table 7. Since the trapping parameters
<figure><img src="image_10.png" /><figcaption>**Figure 6.** comparison between the experimental and simulated TDS spectra after exposure to $D_2$ (a) and after exposure to $D_2 + D_2^+$ (b). The solid line are simulated spectra with the sputtering taken into account, the dashed line are simulated spectra without sputtering taken into account and the symbols are the experimental data. The trapping and sputtering parameters of these simulations are reported in table 6 ($Y_{\text{sput}}^{\text{O}_{0.75\text{ML}}} = 1.500$).</figcaption></figure>

total amount of desorbed D with increasing O coverage is reproduced almost quantitatively.

When the sputtering process is not activated in the simulation, the TDS spectra are more intense as the surface is not partially depleted from its content. When the sputtering process is accounted for, the lowest temperature part of the TDS spectra for the clean and $O_{0.50\text{ML}}$ surfaces is removed. For the $O_{0.75\text{ML}}$ surface, the effect of the sputtering is much less pronounced indicating that the sputtering is not very efficient for this surface. Indeed, $Y_{\text{sput}}$ is lower for this surface and the desorption energies being low, the coverage is low diminishing the sputtered flux (equation (8)).

In order to differentiate the D retention in the bulk from the D retention on the surface, Dunand et al [21] subtracted the TDS signal obtained in both exposure conditions leading to a so-called differential retention. This differential retention is reported in figure 7 and the integrated differential retention is reported in table 7. In the simulation, when the sputtering is turned off, the (integrated) differential retention is very similar for all surfaces and the integrated differential retention corresponds to the bulk retention: thus, if there was no sputtering, the differential retention would give indeed an idea of

<figure><img src="image_11.png" /><figcaption>**Figure 7.** Differential retention ($(D_2 + D_2^+) - (D_2)$) of TDS spectra shown in figure 6.</figcaption></figure>

the bulk retention. However, the experimental results do not show this trend and the differential retention actually decreases and when diminishing the O coverage, even becomes negative for the clean surface: there is less retention when exposing to both $D_2$ and $D_2^+$. The same behavior is observed in the simulation with the sputtering process turned on. Note that, the integrated differential retention is in the same order of magnitude for the simulation and the experiments, the differences being attributed to the cumulative error resulting from the different optimization steps (surface, bulk and sputtering parameters). The simulated differentiated TDS spectra are reproducing only qualitatively the experimental ones for the clean and $O_{0.50\text{ML}}$ surfaces. Indeed, the negative peaks of the differential TDS appear at a lower temperature in the simulation. These negative peaks of the differential TDS appear because the sputtering prevents the surface coverage to build up and reach high $\theta_D$ and thus low $E_{\text{des}}$. It shows that adding the sputtering process, or any similar process removing D from the surface by the incident flux of ions, allows to reproduce the trend observed in the experiments.

### 5.3.2. On the retention dynamics.
In the previous section, we showed that the sputtering process could explain the differential retention behavior observed in the experiment. To investigate more the effect of the sputtering on the D retention dynamics, we look at the evolution of the D retention during the exposure and the storage phase. Figures 8 and 9 show the evolution of the amount of D on the surface (a), in the bulk (b) and the total amount of D (c) for the clean and $O_{0.75\text{ML}}$ surfaces respectively. For each surface, simulations with and without the sputtering process are shown.

When the sputtering process is turned off, the surface is close to saturation (shown by the shadowed area) for both surfaces during the exposure leading to low value of $E_{\text{des}}$. Thus, during the storage phase at room temperature, desorption takes place on top of the detrapping from the trap in the bulk: the total amount of D decreases. It is worth noting that the retention in the bulk is very similar in both case (about $4 \times 10^{18} \text{Dm}^{-2}$) as shown in table 7. Since the trapping parameters
<table><caption>**Table 7.** Integrated differential retention (Dm⁻²) for the three different surfaces.</caption><thead><tr><th></th><th>without sputtering</th><th>with sputtering</th><th>Experiment</th></tr></thead><tbody><tr><td>W(110) clean</td><td>+4.67 × 10<sup>18</sup></td><td>-2.72 × 10<sup>18</sup></td><td>-3.38 × 10<sup>18</sup></td></tr><tr><td>W(110)O<sub>0.50ML</sub></td><td>+3.87 × 10<sup>18</sup></td><td>+4.54 × 10<sup>17</sup></td><td>+1.67 × 10<sup>17</sup></td></tr><tr><td>W(110)O<sub>0.50ML</sub></td><td>+4.05 × 10<sup>18</sup></td><td>+3.21 × 10<sup>18</sup></td><td>+3.30 × 10<sup>18</sup></td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** For the clean surface: time evolution of (a) the amount of D on the surface, (b) the amount of D in the bulk, (c) the total amount of D during the D₂ + D₂⁺ exposure and the storage phase. In (a), the shadowed area shows the saturation of the surface ($\theta_D$ = 1). The conditions during the exposure phase (0–3000 s) and the storage phase (>3000 s) are the one described in table 4. The trapping and sputtering parameters of these simulations are reported in table 6. The sputtering yields is $Y_{\text{sput}(\text{clean})} = 3.453$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** For the O₀.₇₅ML surface: time evolution of (a) the amount of D on the surface, (b) the amount of D in the bulk, (c) the total amount of D during the D₂ + D₂⁺ exposure and the storage phase. In (a), the shadowed area shows the saturation of the surface ($\theta_D$ = 1). The conditions during the exposure phase (0–3000 s) and the storage phase (>3000 s) are the one described in table 4. The trapping and sputtering parameters of these simulations are reported in table 6. The sputtering yields is $Y_{\text{sput}(\text{O}_0.75\text{ML})} = 1.500$.</figcaption></figure>

are identical, the slight difference comes from the quantity of mobile particles in the material which is dependent on the surface state.

When the sputtering process is turned on, there is a significant difference between both surfaces. For the clean surface, the sputtering is very efficient as it halves the surface retention (and the bulk one) at the end of the exposure. During the storage at room temperature, detrapping from the bulk occurs

which increases the surface retention. However, the sputtering process does not allow high coverage/low $E_{\text{des}}$ and the desorption energy during the storage step is about 1.23 eV: no desorption is observed from the clean surface and the total retention almost does not change. A very similar trend is obtained for the O₀.₅₀ML surface (but not shown here). For the O₀.₇₅ML surface, the dynamics is quite different as the desorption energies are much lower. During the exposure step,
<table><caption>**Table 7.** Integrated differential retention (Dm⁻²) for the three different surfaces.</caption><thead><tr><th></th><th>without sputtering</th><th>with sputtering</th><th>Experiment</th></tr></thead><tbody><tr><td>W(110) clean</td><td>+4.67 × 10<sup>18</sup></td><td>-2.72 × 10<sup>18</sup></td><td>-3.38 × 10<sup>18</sup></td></tr><tr><td>W(110)O<sub>0.50ML</sub></td><td>+3.87 × 10<sup>18</sup></td><td>+4.54 × 10<sup>17</sup></td><td>+1.67 × 10<sup>17</sup></td></tr><tr><td>W(110)O<sub>0.50ML</sub></td><td>+4.05 × 10<sup>18</sup></td><td>+3.21 × 10<sup>18</sup></td><td>+3.30 × 10<sup>18</sup></td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Figure 8.** For the clean surface: time evolution of (a) the amount of D on the surface, (b) the amount of D in the bulk, (c) the total amount of D during the D₂ + D₂⁺ exposure and the storage phase. In (a), the shadowed area shows the saturation of the surface ($\theta_D$ = 1). The conditions during the exposure phase (0–3000 s) and the storage phase (>3000 s) are the one described in table 4. The trapping and sputtering parameters of these simulations are reported in table 6. The sputtering yields is $Y_{\text{sput}}(\text{clean}) = 3.453$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 9.** For the O₀.₇₅ML surface: time evolution of (a) the amount of D on the surface, (b) the amount of D in the bulk, (c) the total amount of D during the D₂ + D₂⁺ exposure and the storage phase. In (a), the shadowed area shows the saturation of the surface ($\theta_D$ = 1). The conditions during the exposure phase (0–3000 s) and the storage phase (>3000 s) are the one described in table 4. The trapping and sputtering parameters of these simulations are reported in table 6. The sputtering yields is $Y_{\text{sput}}^{\text{O}_0.\text{75ML}} = 1.500$.</figcaption></figure>

are identical, the slight difference comes from the quantity of mobile particles in the material which is dependent on the surface state.

When the sputtering process is turned on, there is a significant difference between both surfaces. For the clean surface, the sputtering is very efficient as it halves the surface retention (and the bulk one) at the end of the exposure. During the storage at room temperature, detrapping from the bulk occurs

which increases the surface retention. However, the sputtering process does not allow high coverage/low $E_{\text{des}}$ and the desorption energy during the storage step is about 1.23 eV: no desorption is observed from the clean surface and the total retention almost does not change. A very similar trend is obtained for the O₀.₅₀ML surface (but not shown here). For the O₀.₇₅ML surface, the dynamics is quite different as the desorption energies are much lower. During the exposure step,
the surface is also close to the saturation and the desorption energy is 0.95 eV: as soon as the storage step starts, the desorption takes place from the surface and adding the detrapping from the bulk leads to a similar decrease of the total retention as without sputtering. Thus, the sputtering process is very impactful for the clean surface (and for the $O_{0.50ML}$ surface) while it does not impact much the inventory dynamics for the $O_{0.75ML}$ surface.

The effect of the sputtering process for the clean surface could be highlighted because the exposure time was short: thus the surface retention is dominating the total retention (75% in figure 8). The fact that the sputtering process impact less the $O_{0.75ML}$ surface is also due to the smaller amount of available site on the surface: the surface retention for this surface is only 30% (figure 9). Thus, the impact of the sputtering process is the addition of high $E_{des}$ and high surface retention (high $n_{surf}$ and low bulk retention) which explains why the sputtering is more efficient for the clean surface. It also implies that such process would be difficult to highlight with higher fluence irradiation that would increase the bulk retention.

### 5.3.3. On the fuel recycling.
In order to have a quantitative agreement with the integrated differential retention, sputtering yields between 1 and 12 have to be used. It means that about 1 to $12 \times 10^{19} \text{ Dm}^{-2}\text{s}^{-1}$ can be sputtered from the W surface to the plasma and recycled if the surface is saturated. Thus, the flux of D from the surface to the vacuum during the $D_2 + D_2^+$ exposure has another components. To quantify how much this source impacts the D recycling from the surface during the exposure, we introduce the fuel recycling coefficient defined here as:

$$R_D = \frac{\varphi_{des} + \varphi_{sput}}{\varphi_{imp} + \varphi_{gas}}. \qquad (12)$$

It gives how much D atoms are removed from the system (from desorption and sputtering) per D atoms that are introduced in the system (through implantation and sticking).

Figure 10 shows the recycling coefficient $R_D$ for the clean surface and the parameters of table 6 in the simulation with sputtering (solid line) and without sputtering (dashed line). In both cases, the recycling coefficient tends toward 1 expressing that almost all particle introduced in the simulated material are released toward the vacuum/plasma after a sufficient hydrogen fluence.

Without sputtering, different characteristic times can be observed. Such behavior has already been reported in previous studies investigating fuel recycling with rate equation simulations [61, 62]. The first characteristic time is due to the build up of the subsurface inventory, especially the growth of $c_m$ in the implantation zone [62]. In our simulations, the build up of the surface inventory also play an important role as increasing $c_{surf}$ changing $E_{des}$ (figure 2). Once stabilized, the total amount of D grows thanks to migration in the bulk which is slow and evolves with the square root of times [11, 24, 61] as observed in figure 8. In that case, the recycling is only through the desorption of $D_2$ as there is no sputtering activated.

<figure><img src="image_12.png" /><figcaption>Figure 10. Recycling coefficient $R_D$ for the clean surface case during exposure with $D_2 + D_2^+$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 11. Recycling coefficient at the end of the implantation $R_D(t_{imp})$ for the clean surface case during exposure with $D_2 + D_2^+$ for different implantation temperature $T_{imp}$ and incident fluxes from $1 \times 10^{18} \text{ Dm}^{-2}\text{s}^{-1}$ to $1 \times 10^{24} \text{ Dm}^{-2}\text{s}^{-1}$. All the other exposure parameters ($P_{D_2}$, $E_{inc}$, $\theta_{inc}$) are the one reported in table 4.</figcaption></figure>

When the sputtering is activated, the recycling is dominated by the sputtering and it reaches 1 more quickly, explaining also the reduced total amount of retained D atoms. The negligible recycling from $D_2$ desorption comes also from the high $E_{des}$ of the clean surface during the implantation (about 1.23 eV). However, even with the $O_{0.75ML}$ surface, the recycling of D as $D_2$ is only $10^{-3}$.

In these simulations, the implantation is done at 300 K at low $D_2^+$ flux and low $D_2$ pressure which reduces the flux of $D_2$ desorption. The effect of the temperature $T_{imp}$ and the incident flux $\phi_{inc}$ on the processes that contribute to the recycling coefficient is shown in figure 11, keeping the other exposure parameters equal to the values reported in table 4. Because the desorption of $D_2$ is thermally activated, the recycling at low
the surface is also close to the saturation and the desorption energy is 0.95 eV: as soon as the storage step starts, the desorption takes place from the surface and adding the detrapping from the bulk leads to a similar decrease of the total retention as without sputtering. Thus, the sputtering process is very impactful for the clean surface (and for the $O_{0.50ML}$ surface) while it does not impact much the inventory dynamics for the $O_{0.75ML}$ surface.

The effect of the sputtering process for the clean surface could be highlighted because the exposure time was short: thus the surface retention is dominating the total retention (75% in figure 8). The fact that the sputtering process impact less the $O_{0.75ML}$ surface is also due to the smaller amount of available site on the surface: the surface retention for this surface is only 30% (figure 9). Thus, the impact of the sputtering process is the addition of high $E_{des}$ and high surface retention (high $n_{surf}$ and low bulk retention) which explains why the sputtering is more efficient for the clean surface. It also implies that such process would be difficult to highlight with higher fluence irradiation that would increase the bulk retention.

### 5.3.3. On the fuel recycling.
In order to have a quantitative agreement with the integrated differential retention, sputtering yields between 1 and 12 have to be used. It means that about 1 to $12 \times 10^{19} \text{ Dm}^{-2}\text{s}^{-1}$ can be sputtered from the W surface to the plasma and recycled if the surface is saturated. Thus, the flux of D from the surface to the vacuum during the $D_2 + D_2^+$ exposure has another components. To quantify how much this source impacts the D recycling from the surface during the exposure, we introduce the fuel recycling coefficient defined here as:

$$R_D = \frac{\varphi_{des} + \varphi_{sput}}{\varphi_{imp} + \varphi_{gas}}. \qquad (12)$$

It gives how much D atoms are removed from the system (from desorption and sputtering) per D atoms that are introduced in the system (through implantation and sticking).

Figure 10 shows the recycling coefficient $R_D$ for the clean surface and the parameters of table 6 in the simulation with sputtering (solid line) and without sputtering (dashed line). In both cases, the recycling coefficient tends toward 1 expressing that almost all particle introduced in the simulated material are released toward the vacuum/plasma after a sufficient hydrogen fluence.

Without sputtering, different characteristic times can be observed. Such behavior has already been reported in previous studies investigating fuel recycling with rate equation simulations [61, 62]. The first characteristic time is due to the build up of the subsurface inventory, especially the growth of $c_m$ in the implantation zone [62]. In our simulations, the build up of the surface inventory also play an important role as increasing $c_{surf}$ changing $E_{des}$ (figure 2). Once stabilized, the total amount of D grows thanks to migration in the bulk which is slow and evolves with the square root of times [11, 24, 61] as observed in figure 8. In that case, the recycling is only through the desorption of $D_2$ as there is no sputtering activated.

<figure><img src="image_14.png" /><figcaption>Figure 10. Recycling coefficient $R_D$ for the clean surface case during exposure with $D_2 + D_2^+$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 11. Recycling coefficient at the end of the implantation $R_D(t_{imp})$ for the clean surface case during exposure with $D_2 + D_2^+$ for different implantation temperature $T_{imp}$ and incident fluxes from $1 \times 10^{18} \text{ Dm}^{-2}\text{s}^{-1}$ to $1 \times 10^{24} \text{ Dm}^{-2}\text{s}^{-1}$. All the other exposure parameters ($P_{D_2}$, $E_{inc}$, $\theta_{inc}$) are the one reported in table 4.</figcaption></figure>

When the sputtering is activated, the recycling is dominated by the sputtering and it reaches 1 more quickly, explaining also the reduced total amount of retained D atoms. The negligible recycling from $D_2$ desorption comes also from the high $E_{des}$ of the clean surface during the implantation (about 1.23 eV). However, even with the $O_{0.75ML}$ surface, the recycling of D as $D_2$ is only $10^{-3}$.

In these simulations, the implantation is done at 300 K at low $D_2^+$ flux and low $D_2$ pressure which reduces the flux of $D_2$ desorption. The effect of the temperature $T_{imp}$ and the incident flux $\phi_{inc}$ on the processes that contribute to the recycling coefficient is shown in figure 11, keeping the other exposure parameters equal to the values reported in table 4. Because the desorption of $D_2$ is thermally activated, the recycling at low
temperature is dominated by the sputtering while it is dom-
inated by D₂ desorption as soon as the temperature is high
enough. Also, because the sputtering is proportional to θD and
Edes depends also on θD, the transition temperature between
both regimes shifts toward high temperatures with increasing
φinc. In the divertor of a fusion reactor, the surface temperature
can reach 1000 K and fluxes of about 10²² – 10²⁴ Dm⁻²s⁻¹
are expected [1, 2, 6]. Thus, the part of recycled fuel through
sputtering could be important especially in erosion region near
the strike points.

The above results show that the sputtering process changes significantly the nature of the recycling during ion implantation. Guterl et al [13] explained that the kinetic surface model, with φdes as the only term for recycling, can lead to unphysical regimes where the surface is saturated and φdes < φimp: with Edes ≈ 1.5 eV, the maximum fluxes desorbing from a saturated surfaces at about 500 K is φdes (500 K) ≈ 10¹⁷ Dm⁻²s⁻¹ and φdes (1000 K) ≈ 10²⁴ Dm⁻²s⁻¹. Such temperatures and incident fluxes can be encountered in the divertor region [5, 6] with implanted flux higher or similar to φdes. Thus, a model who tries to estimate the D retention in the divertor with such model may predict a linear growth if φdes < φimp, which is not observed in laboratory nor tokamak experiments. Adding first the decrease of Edes with θD and second the sputtering process allows to remove the unphysical regimes and simulate more accurately the recycling behavior of D from the PFCs.

The cross section for the Eley-Rideal recombination with
<5 eV D<sup>-1</sup> atoms is of the order of 10<sup>-21</sup> m<sup>2</sup> [63, 64] which is
taken as input in MHIMS calculations to simulate low energy
D atoms on W surface [11, 24]. The sputtering cross section
(see equation (7)) calculated from the sputtering yields given
in table 6 are around 5 × 10<sup>-19</sup> m<sup>2</sup> which is 2 order of mag-
nitude higher than the abstraction processes for low energy
atom. One still need to understand what are the products of
the sputtering process. As mentioned in section 2.2.2, the sput-
tering processes may have similarities with hot-atom recom-
bination or Eley-Rideal recombination as the reaction pro-
ceeds with adsorbed atom. In addition, the conditions for these
processes are close to the Swift Chemical Sputtering (SCS)
conditions as described in [65]: for the SCS to occur, the
energy of the incident particle need to be high enough to
break the bound with the target but low enough to stay in
the vicinity of the sputtered atom to bind with it. For Eley-
Rideal recombination to occur with ≈100 eV ions, the normal
kinetic energy of the incident atom must be low enough to
allow recombination with the adsorbed D atom, which may
be possible in the strike points region of the divertor where
grazing incidence conditions has been designed for thermal-
dissipation issues. MD simulation of plasma facing materi-
als irradiation by energetic D atoms suggests that the reflec-
tion and recycling of D through D₂ decreases with the incid-
ent energy [66]. Thus, the product of sputtering may be D
atoms with an energy distribution to be determined. The effect
of these sputtered D atoms on the edge plasma differs from
thermal D₂ released via Langmuir-Hinshelwood recombin-
ation or energetic D₂ released via Eley-Rideal or hot atom

recombination. To determine the products of sputtering and
their energy distribution, MD calculations should be done
with accurate interatomic potentials reproducing correctly the
hydrogen surface interactions. Such MD calculations would
also be helpful to parametrized more closely the sputtering
yield used in this work.

6. Conclusions

In this paper, rate equation modelling has been used to sim-
ulate experiments presented in [21], where D₂ and D₂⁺ + D₂
exposure followed by subsequent TDS were realized on a
W(110) sample with various O coverage.

The desorption energies calculated by DFT [17, 19] have
been shifted by just 0.2 eV to reproduce the experimental TDS
spectra. The surface desorption energy determined from the
TDS simulations follow the same trend as the DFT calcula-
tions: the more O is adsorbed on the surface, the easier the
desorption is. The bulk detrapping energies determined from
the simulation of D₂⁺ + D₂ exposure proved to be in the range
of previously determined detrapping energies in bulk W.

In order to reproduce qualitatively the experimental TDS
spectra, and quantitatively the integrated experimental TDS
(retention), the sputtering of D atoms by D₂⁺ has to be con-
sidered. A similar processes has been used to simulated exper-
imental TDS in steel [67] making it important to be taken into
account in surface affected desorption regimes. The sputtering
process is shown to affect significantly the shape and position
of the simulated TDS spectra, the bulk and surface retention
and also the product of the recycling from the wall to the vacu-
um/plasma.

In the experimental conditions of [21], the sputtering is the
main recycling channel and a simple extrapolation to other
exposure condition (higher flux and higher temperature) sug-
gests that in divertor conditions, the sputtering may also be
a significant recycling channel. However, for divertor con-
ditions, the energy may be tens of eV instead of 250 eV as
used here: the evolution of the sputtering yield with energy
should be evaluated. In addition, the nature of the product of
the sputtering process needs to be determined, whether it is
recombined molecules as for the Eley-Rideal and hot atom
recombination or neutral energetic atoms. MD calculations
with accurate potential for surface properties should be able
to assess these open questions.

Acknowledgments

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 opin-
ions 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.
temperature is dominated by the sputtering while it is dom-
inated by D₂ desorption as soon as the temperature is high
enough. Also, because the sputtering is proportional to θD and
Edes depends also on θD, the transition temperature between
both regimes shifts toward high temperatures with increasing
φinc. In the divertor of a fusion reactor, the surface temperature
can reach 1000 K and fluxes of about 10²² – 10²⁴ Dm⁻²s⁻¹
are expected [1, 2, 6]. Thus, the part of recycled fuel through
sputtering could be important especially in erosion region near
the strike points.

The above results show that the sputtering process changes
significantly the nature of the recycling during ion implanta-
tion. Guterl et al [13] explained that the kinetic surface model,
with φdes as the only term for recycling, can lead to unphysical
regimes where the surface is saturated and φdes < φimp: with
Edes ≈ 1.5 eV, the maximum fluxes desorbing from a saturated
surfaces at about 500 K is φdes (500 K) ≈ 10¹⁷ Dm⁻²s⁻¹ and
φdes (1000 K) ≈ 10²⁴ Dm⁻²s⁻¹. Such temperatures and incid-
ent fluxes can be encountered in the divertor region [5, 6]
with implanted flux higher or similar to φdes. Thus, a model
who tries to estimate the D retention in the divertor with such
model may predict a linear growth if φdes < φimp, which is not
observed in laboratory nor tokamak experiments. Adding first
the decrease of Edes with θD and second the sputtering process
allows to remove the unphysical regimes and simulate more
accurately the recycling behavior of D from the PFCs.

The cross section for the Eley-Rideal recombination with
<5 eV D<sup>-1</sup> atoms is of the order of 10<sup>-21</sup> m<sup>2</sup> [63, 64] which is
taken as input in MHIMS calculations to simulate low energy
D atoms on W surface [11, 24]. The sputtering cross section
(see equation (7)) calculated from the sputtering yields given
in table 6 are around 5 × 10<sup>-19</sup> m<sup>2</sup> which is 2 order of mag-
nitude higher than the abstraction processes for low energy
atom. One still need to understand what are the products of
the sputtering process. As mentioned in section 2.2.2, the sput-
tering processes may have similarities with hot-atom recom-
bination or Eley-Rideal recombination as the reaction pro-
ceeds with adsorbed atom. In addition, the conditions for these
processes are close to the Swift Chemical Sputtering (SCS)
conditions as described in [65]: for the SCS to occur, the
energy of the incident particle need to be high enough to
break the bound with the target but low enough to stay in
the vicinity of the sputtered atom to bind with it. For Eley-
Rideal recombination to occur with ≈100 eV ions, the normal
kinetic energy of the incident atom must be low enough to
allow recombination with the adsorbed D atom, which may
be possible in the strike points region of the divertor where
grazing incidence conditions has been designed for thermal-
dissipation issues. MD simulation of plasma facing materi-
als irradiation by energetic D atoms suggests that the reflec-
tion and recycling of D through D₂ decreases with the incid-
ent energy [66]. Thus, the product of sputtering may be D
atoms with an energy distribution to be determined. The effect
of these sputtered D atoms on the edge plasma differs from
thermal D₂ released via Langmuir-Hinshelwood recombin-
ation or energetic D₂ released via Eley-Rideal or hot atom

recombination. To determine the products of sputtering and
their energy distribution, MD calculations should be done
with accurate interatomic potentials reproducing correctly the
hydrogen surface interactions. Such MD calculations would
also be helpful to parametrized more closely the sputtering
yield used in this work.

6. Conclusions

In this paper, rate equation modelling has been used to sim-
ulate experiments presented in [21], where D₂ and D₂⁺ + D₂
exposure followed by subsequent TDS were realized on a
W(110) sample with various O coverage.

The desorption energies calculated by DFT [17, 19] have
been shifted by just 0.2 eV to reproduce the experimental TDS
spectra. The surface desorption energy determined from the
TDS simulations follow the same trend as the DFT calcula-
tions: the more O is adsorbed on the surface, the easier the
desorption is. The bulk detrapping energies determined from
the simulation of D₂⁺ + D₂ exposure proved to be in the range
of previously determined detrapping energies in bulk W.

In order to reproduce qualitatively the experimental TDS
spectra, and quantitatively the integrated experimental TDS
(retention), the sputtering of D atoms by D₂⁺ has to be con-
sidered. A similar processes has been used to simulated exper-
imental TDS in steel [67] making it important to be taken into
account in surface affected desorption regimes. The sputtering
process is shown to affect significantly the shape and position
of the simulated TDS spectra, the bulk and surface retention
and also the product of the recycling from the wall to the vacu-
um/plasma.

In the experimental conditions of [21], the sputtering is the
main recycling channel and a simple extrapolation to other
exposure condition (higher flux and higher temperature) sug-
gests that in divertor conditions, the sputtering may also be
a significant recycling channel. However, for divertor con-
ditions, the energy may be tens of eV instead of 250 eV as
used here: the evolution of the sputtering yield with energy
should be evaluated. In addition, the nature of the product of
the sputtering process needs to be determined, whether it is
recombined molecules as for the Eley-Rideal and hot atom
recombination or neutral energetic atoms. MD calculations
with accurate potential for surface properties should be able
to assess these open questions.

Acknowledgments

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 opin-
ions 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.
The project leading to this publication has received funding from the Excellence Initiative of Aix-Marseille University— A*Midex, a French ‘Investissements d’Avenir’ programme as well as from the ANR under Grant ANR-18-CE05-0012. Centre de Calcul Intensif d’Aix-Marseille is acknowledged for granting access to high performance computing resources.

## Appendix A. Quantitative evaluation of the experimental data

In [21], the desorption rate signals are in arbitrary units. In order to make a direct comparison between our simulations and the experimental results, the experimental TDS signals have been calibrated. This is done by using the TDS signal from the clean W(110) surface exposed to gas only.

According to a thermodynamic model parametrized by DFT and phonon calculations, and validated experimentally [38], the maximum coverage for the exposure reported in table 4 is 1 H/W which represents $1.416 \times 10^{19}$ Hm$^{-2}$. However, after the exposure to D$_2$, the sample is stored for 1 h at 300 K under vacuum and D$_2$ can desorb from the surface before the TDS. According to [39], the coverage of a clean W(110) exposed to H$_2$ at 90 K is also close to 1 H/W (0.94 H/W). Thus, the integrated part of the TDS spectra given in [39] represents 0.94 H/W. An evaluation of what would remain after an exposure at 300 K and a storage phase at 300 K can be made by integrating the TDS spectrum of [39] from $T_1$ (with $T_1 \ge 300$ K) to the final TDS temperature. As experimentally observed [68], the storage at the same temperature as the exposure can lead to desorption of retained D leading to a temperature shift toward high temperature. Thus, for a storage of several hours, one should integrate from $T_1 = 400$ K to account for a possible desorption during the storage. With this method, it is calculated that the quantity of H desorbing in the experiment of [21] after exposure and storage at 300 K can be evaluated between $1.30 \times 10^{19}$ and $1.20 \times 10^{19}$ m$^{-2}$.

(integrating from 300 K or 400 K respectively). A final value of $1.23 \times 10^{19}$ m$^{-2}$ is chosen to scale the integral of the TDS signal of the clean surface exposed to D$_2$ in [21] and thus to scale the desorption flux of all surfaces and exposure conditions.

## Appendix B. Parametric optimization

In order to obtain a quantitative comparison between experimental and simulated desorption spectra, a cost function evaluating the relative error $\epsilon$ is used:

$$ \epsilon(P_{\text{sim}}) = \frac{\sum_i |f_{\text{exp}}(T_i) - f_{\text{sim}}(P_{\text{sim}}, T_i)|}{\sum_i f_{\text{exp}}(T_i)} \quad (13) $$

where $f_{\text{exp}}(T_i)$ is the experimental desorption flux at the temperature $T_i$ and $f_{\text{sim}}(P_{\text{sim}}, T_i)$ is the simulated desorption flux at the same temperature and $P_{\text{sim}}$ is the set of parameters used for the simulation. The description of the parameters in $P_{\text{sim}}$ is presented in section 3.2.

To optimize the set of simulation parameters $P_{\text{sim}}$, the same method as presented in [9, 46] is used. Using an initial guess $P_{\text{sim}}(0)$, $P_{\text{sim}}$ is optimized in order to minimize the cost function $\epsilon$ using the `scipy.optimize.minimize` function from the `scipy` package [69]. An example is shown in figure 12 for the clean W(110) exposed to gas only which shows the comparison with the experimental and simulated spectra (a) and the evolution of the cost function with the iteration (c). The evolutions of $|f_{\text{exp}}(T_i) - f_{\text{sim}}(P_{\text{sim}}, T_i)|$ for the different iterations are also displayed in figure 12(b) to show at which temperature there are the most deviation between the simulation and the experiment. For this particular case, the parameters to optimize are the one defining the evolution of $E_{\text{des}}(\theta_D)$ (equation (9)): $P_{\text{sim}} = \{E_0, \Delta E, \theta_{D0}, \delta\theta_D, \beta\}$ ($\alpha$ is constrained by equation (11)). An optimized set of parameters is achieved after 723 evaluations of the cost function which ends up to be around 8%. The remaining error is mainly due to the small shoulder at 350 K in the experimental data (cf figure 12(a)), as one observes a peak in figure 12(b).
The project leading to this publication has received funding from the Excellence Initiative of Aix-Marseille University— A*Midex, a French ‘Investissements d’Avenir’ programme as well as from the ANR under Grant ANR-18-CE05-0012. Centre de Calcul Intensif d’Aix-Marseille is acknowledged for granting access to high performance computing resources.

## Appendix A. Quantitative evaluation of the experimental data

In [21], the desorption rate signals are in arbitrary units. In order to make a direct comparison between our simulations and the experimental results, the experimental TDS signals have been calibrated. This is done by using the TDS signal from the clean W(110) surface exposed to gas only.

According to a thermodynamic model parametrized by DFT and phonon calculations, and validated experimentally [38], the maximum coverage for the exposure reported in table 4 is 1 H/W which represents $1.416 \times 10^{19}$ Hm$^{-2}$. However, after the exposure to D$_2$, the sample is stored for 1 h at 300 K under vacuum and D$_2$ can desorb from the surface before the TDS. According to [39], the coverage of a clean W(110) exposed to H$_2$ at 90 K is also close to 1 H/W (0.94 H/W). Thus, the integrated part of the TDS spectra given in [39] represents 0.94 H/W. An evaluation of what would remain after an exposure at 300 K and a storage phase at 300 K can be made by integrating the TDS spectrum of [39] from $T_1$ (with $T_1 \ge 300$ K) to the final TDS temperature. As experimentally observed [68], the storage at the same temperature as the exposure can lead to desorption of retained D leading to a temperature shift toward high temperature. Thus, for a storage of several hours, one should integrate from $T_1 = 400$ K to account for a possible desorption during the storage. With this method, it is calculated that the quantity of H desorbing in the experiment of [21] after exposure and storage at 300 K can be evaluated between $1.30 \times 10^{19}$ and $1.20 \times 10^{19}$ m$^{-2}$

(integrating from 300 K or 400 K respectively). A final value of $1.23 \times 10^{19}$ m$^{-2}$ is chosen to scale the integral of the TDS signal of the clean surface exposed to D$_2$ in [21] and thus to scale the desorption flux of all surfaces and exposure conditions.

## Appendix B. Parametric optimization

In order to obtain a quantitative comparison between experimental and simulated desorption spectra, a cost function evaluating the relative error $\epsilon$ is used:

$$ \epsilon(P_{\text{sim}}) = \frac{\sum_i |f_{\text{exp}}(T_i) - f_{\text{sim}}(P_{\text{sim}}, T_i)|}{\sum_i f_{\text{exp}}(T_i)} \quad (13) $$

where $f_{\text{exp}}(T_i)$ is the experimental desorption flux at the temperature $T_i$ and $f_{\text{sim}}(P_{\text{sim}}, T_i)$ is the simulated desorption flux at the same temperature and $P_{\text{sim}}$ is the set of parameters used for the simulation. The description of the parameters in $P_{\text{sim}}$ is presented in section 3.2.

To optimize the set of simulation parameters $P_{\text{sim}}$, the same method as presented in [9, 46] is used. Using an initial guess $P_{\text{sim}}(0)$, $P_{\text{sim}}$ is optimized in order to minimize the cost function $\epsilon$ using the `scipypimize.minimize` function from the `scipy` package [69]. An example is shown in figure 12 for the clean W(110) exposed to gas only which shows the comparison with the experimental and simulated spectra (a) and the evolution of the cost function with the iteration (c). The evolutions of $|f_{\text{exp}}(T_i) - f_{\text{sim}}(P_{\text{sim}}, T_i)|$ for the different iterations are also displayed in figure 12(b) to show at which temperature there are the most deviation between the simulation and the experiment. For this particular case, the parameters to optimize are the one defining the evolution of $E_{\text{des}}(\theta_D)$ (equation (9)): $P_{\text{sim}} = \{E_0, \Delta E, \theta_{D0}, \delta\theta_D, \beta\}$ ($\alpha$ is constrained by equation (11)). An optimized set of parameters is achieved after 723 evaluations of the cost function which ends up to be around 8%. The remaining error is mainly due to the small shoulder at 350 K in the experimental data (cf figure 12(a)), as one observes a peak in figure 12(b).
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Illustration of the parametric optimization for the clean W(110) exposed to D₂ gas only. (a) Comparison between the experimental and simulated TDS spectra. (b) evolution of $|f_{\text{exp}}(T_i) - f_{\text{sim}}(P_{\text{sim}}, T_i)|$ during the optimization process. The spectra obtained with the initial and optimized parameters are displayed in green and red respectively. The spectra obtained with the intermediate parameters are shown in grey. (c) Evolution of the cost function $\epsilon$ with the iteration.</figcaption></figure>

## ORCID iDs

E.A. Hodille https://orcid.org/0000-0002-0859-390X

J. Denis https://orcid.org/0000-0002-1478-0269

Y. Ferro https://orcid.org/0000-0002-3178-937X

M. Minissale https://orcid.org/0000-0001-6331-1402

C. Grisolia https://orcid.org/0000-0002-3038-9593

R. Bisson https://orcid.org/0000-0002-8819-1563

## References

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

[2] Denis J., Bucalossi J., Ciraolo G., Hodille E.A., Pégourié B., Bufferand H., Grisolia C., Loarer T., Marandet Y. and Serre E. 2019 Dynamic modelling of local fuel inventory and desorption in the whole tokamak vacuum vessel for auto-consistent plasma-wall interaction simulations Nucl. Mater. Energy 19 550–7

[3] Philipps V. 2011 Tungsten as material for plasma-facing components in fusion devices J. Nucl. Mater. 415 S2–S9

[4] Grisolia C. 1999 Plasma wall interaction during long pulse operation in Tore Supra J. Nucl. Mater. 266–269 146–52

[5] Hodille E.A. et al 2021 Modelling of hydrogen isotopes trapping, diffusion and permeation in divertor monoblocks under ITER-like conditions Nucl. Fusion 61 126003

[6] Delaporte-Mathurin R. et al 2021 Fuel retention in WEST and ITER divertors based on FESTIM monoblock simulations Nucl. Fusion 61 126001

[7] Hodille E.A., Bonnin X., Bisson R., Angot T., Becquart C.S., Layet J.M. and Grisolia C. 2015 Macroscopic rate equation modeling of trapping/detrapping of hydrogen isotopes in tungsten materials J. Nucl. Mater. 467 424–31

[8] Hodille E.A. et al 2017 Retention and release of hydrogen isotopes in tungsten plasma-facing components: the role of grain boundaries and the native oxide layer from a joint experiment-simulation integrated approach Nucl. Fusion 57 076019

[9] Delaporte-Mathurin R., Hodille E.A., Mougenot J., Charles Y. and Grisolia C. 2021 Parametric optimisation based on TDS experiments for rapid and efficient identification of hydrogen transport materials properties Nucl. Mater. Energy 27 100984

[10] Markelj S., Ogorodnikova O.V., Pelicon P., Schwarz-Selinger T. and Čadež I. 2013 Temperature dependence of D atom adsorption on polycrystalline tungsten Appl. Surf. Sci. 282 478–86

[11] Hodille E.A., Založnik A., Markelj S., Schwarz-Selinger T., Becquart C.S., Bisson R. and Grisolia C. 2017 Simulations of atomic deuterium exposure in self-damaged tungsten Nucl. Fusion 57 056002

[12] Založnik A., Markelj S., Schwarz-Selinger T. and Schmid K. 2017 Deuterium atom loading of self-damaged tungsten at different sample temperatures J. Nucl. Mater. 496 1–8

[13] Guterl J., Smirnov R.D. and Snyder P. 2019 Effects of surface processes on hydrogen outgassing from metal in desorption experiments Nucl. Fusion 59 096042

[14] Schmid K. and Zibrov M. 2021 On the use of recombination rate coefficients in hydrogen transport calculations Nucl. Fusion 61 086008

[15] Hodille E.A., Payet M., Marascu V., Peillon S., Mougenot J., Ferro Y., Delaporte-Mathurin R., Leblond F., Bernard E. and Grisolia C. 2021 Modelling tritium adsorption and desorption from tungsten dust particles with a surface kinetic model Nucl. Fusion 61 086030

[16] Pick M.A. and Sonnenberg K. 1985 A model for atomic hydrogen-metal interactions—application to recycling, recombination and permeation J. Nucl. Mater. 131 208–20

[17] Ajmalghan M., Piazza Z.A., Hodille E.A. and Ferro Y. 2019 Surface coverage dependent mechanisms for the absorption and desorption of hydrogen from the W(110) and w(100) surfaces: a density functional theory investigation Nucl. Fusion 59 106022

[18] Piazza Z.A., Ajmalghan M., Ferro Y. and Kolasinski R.D. 2018 Saturation of tungsten surfaces with hydrogen: a density functional theory study complemented by low energy ion scattering and direct recoil spectroscopy data Acta Mater. 145 388–98

[19] Ferro Y., Augustin Hodille E., Denis J., Piazza Z.A. and Ajmalghan M. 2023 Hydrogen and oxygen on tungsten (110) surface: adsorption, absorption and desorption investigated by density functional theory Nucl. Fusion 63 036017

[20] Bergstrom Z.J., Li C., Samolyuk G.D., Uberuaga B.P. and Wirth B.D. 2019 Hydrogen interactions with low-index surface orientations of tungsten J. Phys.: Condens. Matter 31 255002

[21] Dunand A., Minissale M., Faure J.-B., Gallais L., Angot T. and Bisson R. 2022 Surface oxygen versus native oxide on
<figure><img src="image_16.png" /><figcaption>**Figure 12.** Illustration of the parametric optimization for the clean W(110) exposed to D₂ gas only. (a) Comparison between the experimental and simulated TDS spectra. (b) evolution of $|f_{\text{exp}}(T_i) - f_{\text{sim}}(P_{\text{sim}}, T_i)|$ during the optimization process. The spectra obtained with the initial and optimized parameters are displayed in green and red respectively. The spectra obtained with the intermediate parameters are shown in grey. (c) Evolution of the cost function $\epsilon$ with the iteration.</figcaption></figure>

## ORCID iDs

E.A. Hodille https://orcid.org/0000-0002-0859-390X

J. Denis https://orcid.org/0000-0002-1478-0269

Y. Ferro https://orcid.org/0000-0002-3178-937X

M. Minissale https://orcid.org/0000-0001-6331-1402

C. Grisolia https://orcid.org/0000-0002-3038-9593

R. Bisson https://orcid.org/0000-0002-8819-1563

## References

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

[2] Denis J., Bucalossi J., Ciraolo G., Hodille E.A., Pégourié B., Bufferand H., Grisolia C., Loarer T., Marandet Y. and Serre E. 2019 Dynamic modelling of local fuel inventory and desorption in the whole tokamak vacuum vessel for auto-consistent plasma-wall interaction simulations Nucl. Mater. Energy 19 550–7

[3] Philipps V. 2011 Tungsten as material for plasma-facing components in fusion devices J. Nucl. Mater. 415 S2–S9

[4] Grisolia C. 1999 Plasma wall interaction during long pulse operation in Tore Supra J. Nucl. Mater. 266–269 146–52

[5] Hodille E.A. et al 2021 Modelling of hydrogen isotopes trapping, diffusion and permeation in divertor monoblocks under ITER-like conditions Nucl. Fusion 61 126003

[6] Delaporte-Mathurin R. et al 2021 Fuel retention in WEST and ITER divertors based on FESTIM monoblock simulations Nucl. Fusion 61 126001

[7] Hodille E.A., Bonnin X., Bisson R., Angot T., Becquart C.S., Layet J.M. and Grisolia C. 2015 Macroscopic rate equation modeling of trapping/detrapping of hydrogen isotopes in tungsten materials J. Nucl. Mater. 467 424–31

[8] Hodille E.A. et al 2017 Retention and release of hydrogen isotopes in tungsten plasma-facing components: the role of grain boundaries and the native oxide layer from a joint experiment-simulation integrated approach Nucl. Fusion 57 076019

[9] Delaporte-Mathurin R., Hodille E.A., Mougenot J., Charles Y. and Grisolia C. 2021 Parametric optimisation based on TDS experiments for rapid and efficient identification of hydrogen transport materials properties Nucl. Mater. Energy 27 100984

[10] Markelj S., Ogorodnikova O.V., Pelicon P., Schwarz-Selinger T. and Čadež I. 2013 Temperature dependence of D atom adsorption on polycrystalline tungsten Appl. Surf. Sci. 282 478–86

[11] Hodille E.A., Založnik A., Markelj S., Schwarz-Selinger T., Becquart C.S., Bisson R. and Grisolia C. 2017 Simulations of atomic deuterium exposure in self-damaged tungsten Nucl. Fusion 57 056002

[12] Založnik A., Markelj S., Schwarz-Selinger T. and Schmid K. 2017 Deuterium atom loading of self-damaged tungsten at different sample temperatures J. Nucl. Mater. 496 1–8

[13] Guterl J., Smirnov R.D. and Snyder P. 2019 Effects of surface processes on hydrogen outgassing from metal in desorption experiments Nucl. Fusion 59 096042

[14] Schmid K. and Zibrov M. 2021 On the use of recombination rate coefficients in hydrogen transport calculations Nucl. Fusion 61 086008

[15] Hodille E.A., Payet M., Marascu V., Peillon S., Mougenot J., Ferro Y., Delaporte-Mathurin R., Leblond F., Bernard E. and Grisolia C. 2021 Modelling tritium adsorption and desorption from tungsten dust particles with a surface kinetic model Nucl. Fusion 61 086030

[16] Pick M.A. and Sonnenberg K. 1985 A model for atomic hydrogen-metal interactions—application to recycling, recombination and permeation J. Nucl. Mater. 131 208–20

[17] Ajmalghan M., Piazza Z.A., Hodille E.A. and Ferro Y. 2019 Surface coverage dependent mechanisms for the absorption and desorption of hydrogen from the W(110) and w(100) surfaces: a density functional theory investigation Nucl. Fusion 59 106022

[18] Piazza Z.A., Ajmalghan M., Ferro Y. and Kolasinski R.D. 2018 Saturation of tungsten surfaces with hydrogen: a density functional theory study complemented by low energy ion scattering and direct recoil spectroscopy data Acta Mater. 145 388–98

[19] Ferro Y., Augustin Hodille E., Denis J., Piazza Z.A. and Ajmalghan M. 2023 Hydrogen and oxygen on tungsten (110) surface: adsorption, absorption and desorption investigated by density functional theory Nucl. Fusion 63 036017

[20] Bergstrom Z.J., Li C., Samolyuk G.D., Uberuaga B.P. and Wirth B.D. 2019 Hydrogen interactions with low-index surface orientations of tungsten J. Phys.: Condens. Matter 31 255002

[21] Dunand A., Minissale M., Faure J.-B., Gallais L., Angot T. and Bisson R. 2022 Surface oxygen versus native oxide on
tungsten: contrasting effects on deuterium retention and release *Nucl. Fusion* **62** 054002

[22] Rodríguez-Fernández A., Bonnet L., Larrégaray P. and Díez Muño R. 2023 How adsorbed oxygen atoms inhibit hydrogen dissociation on tungsten surfaces *J. Phys. Chem. Lett.* **14** 1246–52

[23] Nahm T.-U. and Gomer R. 1997 The adsorption of hydrogen on W(110) and Fe covered W(110) surfaces *Surf. Sci.* **375** 281–92

[24] Hodille E.A., Markelj S., Pecovnik M., Ajmalghan M., Piazza Z.A., Ferro Y., Schwarz-Selinger T. and Grisolia C. 2020 Kinetic model for hydrogen adsorption in tungsten with coverage dependent surface mechanisms *Nucl. Fusion* **60** 106011

[25] Montupet-Leblond F., Corso L., Payet M., Delaporte-Mathurin R., Bernard E., Charles Y., Mougenot J., Vartanian S., Hodille E.A. and Grisolia C. 2021 Permeation and trapping of hydrogen in Eurofer97 *Nucl. Mater. Energy* **29** 101062

[26] Hodille E.A., Ferro Y., Fernandez N., Becquart C.S., Angot T., Layet J.M., Bisson R. and Grisolia C. 2016 Study of hydrogen isotopes behavior in tungsten by a multi trapping macroscopic rate equation model *Phys. Scr.* **2016** 014011

[27] Pečovnik M., Hodille E.A., Schwarz-Selinger T., Grisolia C. and Markelj S. 2020 New rate equation model to describe the stabilization of displacement damage by hydrogen atoms during ion irradiation in tungsten *Nucl. Fusion* **60** 036024

[28] Ziegler J.F., Ziegler M.D. and Biersack J.P. 2010 SRIM—the stopping and range of ions in matter (2010) *Nucl. Instrum. Methods Phys. Res. B* **268** 1818–23

[29] Hindmarsh A.C. *et al* 1983 Odepack, a systemaized collection of ode solvers *Scientific Computing (IMACS Transactions on Scientific Computation vol 1)* ed R.S. Stepleman (North Holland) pp 55–64

[30] Radhakrishnan K. and Hindmarsh A.C. 1993 Description and use of Isode, the livemore solver for ordinary differential equations (available at: https://computing.lnl.gov/sites/default/files/ODEPACK_pub2_u113855.pdf)

[31] McNabb A. and Foster P.K. 1963 A new analysis of the diffusion of hydrogen in rion and ferritic steels *J. Trans. Metall. Soc.* **AIME** **227** 618–27

[32] Fernandez N., Ferro Y. and Kato D. 2015 Hydrogen diffusion and vacancies formation in tungsten: density functional theory calculations and statistical models *Acta Mater.* **94** 307–18

[33] Heinola K. and Ahlgren T. 2010 Diffusion of hydrogen in bcc tungsten studied with first principle calculations *J. Appl. Phys.* **107** 113531

[34] Vineyard G.H. 1957 Frequency factors and isotope effects in solid state rate processes *J. Phys. Chem. Solids* **3** 121–7

[35] Wang Z., Sang C. and Wang D. 2020 Modelling of deuterium retention and outgassing in self-damaged tungsten under low-energy atomic d flux irradiation: the effects of surface processes *J. Nucl. Mater.* **540** 152390

[36] Johnson D.F. and Carter E.A. 2010 Hydrogen in tungsten: absorption, diffusion, vacancy trapping and decohesion *J. Mater. Res.* **25** 315–27

[37] Piazza Z.A., Ajmalghan M., Kolasinski R.D. and Ferro Y. 2020 A density functional theory based thermodynamic model of hydrogen coverage on the W(110) surface *Phys. Scr.* **2020** 014025

[38] Piazza Z.A., Kolasinski R.D., Ajmalghan M., Hodille E.A. and Ferro Y. 2021 Predictive atomistic model for hydrogen adsorption on metal surfaces: comparison with low-energy ion beam analysis on tungsten *J. Phys. Chem. C* **125** 16086–96

[39] Whitten J.E. and Gomer R. 1998 The coadsorption of H and O on the W(110) plane *Surf. Sci.* **409** 16–26

[40] Frauenfelder R. 1969 Solution and diffusion of hydrogen in tungsten *J. Vac. Sci. Technol.* **6** 388–97

[41] Martinazzo R., Assoni S., Marinoni G. and Tantardini G.F. 2004 Hot-atom versus Eley–Rideal dynamics in hydrogen recombination on Ni(100). I. The single-adsorbate case *J. Chem. Phys.* **120** 8761–71

[42] Tamm P.W. and Schmidt L.D. 1971 Binding states of hydrogen on tungsten *J. Chem. Phys.* **54** 4775–87

[43] Matveev D., Wensing M., Ferry L., Virot F., Barrachin M., Ferro Y. and Linsmeier C. 2018 Reaction-diffusion modeling of hydrogen transport and surface effects in application to single-crystalline be *Nucl. Instrum. Methods Phys. Res. B* **430** 23–30

[44] Oda T. 2023 Steady-state tritium inventory in plasma-facing tungsten for fusion reactors: an effective calculation method and implications of calculation results *J. Nucl. Mater.* **577** 154294

[45] Kulagin V., Gasparyan Y. and Degtyarenko N. 2022 Numerical estimation of the atomic fraction during laser-induced desorption of hydrogen from tungsten and beryllium *Fusion Eng. Des.* **184** 113287

[46] Drexler A., Depover T., Verbeken K. and Ecker W. 2019 Model-based interpretation of thermal desorption spectra of Fe-C-Ti alloys *J. Alloys Compd.* **789** 647–57

[47] Lu G.-H., Zhou H.-B. and Becquart C.S. 2014 A review of modelling and simulation of hydrogen behaviour in tungsten at different scales *Nucl. Fusion* **54** 086001

[48] De Backer A., Mason D.R., Domain C., Nguyen-Manh D., Marinica M.-C., Ventelon L., Becquart C.S. and Dudarev S.L. 2017 Hydrogen accumulation around dislocation loops and edge dislocations: from atomistic to mesoscopic scales in bcc tungsten *Phys. Scr.* **2017** 014073

[49] Bakaev A., Grigorev P., Terentyev D., Bakaeva A., Zhurkin E.E. and Mastrikov Y.A. 2017 Trapping of hydrogen and helium at dislocations in tungsten: an *ab initio* study *Nucl. Fusion* **57** 126040

[50] Butz R. and Wagner H. 1977 Diffusion of oxygen on tungsten (110) *Surf. Sci.* **63** 448–59

[51] Uebing C. and Gomer R. 1997 The diffusion of oxygen on W(110) the influence of the p(2 × 1) ordering *Surf. Sci.* **381** 33–48

[52] Schmid K., Rieger V. and Manhard A. 2012 Comparison of hydrogen retention in W and W/Ta alloys *J. Nucl. Mater.* **426** 247–53

[53] Ogorodnikova O.V., Roth J. and Mayer M. 2003 Deuterium retention in tungsten in dependence of the surface conditions *J. Nucl. Mater.* **313–316** 469–77

[54] Zhou H.-B., Liu Y.-L., Jin S., Zhang Y., Luo G.-N. and Lu G.-H. 2010 Investigating behaviours of hydrogen in a tungsten grain boundary by first principles: from dissolution and diffusion to a trapping mechanism *Nucl. Fusion* **50** 025016

[55] Bakaev A., Terentyev D. and Zhurkin E.E. 2020 *Ab initio* study of the stability of h-he clusters at lattice defects in tungsten *Nucl. Instrum. Methods Phys. Res. B* **478** 269–73

[56] Heinola K., Ahlgren T., Nordlund K. and Keinonen J. 2010 Hydrogen interaction with point defects in tungsten *Phys. Rev. B* **82** 094102

[57] Gao L., Wilde M., Manhard A., von Toussaint U. and Jacob W. 2020 Hydrogen atom-ion synergy in surface lattice modification at sub-threshold energy *Acta Mater.* **201** 55–62

[58] Hodille E.A., Fernandez N., Piazza Z.A., Ajmalghan M. and Ferro Y. 2018 Hydrogen supersaturated layers in H/D plasma-loaded tungsten: a global model based on thermodynamics, kinetics and density functional theory data *Phys. Rev. Mater.* **2** 093802

[59] Nicholson E., Davis J.W. and Veer Singh C. 2023 Nanoindentation yield load reduction in tungsten caused by
tungsten: contrasting effects on deuterium retention and release *Nucl. Fusion* **62** 054002

[22] Rodríguez-Fernández A., Bonnet L., Larrégaray P. and Díez Muño R. 2023 How adsorbed oxygen atoms inhibit hydrogen dissociation on tungsten surfaces *J. Phys. Chem. Lett.* **14** 1246–52

[23] Nahm T.-U. and Gomer R. 1997 The adsorption of hydrogen on W(110) and Fe covered W(110) surfaces *Surf. Sci.* **375** 281–92

[24] Hodille E.A., Markelj S., Pecovnik M., Ajmalghan M., Piazza Z.A., Ferro Y., Schwarz-Selinger T. and Grisolia C. 2020 Kinetic model for hydrogen adsorption in tungsten with coverage dependent surface mechanisms *Nucl. Fusion* **60** 106011

[25] Montupet-Leblond F., Corso L., Payet M., Delaporte-Mathurin R., Bernard E., Charles Y., Mougenot J., Vartanian S., Hodille E.A. and Grisolia C. 2021 Permeation and trapping of hydrogen in Eurofer97 *Nucl. Mater. Energy* **29** 101062

[26] Hodille E.A., Ferro Y., Fernandez N., Becquart C.S., Angot T., Layet J.M., Bisson R. and Grisolia C. 2016 Study of hydrogen isotopes behavior in tungsten by a multi trapping macroscopic rate equation model *Phys. Scr.* **2016** 014011

[27] Pečovnik M., Hodille E.A., Schwarz-Selinger T., Grisolia C. and Markelj S. 2020 New rate equation model to describe the stabilization of displacement damage by hydrogen atoms during ion irradiation in tungsten *Nucl. Fusion* **60** 036024

[28] Ziegler J.F., Ziegler M.D. and Biersack J.P. 2010 SRIM—the stopping and range of ions in matter (2010) *Nucl. Instrum. Methods Phys. Res. B* **268** 1818–23

[29] Hindmarsh A.C. *et al* 1983 Odepack, a systemaized collection of ode solvers *Scientific Computing (IMACS Transactions on Scientific Computation vol 1)* ed R.S. Stepleman (North Holland) pp 55–64

[30] Radhakrishnan K. and Hindmarsh A.C. 1993 Description and use of Isode, the livemore solver for ordinary differential equations (available at: https://computing.lnl.gov/sites/default/files/ODEPACK_pub2_u113855.pdf)

[31] McNabb A. and Foster P.K. 1963 A new analysis of the diffusion of hydrogen in rion and ferritic steels *J. Trans. Metall. Soc.* **AIME** **227** 618–27

[32] Fernandez N., Ferro Y. and Kato D. 2015 Hydrogen diffusion and vacancies formation in tungsten: density functional theory calculations and statistical models *Acta Mater.* **94** 307–18

[33] Heinola K. and Ahlgren T. 2010 Diffusion of hydrogen in bcc tungsten studied with first principle calculations *J. Appl. Phys.* **107** 113531

[34] Vineyard G.H. 1957 Frequency factors and isotope effects in solid state rate processes *J. Phys. Chem. Solids* **3** 121–7

[35] Wang Z., Sang C. and Wang D. 2020 Modelling of deuterium retention and outgassing in self-damaged tungsten under low-energy atomic d flux irradiation: the effects of surface processes *J. Nucl. Mater.* **540** 152390

[36] Johnson D.F. and Carter E.A. 2010 Hydrogen in tungsten: absorption, diffusion, vacancy trapping and decohesion *J. Mater. Res.* **25** 315–27

[37] Piazza Z.A., Ajmalghan M., Kolasinski R.D. and Ferro Y. 2020 A density functional theory based thermodynamic model of hydrogen coverage on the W(110) surface *Phys. Scr.* **2020** 014025

[38] Piazza Z.A., Kolasinski R.D., Ajmalghan M., Hodille E.A. and Ferro Y. 2021 Predictive atomistic model for hydrogen adsorption on metal surfaces: comparison with low-energy ion beam analysis on tungsten *J. Phys. Chem. C* **125** 16086–96

[39] Whitten J.E. and Gomer R. 1998 The coadsorption of H and O on the W(110) plane *Surf. Sci.* **409** 16–26

[40] Frauenfelder R. 1969 Solution and diffusion of hydrogen in tungsten *J. Vac. Sci. Technol.* **6** 388–97

[41] Martinazzo R., Assoni S., Marinoni G. and Tantardini G.F. 2004 Hot-atom versus Eley–Rideal dynamics in hydrogen recombination on Ni(100). I. The single-adsorbate case *J. Chem. Phys.* **120** 8761–71

[42] Tamm P.W. and Schmidt L.D. 1971 Binding states of hydrogen on tungsten *J. Chem. Phys.* **54** 4775–87

[43] Matveev D., Wensing M., Ferry L., Virot F., Barrachin M., Ferro Y. and Linsmeier C. 2018 Reaction-diffusion modeling of hydrogen transport and surface effects in application to single-crystalline be *Nucl. Instrum. Methods Phys. Res. B* **430** 23–30

[44] Oda T. 2023 Steady-state tritium inventory in plasma-facing tungsten for fusion reactors: an effective calculation method and implications of calculation results *J. Nucl. Mater.* **577** 154294

[45] Kulagin V., Gasparyan Y. and Degtyarenko N. 2022 Numerical estimation of the atomic fraction during laser-induced desorption of hydrogen from tungsten and beryllium *Fusion Eng. Des.* **184** 113287

[46] Drexler A., Depover T., Verbeken K. and Ecker W. 2019 Model-based interpretation of thermal desorption spectra of Fe-C-Ti alloys *J. Alloys Compd.* **789** 647–57

[47] Lu G.-H., Zhou H.-B. and Becquart C.S. 2014 A review of modelling and simulation of hydrogen behaviour in tungsten at different scales *Nucl. Fusion* **54** 086001

[48] De Backer A., Mason D.R., Domain C., Nguyen-Manh D., Marinica M.-C., Ventelon L., Becquart C.S. and Dudarev S.L. 2017 Hydrogen accumulation around dislocation loops and edge dislocations: from atomistic to mesoscopic scales in bcc tungsten *Phys. Scr.* **2017** 014073

[49] Bakaev A., Grigorev P., Terentyev D., Bakaeva A., Zhurkin E.E. and Mastrikov Y.A. 2017 Trapping of hydrogen and helium at dislocations in tungsten: an *ab initio* study *Nucl. Fusion* **57** 126040

[50] Butz R. and Wagner H. 1977 Diffusion of oxygen on tungsten (110) *Surf. Sci.* **63** 448–59

[51] Uebing C. and Gomer R. 1997 The diffusion of oxygen on W(110) the influence of the $p(2 \times 1)$ ordering *Surf. Sci.* **381** 33–48

[52] Schmid K., Rieger V. and Manhard A. 2012 Comparison of hydrogen retention in W and W/Ta alloys *J. Nucl. Mater.* **426** 247–53

[53] Ogorodnikova O.V., Roth J. and Mayer M. 2003 Deuterium retention in tungsten in dependence of the surface conditions *J. Nucl. Mater.* **313–316** 469–77

[54] Zhou H.-B., Liu Y.-L., Jin S., Zhang Y., Luo G.-N. and Lu G.-H. 2010 Investigating behaviours of hydrogen in a tungsten grain boundary by first principles: from dissolution and diffusion to a trapping mechanism *Nucl. Fusion* **50** 025016

[55] Bakaev A., Terentyev D. and Zhurkin E.E. 2020 *Ab initio* study of the stability of h-he clusters at lattice defects in tungsten *Nucl. Instrum. Methods Phys. Res. B* **478** 269–73

[56] Heinola K., Ahlgren T., Nordlund K. and Keinonen J. 2010 Hydrogen interaction with point defects in tungsten *Phys. Rev. B* **82** 094102

[57] Gao L., Wilde M., Manhard A., von Toussaint U. and Jacob W. 2020 Hydrogen atom-ion synergy in surface lattice modification at sub-threshold energy *Acta Mater.* **201** 55–62

[58] Hodille E.A., Fernandez N., Piazza Z.A., Ajmalghan M. and Ferro Y. 2018 Hydrogen supersaturated layers in H/D plasma-loaded tungsten: a global model based on thermodynamics, kinetics and density functional theory data *Phys. Rev. Mater.* **2** 093802

[59] Nicholson E., Davis J.W. and Veer Singh C. 2023 Nanoindentation yield load reduction in tungsten caused by
hydrogen ion irradiation near the threshold displacement energy *Nucl. Instrum. Methods Phys. Res. B* **544** 165105

[60] Dunand A., Minissale M., Angot T. and Bisson R. 2023 Flux dependence of helium retention in clean W(110): experimental evidence for He self-trapping *Nucl. Mater. Energy* **34** 101324

[61] Schmid K. 2016 Diffusion-trapping modelling of hydrogen recycling in tungsten under ELM-like heat loads *Phys. Scr.* **2016** 014025

[62] Denis J., Hodille E.A., Marandet Y. and Ferro Y. 2022 Analytical model of hydrogen inventory saturation in the subsurface of the wall material and comparison to reaction-diffusion simulations *J. Nucl. Mater.* **570** 153972

[63] Pétuya R., Nosir M.A., Crespos C., Díez Muiño R. and Larrégaray P. 2015 Isotope effects in Eley–Rideal and hot-atom abstraction dynamics of hydrogen from tungsten (100) and (110) surfaces *J. Phys. Chem. C* **119** 15325–32

[64] Galparsoro O., Pétuya R., Busnengo F., Iñaki Juaristi J., Crespos C., Alducin M. and Larregaray P. 2016 Hydrogen abstraction from metal surfaces: when electron–hole pair excitations strongly affect hot-atom recombination *Phys. Chem. Chem. Phys.* **18** 31378–83

[65] Nordlund K., Björkas C., Vörtler K., Meinander A., Lasa A., Mehine M. and Krasheninnikov A.V. 2011 Mechanism of swift chemical sputtering: comparison of Be/C/W dimer bond breaking *Nucl. Instrum. Methods Phys. Res. B* **269** 1257–61

[66] Hodille E.A., Byggmästar J., Safi E. and Nordlund K. 2019 Molecular dynamics simulation of beryllium oxide irradiated by deuterium ions: sputtering and reflection *J. Phys.: Condens. Matter* **31** 185001

[67] Schmid K., Schwarz-Selinger T. and Theodorou A. 2023 Modeling intrinsic- and displacement-damage-driven retention in EUROFER *Nucl. Mater. Energy* **36** 101494

[68] Bisson R., Markelj S., Mourey O., Ghiorghiu F., Achkasov K., Layet J.-M., Roubin P., Cartry G., Grisolia C. and Angot T. 2015 Dynamic fuel retention in tokamak wall materials: an *in situ* laboratory study of deuterium release from polycrystalline tungsten at room temperature *J. Nucl. Mater.* **467** 432–8

[69] Virtanen P. et al (SciPy 1.0 Contributors) 2020 SciPy 1.0: fundamental algorithms for scientific computing in python *Nat. Methods* **17** 261–72
hydrogen ion irradiation near the threshold displacement energy *Nucl. Instrum. Methods Phys. Res. B* **544** 165105

[60] Dunand A., Minissale M., Angot T. and Bisson R. 2023 Flux dependence of helium retention in clean W(110): experimental evidence for He self-trapping *Nucl. Mater. Energy* **34** 101324

[61] Schmid K. 2016 Diffusion-trapping modelling of hydrogen recycling in tungsten under ELM-like heat loads *Phys. Scr.* **2016** 014025

[62] Denis J., Hodille E.A., Marandet Y. and Ferro Y. 2022 Analytical model of hydrogen inventory saturation in the subsurface of the wall material and comparison to reaction-diffusion simulations *J. Nucl. Mater.* **570** 153972

[63] Pétuya R., Nosir M.A., Crespos C., Díez Muiño R. and Larrégaray P. 2015 Isotope effects in Eley–Rideal and hot-atom abstraction dynamics of hydrogen from tungsten (100) and (110) surfaces *J. Phys. Chem. C* **119** 15325–32

[64] Galparsoro O., Pétuya R., Busnengo F., Iñaki Juaristi J., Crespos C., Alducin M. and Larregaray P. 2016 Hydrogen abstraction from metal surfaces: when electron–hole pair excitations strongly affect hot-atom recombination *Phys. Chem. Chem. Phys.* **18** 31378–83

[65] Nordlund K., Björkas C., Vörtler K., Meinander A., Lasa A., Mehine M. and Krasheninnikov A.V. 2011 Mechanism of swift chemical sputtering: comparison of Be/C/W dimer bond breaking *Nucl. Instrum. Methods Phys. Res. B* **269** 1257–61

[66] Hodille E.A., Byggmästar J., Safi E. and Nordlund K. 2019 Molecular dynamics simulation of beryllium oxide irradiated by deuterium ions: sputtering and reflection *J. Phys.: Condens. Matter* **31** 185001

[67] Schmid K., Schwarz-Selinger T. and Theodorou A. 2023 Modeling intrinsic- and displacement-damage-driven retention in EUROFER *Nucl. Mater. Energy* **36** 101494

[68] Bisson R., Markelj S., Mourey O., Ghiorghiu F., Achkasov K., Layet J.-M., Roubin P., Cartry G., Grisolia C. and Angot T. 2015 Dynamic fuel retention in tokamak wall materials: an *in situ* laboratory study of deuterium release from polycrystalline tungsten at room temperature *J. Nucl. Mater.* **467** 432–8

[69] Virtanen P. et al (SciPy 1.0 Contributors) 2020 SciPy 1.0: fundamental algorithms for scientific computing in python *Nat. Methods* **17** 261–72