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

Rate equation modelling is performed to simulate D2 and D2+D2+ 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 D2 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 D2+D2+ 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.


Introduction
In a fusion reactor, Plasma Facing Components (PFCs) are exposed to high heat flux (≈10 MWm −2 ) carried in part by an intense fluxes of hydrogen (H) isotopes [1,2].The most 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][8][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][11][12][13][14][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][18][19][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 2 and D + 2 + D 2 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 2 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 + 2 + D 2 exposure shows a reduction of the total retention which suggests the existence of another process which removes particles during D + 2 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 + 2 ) to reproduce the difference of retention between D 2 and D + 2 + D 2 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 2 exposure [25] and to which the sputtering process is added.

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][25][26][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 ∂Ω.
In this work, we simulate a single crystal of W with W(110) surface exposed to D + 2 ions and D 2 molecules.D + 2 ions are implanted in the bulk through a volume source S (m −3 s −1 ) characterized by an incident flux ϕ inc (m −2 s −1 ), an incident energy E inc (eV), and an incident angle α inc : where r(E inc , α inc ) (dimensionless) is the reflection of ions at the surface and f(E inc , α inc ) (m −1 ) the spatial distribution of the implanted D in the material.The implanted flux ϕ imp = (1 − r(E inc , α inc )) ϕ inc is the part of the incident flux which stays in the material, diffuses and can be trapped.Both r(E inc , α inc ) and f(E inc , α inc ) are evaluated with SRIM [28] while ϕ inc is defined by the experimental parameters of [21].D 2 molecules interact with the surface as a flux of molecules Γ D2 (D 2 m −2 s −1 ) calculated from the D 2 pressure p D2 (Pa) and the temperature of the gas T gas (K): with m D2 (kg) the mass of the molecule impinging the wall and k B = 8.6175 × 10 −5 eV K −1 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].

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 −3 ) and the trapped D in the trapping site of type i; with the concentration c t,i (m −3 ).The equations governing the spatial and time evolutions of both quantities are given here: 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/ √ 2 [34].Since DFT calculates that H sits in tetrahedral position in bcc W [32,33], n IS = 6ρ W .The trapping and detrapping pre-exponential frequencies are taken equal to 10 13 s −1 based on phonon calculations and harmonic transition state theory [32,33].

D(T)
Interstitial D diffusivity (m 2 s −1 ) Concentration of interstitial site (m −3 ) 6ρ W [32,33] ρ W W atomic concentration (m −3 ) 6.3382 × 10 28 m −3 Table 2. Expression of the different flux φ P in equations ( 5) and ( 6) which are reported in [11,24,25].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 2 gas [25].It considers two types of D atoms: the adsorbed D atoms with the concentration c surf (m −2 ) and the interstitial D atoms at the boundary ∂Ω, i.e. just beneath the surface, with the concentration c m (∂Ω) (m −3 ).The equations governing the spatial and time evolutions of both quantities are given here: with λ = 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-handside 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 −1 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 θ O : clean (θ O = 0), 0.50 ML of O (θ O = 0.5) and 0.75 ML of O (θ 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 surf (clean) = 1.416 × 10 19 m −2 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 θ D + θ O = 1.Thus, n surf (θ O ) = n surf (clean)(1 − θ O ) and the values for the three surfaces are reported in table 3.

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 φ sput in equation (5).It represents the sputtering of adsorbed D (D ads ) by incident D ions (D inc ) from the plasma which is schematically described in figure 1.
Assuming the sputtering reaction D inc + D ads is an elementary process, its kinetic is proportional to the concentration of its reactants: c surf is the concentration of D ads and the incident flux ϕ inc sets the amount of D inc reacting with the surface.For a fully covered surface, the frequency of the sputtering events is σ sput ϕ inc with σ sput (m 2 ) the sputtering cross section.It yields: 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 2 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 Other possible products may be obtained such as the recombination an adsorbed atom with a sputtered atom hoping on the surface.
behavior of D atoms and D 2 molecules in the edge plasma are different in a fusion reactor.The sputtering process is also quantified by the sputtering yield Y sput (dimensionless) which takes the following form: Y sput = σ sput n surf .In that case, the sputtering flux is given as: Thus, Y sput is the amount of sputtered D adsorbed atoms per incident D ions in the case of a saturated surface (θ D = 1).In the next, Y sput is a free parameter to reproduce the experimental results of Dunand et al [21].

2.2.3.
Coverage dependent E 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 des (θ D ) [15,24,[43][44][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 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 des (θ D ) is: In equation ( 9), E FD (θ D ) is the part already shown in [15,24] and the second part is the exponential decrease where α (dimensionless) indicates the amplitude of the drop of E des and β (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 DFT oversat , is in the range of 0.7-0.8eV for clean surface [17] and 0.6 eV for W(110) with pre-adsorbed O [19].These values are used to calculate α as:  9) and the solid lines are the optimized E des (θ 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.
The remaining parameters of E des are obtained by fitting the DFT data from [17,19].The obtained evolutions of E des (θ 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.

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 D2 and the incidence flux (if D + 2 are present) is ϕ 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 −1 .Table 4. First row: experimental parameters for the exposure to D 2 (only p D2 ), implantation (p D2 and ϕ inc ) and storage (between the exposure and the TDS) in Dunand et al [21].Second row: parameters calculated for the MHIMS simulations. [21]

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 surf sim = {E 0 , ∆E, θ D0 , δθ D , β} which set the evolution of the surface energies E des (θ 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 trap,i sim = {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 sput sim = {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 O 0.75ML sput 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 O 0.75ML sput .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 O 0.75ML sput ).

Surface parameters
Figure 3 shows the results of the optimization of the surface parameters, done with the TDS spectra after exposure to D 2 only, for (a) the clean surface, (b) the O 0.50ML surface and (c) the O 0.75ML surface.The cost function ϵ 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 0.50ML surface, the relative error with the optimized surface parameters is below 9%.For the O 0.75ML 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 (θ 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.4eV range.For the oxygen covered surfaces, the optimized desorption energies at the peak of the TDS spectra  3).α is calculated using equation ( 11) using E DFT oversat equal to 0.8 eV for the clean surface [17]

Bulk parameters
The bulk parameters are determined with the data for the O 0.75ML surface exposed to D 2 and D + 2 .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 O 0.75ML sput , we arbitrarily chose the ones for Y O 0.75ML sput = 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 O 0.75ML sput , 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 particle in the bulk, hence reducing the trapping rate: to compensate this lost of trapped particles, the trap concentration has to be increased.

Sputtering yield
To determine the sputtering yield for the clean surface and the O 0.50ML surface, the trapping parameters for each value of Y O 0.75ML sput 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 sput to match the experimental integrated retention 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.50ML 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 θ D (equation ( 8)).The desorption energies decreases with the coverage of O (figure 2) which means θ D on the O 0.50ML surface will be smaller than on the clean surface: to remove the same amount of particles with the sputtering, Y sput has to be higher.

Surface energies
During the optimization of the surface parameters, differences between E 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 des as the hydrogen surface coverage decreases during the TDS.However, the desorption energy distribution E des (θ 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.4eV range.For the O 0.50ML 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 0.75ML 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 des at low coverage up to 1.28 eV.It should be noted that with the energy calculated by DFT, the O 0.75ML 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 des rise up to about 1.3 eV for θ D < 0.2.This behavior possibly comes from the structure of O on the O 0.50ML and O 0.75ML surfaces observed in calculations [19] and experiments [21].For the O 0.50ML 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 0.75ML 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.

Nature of the identified defects
The contribution of the two traps can be identified on the TDS spectrum especially for the O 0.75ML 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 dt ≈ 0.9 eV).The high temperature shoulder from 500 K to 600 K corresponds to the detrapping from the bulk trap 2 (E dt ≈ 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.2eV could corresponds to impurities [47] or edge dislocation lines [48,49], assuming E dt = E bind + E diff , E bind the binding energy of H with the defect.
In principle, these detrapping energies could also correspond to trapping at grain boundaries [54,55] or monovacancies [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 ioninduced 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 + 2 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 −1 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 2 + D + 2 exposed W samples are not influenced by the presence of O atoms prior to D + 2 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 2 + D + 2 exposure.Previous TDS simulations [7,52,53] suggests that the concentration of intrinsic traps in the bulk is around 10 −5 − 10 −3 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.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.50ML surfaces is removed.For the O 0.75ML surface, the effect of the sputtering is much less pronounced indicating that the sputtering is not very efficient for this surface.Indeed, Y 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 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.50ML 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 θ D and thus low E 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.

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.75ML 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 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 × 10 18 Dm −2 ) as shown in table 7. Since the trapping parameters  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 Ysput(clean) = 3.453.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 2 exposure and the storage phase.In (a), the shadowed area shows the saturation of the surface (θ 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 O 0.75ML sput ) = 1.500.which increases the surface retention.However, the sputtering process does not allow high coverage/low E 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 0.50ML surface (but not shown here).For the O 0.75ML 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.
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.4.
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 ϕ 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 dominated by D 2 desorption as soon as the temperature is high enough.Also, because the sputtering is proportional to θ D and E des 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 22 − 10 24 Dm −2 s −1 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 E des ≈ 1.5 eV, the maximum fluxes desorbing from a saturated surfaces at about 500 K is φ des (500 K) ≈ 10 17 Dm −2 s −1 and φ des (1000 K) ≈ 10 24 Dm −2 s −1 .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 E des 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 −1 atoms is of the order of 10 −21 m 2 [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 −19 m 2 which is 2 order of magnitude 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 sputtering processes may have similarities with hot-atom recombination or Eley-Rideal recombination as the reaction proceeds 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 thermaldissipation issues.MD simulation of plasma facing materials irradiation by energetic D atoms suggests that the reflection and recycling of D through D 2 decreases with the incident 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 2 released via Langmuir-Hinshelwood recombination or energetic D 2 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.

Conclusions
In this paper, rate equation modelling has been used to simulate experiments presented in [21], where D 2 and D + 2 + D 2 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 calculations: the more O is adsorbed on the surface, the easier the desorption is.The bulk detrapping energies determined from the simulation of D + 2 + D 2 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 + 2 has to be considered.A similar processes has been used to simulated experimental 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 vacuum/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) suggests that in divertor conditions, the sputtering may also be a significant recycling channel.However, for divertor conditions, 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.

Figure 2 .
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 DFT des which fit the DFT data with equation (9) and the solid lines are the optimized E des (θ 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.

Figure 3 .
Figure 3.Comparison between the experimental and simulated TDS spectra after exposure to D 2 only for (a) the clean surface, (b) the O 0.50ML surface and (c) the O 0.75ML surface.The dashed line are the simulated spectra obtained with E des (θ D ) derived from the DFT calculations and the solid lines are the results of the optimization on the surface parameters (solid line in figure 2).

Figure 4 .
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 O 0.75ML sput for the O 0.75ML surface.The reference values are taken for Y O 0.75ML sput = 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.

.Figure 5 .
Figure 5. Sensitivity study of the optimized value of Y clean sput and Y O 0.50ML sput with the trapping parameters obtained for value of Y O 0.75ML sput

Figure 6
Figure6shows the comparison between the experimental and simulated TDS spectra after exposure to D 2 (a) and D 2 + D + 2 (b) for the three different surfaces.To highlight the effect of the sputtering process, simulations without the sputtering process are also shown in figure6(b).As mentioned earlier, only one parameter is changed to simulate the TDS spectra, Y sput .For the clean surface and the O 0.50ML surface, the position of the TDS peak is not reproduced quantitatively, however the evolution of the

Figure 6 .
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 O 0.75ML sput = 1.500).

Figure 8 .
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 2 + D +2 exposure and the storage phase.In (a), the shadowed area shows the saturation of the surface (θ D = 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 Ysput(clean) = 3.453.

Figure 9 .
Figure 9.For the O 0.75ML 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 2 + D +2 exposure and the storage phase.In (a), the shadowed area shows the saturation of the surface (θ 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 O 0.75ML

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×10 19 Dm −2 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 = φ des + φ sput ϕ imp + φ gas .(

Figure 10 .
Figure 10.Recycling coefficient R D for the clean surface case during exposure with D 2 + D + 2 .

Figure 11 .
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 × 10 18 Dm −2 s −1 to 1 × 10 24 Dm −2 s −1 .All the other exposure parameters (P D2 , E inc , θ inc ) are the one reported in table 4.

Figure 12 .
Figure 12.Illustration of the parametric optimization for the clean W(110) exposed to D 2 gas only.(a) Comparison between the experimental and simulated TDS spectra.(b) evolution of | fexp(T i ) − f sim (P 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 ϵ with the iteration.

Table 1 .
Parameters of the bulk model.The W atomic concentration is calculated considering a lattice constant of 3.16 Å.

Table 3 .
Parameters to express the D fluxes at the surface.

Table 6 .
[19]0.6 eV for the surface with adsorbed O[19].Optimized value of the trapping parameters and the sputtering yield of the clean surface and the O 0.50ML surface assuming Y O 0.75ML sput = 1.500.

Table 7 .
Integrated differential retention (Dm −2 ) for the three different surfaces.