E A Hodille

J Byggmästar

E Safi

K Nordlund

Molecular dynamics simulation of beryllium oxide irradiated by deuterium ions: sputtering and reflection

Keywords: Plasma-wall interactions, Molecular dynamics, Beryllium oxide, Deuterium

The sputtering and reflection properties of wurtzite beryllium oxide (BeO) subjected to deuterium (D) ions bombardment at 300 K with ion energy between 10 eV and 200 eV is studied by classical molecular dynamics. Cumulative irradiations of wurtzite BeO show a D concentration threshold above which an 'unphysical dramatic' sputtering is observed. From the cumulative irradiations, simulation cells with different D concentrations are used to run non-cumulative irradiations at different concentrations. Using a D concentration close to the experimentally determined saturation concentration (0.12 atomic fraction), the simulations are able to reproduce accurately the experimental sputtering yield of BeO materials. The processes driving the sputtering of beryllium (Be) and oxygen (O) atoms as molecules are subsequently determined. At low irradiation energy, between 10 eV and 80 eV, swift chemical sputtering is dominant and produces mostly ODz molecules. At high energy, the sputtered molecules are mostly BexOy molecules (mainly BeO dimer). Four different processes are associated to the formation of such molecules: the physical sputtering of BeO dimer, the delayed swift chemical sputtering not involving D ions and the detachmentinduced sputtering. The physical sputtering of BeO dimer can be delayed if the sputtering event implies two interactions with the incoming ion (first interaction in its way in the material, the other in its way out if it is backscattered). The detachment-induced sputtering is a characteristic feature of the 'dramatic' sputtering and is mainly observed when the concentration of D is close to the threshold leading to this sputtering regime.

Introduction

Beryllium (Be) materials have been chosen as materials for the first wall of current tokamak JET [START_REF] Philipps | Overview of the jet iter-like wall project[END_REF][START_REF] Brezinsek | Plasma-surface interaction in the be/w environment: Conclusions drawn from the jet-ilw for iter[END_REF] and the future tokamak ITER [START_REF] Loarte | the ITPA Scrape-off Layer, and Divertor Physics Topical Group. Chapter 4: Power and particle control[END_REF]4]. During their lifetime in the tokamak, Be plasma-facing materials are exposed to plasma (mostly hydrogen isotopes and mainly deuterium (D) during non-nuclear operations). Through their interactions with the plasma, the Be materials can be sputtered which limits the in-vessel lifetime [START_REF] Brezinsek | Study of physical and chemical assisted physical sputtering of beryllium in the jet iter-like wall[END_REF]. Furthermore, the sputtered material will migrate in the edge plasma and can be redeposited, leading to the growth of a deposited layer containing a high amount of fuel particles. From JET experiments, this co-deposition process is expected to be the main media for fuel retention in ITER [START_REF] Brezinsek | Fuel retention studies with the iter-like wall in jet[END_REF][START_REF] Heinola | Fuel retention in jet iter-like wall from post-mortem analysis[END_REF]. In addition to sputtering, D ions can be reflected back to the plasma. This reflection can affect the overall stability of the plasma depending on the way these atoms are reflected. The impact of these different processes on the plasma operations can be estimated via tokamakscale simulations. For instance, the migration of impurities in the edge plasma is simulated by ERO [START_REF] Kirschner | Simulation of the plasma-wall interaction in a tokamak with the monte carlo code ero-textor[END_REF] (or its upgrade ERO2.0 [START_REF] Romazanov | First ero2.0 modeling of be erosion and non-local transport in jet iter-like wall[END_REF]) and WALLDYN [START_REF] Schmid | An integrated model of impurity migration and wall composition dynamics for tokamaks[END_REF][START_REF] Schmid | Implementation of a diffusion convection surface evolution model in walldyn[END_REF][START_REF] Schmid | Quantitative modeling of fuel retention in the jet-c and jet-ilw wall configurations by walldyn and predictions for iter[END_REF] to estimate the growth and the position in the tokamak of co-deposited layers. On the other hand, particle and energy transport codes such as SOLPS [START_REF] Wiesen | The new solps-iter code package[END_REF] and SOLEDGE2D [START_REF] Bufferand | Numerical modelling for divertor design of the west device with a focus on plasmawall interactions[END_REF] are used to estimate the heat and particle fluxes deposited on the plasmafacing materials, which can then be used for divertor design. Such codes, being at the tokamak scales, cannot simulate the complex atomic-scale physic and chemistry happening at the surface of the materials. This is why, they usually rely on binary collision approximation (BCA) or molecular dynamics (MD) inputs. Despite its low computational cost compared to MD, BCA results are not fully reliable when chemistry is involved, which is the case with Be materials [START_REF] Björkas | Molecules can be sputtered also from pure metals: sputtering of beryllium hydride by fusion plasmawall interactions[END_REF][START_REF] Safi | Atomistic simulations of the effect of reactor-relevant parameters on be sputtering[END_REF][START_REF] Safi | Multiscale modelling to relate beryllium surface temperature, deuterium concentration and erosion in fusion reactor environment[END_REF]. Thus, MD can be a better method to give physically accurate input data to these large-scale codes.

Concerning MD data for the interactions between pure Be and D, different works have already been carried out [START_REF] Björkas | Molecules can be sputtered also from pure metals: sputtering of beryllium hydride by fusion plasmawall interactions[END_REF][START_REF] Safi | Atomistic simulations of the effect of reactor-relevant parameters on be sputtering[END_REF][START_REF] Safi | Multiscale modelling to relate beryllium surface temperature, deuterium concentration and erosion in fusion reactor environment[END_REF] and successfully used in ERO [START_REF] Björkas | Molecules can be sputtered also from pure metals: sputtering of beryllium hydride by fusion plasmawall interactions[END_REF][START_REF] Lasa | and JET Contributors 8. Ero modeling and sensitivity analysis of locally enhanced beryllium erosion by magnetically connected antennas[END_REF]. However, pure Be is easy to oxidize and a beryllium oxide (BeO) layer is expected to form on top of plasma-facing components. Such a layer has indeed been observed in JET samples by Raman spectroscopy [START_REF] Kumar | Identification of beo and beoxdy in melted zones of the jet be limiter tiles: Raman study using comparison with laboratory samples[END_REF] and large concentration of oxygen (O) atoms is observed in post-mortem analysis of Be co-deposited layers on the JET divertor [START_REF] Heinola | Fuel retention in jet iter-like wall from post-mortem analysis[END_REF]. Therefore, the sputtering and reflection properties of BeO irradiated by D are also needed to complete the picture of the plasma/wall interactions for Be materials. To that end, MD simulations have been carried out that simulate the irradiation of wurtzite BeO by D ions at 300 K. The present paper reports the results of these simulations, focusing on the sputtering and reflection yields. Particular care has been devoted to the study of the chemical products of sputtering/reflection events.

Method

The irradiation of Beryllium Oxide (BeO) is modelled using classical MD. This atomistic scale technique has the suitable time scale to study bond breaking and formation, that can happen when a material is irradiated with energetic incident particles. The simulations presented here have been performed using the MD code PARCAS [START_REF] Nordlund | Defect production in collision cascades in elemental semiconductors and fcc metals[END_REF]. The snapshot and trajectories of the MD simulation results presented here are done using the OVITO software [START_REF] Stukowski | Visualization and analysis of atomistic simulation data with ovitothe open visualization tool[END_REF].

Interatomic potentials

To give relevant results, any MD simulations require accurate interatomic potentials which express the interaction between each atom of the simulation cell. The potential formalism used in this study has been developed by Tersoff [START_REF] Tersoff | New empirical model for the structural properties of silicon[END_REF][START_REF] Tersoff | New empirical approach for the structure and energy of covalent systems[END_REF] and is presented in Appendix A. For this study, the Be-Be and the Be-H potential are taken from [START_REF] Björkas | Interatomic potentials for the Be-C-H system[END_REF] (version II), the H-H potential is taken from [START_REF] Brenner | Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond films[END_REF], the O-O is taken from [START_REF] Erhart | Analytic bond-order potential for atomistic simulations of zinc oxide[END_REF] and the Be-O is taken from [START_REF] Byggmästar | Analytical bond order potential for simulations of BeO 1D and 2D nanostructures and plasma-surface interactions[END_REF]. The parameters of these potentials are presented in the Appendix A (table A1). To this date, no O-H potential exist within the formalism used in this study. Indeed, the O-H interaction can be quite complex when a large number of water molecules are considered (especially through hydrogen bonds) and it usually requires supplementary reactive terms [START_REF] Ni | A reactive empirical bond order (REBO) potential for hydrocarbonoxygen interactions[END_REF]. However, since we are interested in plasma-surface interaction, we do not expect liquid water to form on the surface of plasma-facing materials. In this context, an analytical bond order potential (ABOP) has been developed within the formalism used for the other potentials. The parameters of this potential as well as the testing properties are presented in Appendix A. It has been built so that it describes accurately the O-H dimer bond calculated by density functional theory (DFT) [START_REF] Ni | A reactive empirical bond order (REBO) potential for hydrocarbonoxygen interactions[END_REF] and the behavior of H as interstitial in BeO (perfect wurtzite) calculated by DFT [START_REF] Marinopoulos | Defect levels and hyperfine constants of hydrogen in beryllium oxide from hybrid-functional calculations and muonium spectroscopy[END_REF][START_REF] Hodille | Hydrogen in beryllium oxide investigated by DFT: on the relative stability of charged-state atomic versus molecular hydrogen[END_REF]. Especially, recent DFT calculations predict the neutral H 2 molecules to be the most stable configuration of interstitial H in BeO [START_REF] Hodille | Hydrogen in beryllium oxide investigated by DFT: on the relative stability of charged-state atomic versus molecular hydrogen[END_REF] which is well reproduced by the presently used potential (see table A2 in appendix Appendix A).

Simulation set-up

The BeO cells used in this study are irradiated with deuterium (D) ions at different energies at 300 K. Two types of simulations are carried out:

-Cumulative irradiations: the final state of a bombardment event is used as the initial frame of the next irradiation event. -Non-cumulative irradiations: the initial frame of a bombardment event is always the same cell.

The cumulative irradiations allow to build up a concentration in the material as it would occur in a material that is irradiated for a long period of time. However, because of the time/space limitation of MD, the fluxes of incident ions are ∝ 10 28 m -2 s -1 which is 7-9 orders of magnitude higher than what is estimated for the ITER Be first wall [START_REF] Roth | Recent analysis of key plasma wall interactions issues for ITER[END_REF]. Thus, the concentration of D resulting from cumulative irradiations can be unrealistically high. Non-cumulative irradiations allow to get rid of the high-flux problem since each irradiation event is independent from the others. However, such irradiations on pure BeO lack the possibility for the incident D ions to interact with any D that would be previously implanted. BeO crystallizes in the wurtzite (hexagonal symmetry) structure under standard conditions. DFT calculations by Freeman et al. [START_REF] Claeyssens | Growth of ZnO thin filmsexperiment and theory[END_REF][START_REF] Freeman | Graphitic nanofilms as precursors to wurtzite films: Theory[END_REF] suggest that wurtzite films, like BeO, terminate with {0001} surface thanks to graphitic nanofilms as precursors. Furthermore, experiments and simulations suggest that the growth of an oxide layer on top of Be is limited by the diffusion of Be through the oxide layer [START_REF] Byggmästar | Analytical bond order potential for simulations of BeO 1D and 2D nanostructures and plasma-surface interactions[END_REF][START_REF] Roth | Oxidation of beryllium and exposure of beryllium oxide to deuterium plasmas in PISCES B[END_REF]. Therefore, the topmost atoms are more likely to be oxygen atoms. Thus, it has been decided to simulate the irradiation of (000 1) surface, where the topmost atoms are O.

The BeO cells used in this study are thus perfect wurtzite crystal with a 25×24 Å2 rectangular cross section in the (x,y) plane elongated in the z direction. The BeO wurtzite cells are initially relaxed at the desired temperature for 50 ps with periodic boundary conditions in x,y and z directions. The (000 1) surface is then opened in the z direction (the periodic boundary condition are kept in x and y directions). The irradiation are done at normal incident ion in the z direction with energies of 10 eV, 30 eV, 50eV, 80 eV, 140 eV and 200 eV. To avoid channeling, the incident angle is tilted by 7 • from the normal incident and the cell thickness in the z direction is different for different energies. Namely, above 80 eV, Z=104 Å and below 80 eV, Z=68 Å.

In both cumulative and non-cumulative irradiations, each single bombardment starts by creating a recoil D ion 20 Å above the surface. The irradiations are done at 300 K. Each event lasts 7 ps that is divided into two parts depending on the temperature control applied. During the first 5 ps, the temperature is controlled by using the Berendsen's method [START_REF] Berendsen | Molecular dynamics with coupling to an external bath[END_REF] at the borders and the bottom of the cell. The two bottom layers are fixed. The temperature control during this period is not applied on the entire cell because it would rescale the velocity of the atoms in the region of interest for the sputtering, which is not desired. Such temperature control is also used, in this kind of irradiation simulations, as a damping mechanism for elastic waves [START_REF] Vörtler | Molecular Dynamics Simulations of Deuterium Trapping and Re-emission in Tungsten Carbide[END_REF]. It has been shown to have no impact on the sputtering results for system in the size range considered here [START_REF] Vörtler | Molecular Dynamics Simulations of Deuterium Trapping and Re-emission in Tungsten Carbide[END_REF]. Applying the temperature control only at the boundaries leads to an increase of the temperature to 325 K after 5 ps with an ion energy of 200 eV. Accumulating many impacts would then lead to an undesired increase of the temperature up to 500 K. Thus, to ensure that the temperature at the beginning of each bombardment is 300 K, during the remaining 2 ps, the temperature of all atoms (except the fixed ones) is quenched to 300K. To simulate a uniform bombardment, for both cumulative and non-cumulative irradiations, the simulation box is randomly shifted over the periodic boundaries in the x and y directions before each impact. The position of the recoil is always the center of simulation box in the (x,y) plane. The box is shifted instead of the position of the recoil in the (x,y) plane in order to prevent the incident D from hitting the surface in a part of the material which is temperature controlled. For the cumulative irradiations, the incident fluxes are about 2.4 × 10 28 m -2 s -1 and 3000 cumulative impacts are simulated.

Simulation Results and Discussion

Sputtering of Be and O atoms

The sputtering yield of X ∈ {Be + O, Be, O} is called in the next Y X and expresses the number of sputtered atom per incident ions. As explained in [START_REF] Träskelin | Molecular dynamics simulations of hydrogen bombardment of tungsten carbide surfaces[END_REF], in cumulative irradiations, after the i th impact, Y X is estimated by dividing the total number of sputtered X atoms by the number of impacts i. The standard deviation over all individual bombardments are used to provide error bars. For the non-cumulative irradiation, Y X is estimated by averaging the number of sputtered X atoms over incoming D ions. The MD sputtering yields are compared to experimental data that have been obtained with oxidized Be samples [START_REF] Hirooka | Effect of impurities on the erosion behavior of beryllium under steady-state deuterium plasma bombardment[END_REF][START_REF] Roth | Beryllium selfsputtering: An interpolation of data for D, He, Ne and Ar[END_REF] and sintered solid BeO [START_REF] Roth | Sputtering of Be and BeO by light ions[END_REF] at room temperature. 3.1.1. Unphysical dramatic sputtering in cumulative irradiations. Figure 1 shows snapshots of the BeO cell with the (000 1) bombarded at 300 K by 200 eV D ions after different numbers of cumulative impacts. Up to 750 impacts, the simulation cell stays compact even though the cell expands toward the vacuum, and the implanted zone is amorphized. However, as the concentration of deuterium grows beneath the surface, the implanted layer is less and less stable leading to detachment of the top of the cell. Similar detachment, called blistering, has also been observed in MD cumulative irradiation of WC [START_REF] Vörtler | Molecular Dynamics Simulations of Deuterium Trapping and Re-emission in Tungsten Carbide[END_REF] or Fe 3 C [START_REF] Safi | Atomistic simulations of deuterium irradiation on iron-based alloys in future fusion reactors[END_REF]. We note, however, that this process of a layer being removed after hydrogen irradiation from a smooth surface has also been observed experimentally for higher ion energies [START_REF] Bister | Ranges of the 0.32 mev h+ and 0.72 mev h+ 2 ions in si and ge[END_REF][START_REF] Rauhala | Exfoliation of gaas caused by mev1h and 4he ion implantation at 100, 110 axial and random orientations[END_REF], and is the basis of the industrial SMARTCUT(TM) process [START_REF] Bruel | Process for the production of thin semiconductor material films[END_REF]. In the WC and Fe 3 C cases, the blistering is initiated by the accumulation of D that eventually form a bubble. According to Vörtler and Nordlund [START_REF] Vörtler | Molecular Dynamics Simulations of Deuterium Trapping and Re-emission in Tungsten Carbide[END_REF], the detachment of the upper part is linked to the D 2 gas pressure. In figure 1, no clear bubble is formed (with an accumulation of D 2 like in [START_REF] Vörtler | Molecular Dynamics Simulations of Deuterium Trapping and Re-emission in Tungsten Carbide[END_REF]) but still the detachment of the upper part of the simulation cell is observed. This detached part disintegrate eventually leading to an "unphysical dramatic" sputtering. The detachment shown in figure 1 is an extreme case, since the layer detachment requires the implantation depth to be deep enough (so the ion energy to be high enough). Such detachment has been observed for 140 eV and 80 eV as well. Despite showing no detachment, 10 eV and 30 eV cumulative irradiations with 3000 ions also lead to an "unphysical dramatic" sputtering. The sputtering yields in that regime are reported in figure 2. These values are orders of magnitude higher than the experimental ones. In addition, above 80 eV, the sputtering yield is close to 1 (even higher than one for 200 eV) which is an artifact of the detachment described in figure 1. It has to be noted that the dramatic unphysical sputtering come together with and dramatic unphysical reflection (or re-emission) of D in several forms (atoms, molecules) as it will be described in section 3.2. This dramatic sputtering is clearly an effect caused by the extreme flux limited by the MD scales.

To investigate the conditions leading to such [START_REF] Roth | Beryllium selfsputtering: An interpolation of data for D, He, Ne and Ar[END_REF] high sputtering, the profiles of Be, O and D atoms are calculated. For the 200 eV irradiation, these profiles are given on figure 3 after 350, 750 and 1150 impacts. After 350 impacts, the concentration of Be/O in the implantation zone (between 0 and 50 Å) stays equivalent to the concentration of Be/O atoms in the unimplanted part. After 750 ions, The amount of D has increased compared to 350 impacts in the implantation zone. Calling N X∈{Be,O,D} the amount of atom of species X, one can define the concentration of D, in atomic fraction (at.fr.), as:

c D = N D N Be + N O + N D . (1)
A significant drop of the concentration of Be/O atoms is observed where c D is the highest. The cause of this drop is the implantation of D atoms which stretches the material, which is clearly seen by the motion of the position of the surface. The part close to the surface, where few D atoms have stopped, experienced less stretching and thus a smaller decrease of the concentration of Be/O. After 1150 ions, the amount of D in the implantation zone is about the same as after 750 impacts, despite the 400 D ions that might have been added to the material. Thus, to understand the condition under which the dramatic sputtering occurs, one can look at the evolution of the maximum concentration of D with the number of impacts (figure 4(a)).

For any energies, the maximum concentration of D exhibits two trends. It first increases linearly and then saturates around 0.3 at.fr. which represents the situation where there is an equivalent amount of Be, O and D atoms (N D ≈ N Be ≈ N O ). For 10 eV, the saturation concentration is higher than 0.3 and between 0.4 and 0.5. For that energy, the ions are essentially adsorbed on the surface and some are implanted below the first atomic layer. In that case, the D form O-D at the surface and the saturation occurs when N D ≈ N O (since the topmost atoms are only O for the (000 1) surface).

Once the saturation level is reached, to sustain the high flux (imposed by the MD limitation on time/space scale), each impact will tend to remove previously implanted D atoms. At high energy, the D atoms are displaced either closer to the surface or deeper in the bulk but cannot reach the surface immediately. The stretching and subsequent detachment described in both figure 1 and figure 3 follows. At low energy, the implanted atoms are closer to the surface: the gap between the peak of hydrogen concentration and the surface is not large enough to have a few-Å-thick part of cell fly away. Thus, the replacement is not followed by a big detachment as in figure 1. However, since the D atoms can be bounded to O, Be or D inside BeO, the replaced D take away some atoms from the materials leading to a large number of sputtering events: this will be further explained as the detachment-induced sputtering. This is the starting point of the unphysical dramatic sputtering. Thus, unlike WC and Fe 3 C cases, the detachment is not induced by a high pressure in a subsurface bubble but by a increase of the D concentration below the surface. The high formation energy of interstitial H in BeO suggests (see table A2) weak interactions between D and the atoms of BeO. Thus, the accumulation of D below the surface most likely induces amorphization by weakening the overall binding energy in that region. This amorphization stretches the material and eventually lead to the detachment and dramatic sputtering. The difference with WC and Fe 3 C might be because bubble or blister are harder to form in BeO. Indeed, experiments by Roth et al. [START_REF] Roth | Deuterium retention and out-gassing from beryllium oxide on beryllium[END_REF] do not report any blistering or striking surface modification of oxidized Be samples while Be samples do exhibit blisters. Vörtler et al. [START_REF] Vörtler | Molecular Dynamics Simulations of Deuterium Trapping and Re-emission in Tungsten Carbide[END_REF] stopped the simulations when the blistering occurred. Thus, in order to get rid of the unphysical dramatic sputtering, a concentration threshold has to be defined so that only the bombardment occurring before reaching this concentration threshold are considered. As the concentration of D is fluctuating in the saturation part (figure 4(a)), the condition that limits the number of considered impacts is 0.25 ≤ c max D ≤ 0.3 atomic fraction (at.fr.). The amount of impacts needed to reach this condition is reported in figure 4(b) as a function of the ion energy. It increases with the energy being 350 atoms at 10 eV and 750 at 200 eV. Indeed, as the energy increases, the implantation range increase as well as the spreading of the D atoms around this range. It means that at a given depth, it requires more impacts to reach a given concentration for a higher energy. This is clearly shown on figure 4 Compared to the data including all 3000 impacts, the "limited" cumulative data are much closer to the experimental data, especially the high energy ones.

However, the low energy data (below 80 eV) still remain about 2-5 times higher than the experimental one. In addition, the error bars are quite large (especially considering the logarithmic scale). This poor estimation of the sputtering yield may be explained by the very small amount of considered atoms. Indeed, considering 350 impacts (at 10 eV) leads to only 1 sputtering event which is obviously not enough to have a good estimate of the sputtering yield. To increase the accuracy of the estimate of the sputtering yield, it has been decided to run non-cumulative irradiations. To have an understanding of the impact of the D concentration on the sputtering yield, different simulations are considered with different D concentrations as proposed in [START_REF] Salonen | Swift chemical sputtering of amorphous hydrogenated carbon[END_REF]. We considered first the situation where there is no D atom in the BeO (c D = 0 at.fr.). Then, another D concentration considered is the saturation concentration obtained with the MD simulation on figure 4(a). For this concentration, we take the cell resulting from the cumulative irradiation after the amount of impacts given by figure 4(b). Experiments [START_REF] Roth | Oxidation of beryllium and exposure of beryllium oxide to deuterium plasmas in PISCES B[END_REF][START_REF] Roth | Deuterium retention and out-gassing from beryllium oxide on beryllium[END_REF][START_REF] Behrisch | Implantation of 5kev deuterium in beo[END_REF][START_REF] Mayer | Codeposition of hydrogen with beryllium, carbon and tungsten[END_REF] suggest that the saturation concentration of D in BeO at room temperature is about 0.12 at.fr.. This saturation concentration is different to the one obtained in the MD simulations because in the MD cumulative irradiations, the extreme flux leads to an extreme saturation that cannot be reached under experimental conditions: in the simulations, the diffusion that could spread the D in the bulk decreasing the D concentration is basically non-visible due to the time scale considered. Thus, we also considered the experimental saturation concentration of 0.12 at.fr.. The number of cumulated impacts required to reach such concentration is reported on figure 4(b). As for the saturation concentration, we take the cell resulting from the cumulative irradiation limited to this number of impacts as initial cell for the noncumulative irradiations.

Effect of the D concentration on the sputtering yield.

Figure 5 shows the comparison between the sputtering obtained by the MD simulations and the experimental data. First, one can clearly see that the sputtering yields increase with the D concentration: there is a bit more than one order of magnitude between the two extrema. One can also note that these two extrema underestimate (c D = 0 at.fr.) or overestimate (0.25 ≤ c D ≤ 0.3 at.fr.) the experimental sputtering yield. On the other hand, the sputtering yields obtained with the experimental saturation limit (c D ≈ 0.12 at.fr.) reproduce well the experimental sputtering yields. One can also see that the sputtering yields obtained with the limited cumulative irradiations range between the c-limited nc cD=0 at.fr. nc cD≈0.12 at.fr. nc 0.25≤cD≤0.3 at.fr. sputtering yield of the non-cumulative irradiations. It makes sense as in the non-cumulative irradiations, the concentrations of D are frozen to a given value, while the concentrations of D increase progressively in the limited cumulative irradiations. However, as already pointed out, the error bars are too large for limited cumulative irradiations due to a small number of impacts. They are much smaller for noncumulative irradiations as the number of impacts has been increased. At 10 eV, having c D > 0 impacts greatly the sputtering yield. Indeed, for the BeO cell without D, no sputtering events are reported at this energy while there are for BeO cell containing D. Only O atoms are sputtered at 10 eV and a condition for them to be sputtered is that they have to be bounded to D atoms on the surface. Thus, their potential energy is lower than if they would be on a perfect (000 1), which require a smaller kinetic energy to sputter them. A closer investigation of these sputtering events will be made in section 3.1.3, when the chemical sputtering will be investigated.

From the low (or non-existent) sputtering yield at 10 eV, one can estimate that the threshold for sputtering of BeO is between 30 eV and 10 eV. Such threshold agrees with the one determined experimentally by Roth et al. [START_REF] Roth | Sputtering of Be and BeO by light ions[END_REF] which is 22 eV for D on BeO. One can note that, at 10 eV, a sputtered O atoms would be easily replaced even under a very small partial pressure of O 2 , making the sputtering yield at 10 eV difficult to determine experimentally.

Figure 6 reports the evolution with the ion energy of the ratio Y O /Y Be+O for the different D concentrations in the non-cumulative irradiations. This ratio gives the average proportion of O atoms sputtered per sputtering events. If it is above 0.5, O atoms are preferentially sputtered while if it is below 0.5, Be are preferentially sputtered. For all energies (above 10 eV) and concentrations, the ratio Y Be /Y O spans around 0.5 between 0.3 and 0.7. Roth et al. [START_REF] Roth | Oxidation of beryllium and exposure of beryllium oxide to deuterium plasmas in PISCES B[END_REF] reports an equal sputtering of Be and O at room temperature while irradiating a BeO layer by 30-50 eV D ions. In figure 6, if one considers c D ≈ 0.12 at.fr., the simulation data are consistent with the observation of Roth et al. Indeed, even though the sputtering at 30 eV seems to be largely dominated by the sputtering of O, at 50 eV, Be and O are sputtered rather equivalently. Considering that figure 5 reports a sputtering yield about 5 times higher at 50 eV than at 30 eV, one could expect the sputtering yield of 50 eV to dominate when irradiating BeO with 30 -50 eV. Roth et al. reported a sputtering yield of 3×10 -3 in their experiment. Doing the average of 30 eV and 50 eV, one obtains 4.1×10 -3 which agrees quite well considering we only took the extrema of the energy used by Roth et al.

Looking in more detail at the evolution of the ratio with the different energies one can distinguish different features for each concentration. First, at 10 eV, as already pointed out (Y O /Y Be+O = 1), only O atoms are sputtered for non-zero D concentrations. This is in agreement with pure Be sputtering studies [START_REF] Björkas | Molecules can be sputtered also from pure metals: sputtering of beryllium hydride by fusion plasmawall interactions[END_REF][START_REF] Safi | Atomistic simulations of the effect of reactor-relevant parameters on be sputtering[END_REF][START_REF] Safi | Multiscale modelling to relate beryllium surface temperature, deuterium concentration and erosion in fusion reactor environment[END_REF] that report negligible sputtering at 10-20 eV. In addition, the topmost atoms of the (000 1) surface are O. Thus, 10 eV ions that only interact with the atoms close to the surface mainly interact with O atoms. On the other hand, because of the chemical affinity between O and H, the incident D ions can easily be bound to O and remove it from the subsurface, forming O-D dimer or O 2 D molecules. This chemical sputtering will be investigated further in section 3.1.3. At 30 eV, even though some Be atoms might be sputtered, the sputtering is still dominated by the chemical sputtering of O atoms. For higher energy, one can see the effect of the concentration on the sputtering regime. For 0.25 ≤ c D ≤ 0.3 at.fr., the sputtering is dominated by the sputtering of O atoms especially at low energy (above 140 eV, both O and Be sputtering yields equilibrate). In contrast, for c D = 0 at.fr., the sputtering seems to be dominated by Be sputtering once the energy is above 50 eV. For the intermediate concentration determined as the experimental D saturation, the sputtering regime is between these two extrema, with a evolution similar to the case with a c D = 0 at.fr. (with a dip in the O sputtered fraction around 50-80 eV). Therefore, the sputtering of O seems to be positively correlated with the D concentration inside the material.

Y O /Y
Be+O nc cD=0 at.fr. nc cD≈0.12 at.fr. nc 0.25≤cD≤0.3 at.fr.

Chemical analysis of the sputtering yield.

To investigate further the dependence of the sputtering yield of O with the D concentration and to understand the trends reported in figure 6, the sputtering events leading to the formation of chemically bounded clusters are reported and analyzed. From this analysis, the sputtering yield of molecules Y chem Be+O is calculated as the average of the number of sputtered Be or O atoms in molecules (it is the same way of calculating the total sputtering yields Y Be+O presented in section 3.1). Figure 7 reports the evolution with the ion energy of the fraction of this sputtering relatively to the total sputtering Y chem Be+O /Y Be+O . As for the sputtering of O atoms, the fraction of sputtered atoms inside molecules is positively correlated with the concentration of D in the simulation cell, especially between 80 eV and 200 eV. Another trend is seen in figure 7 showing a negative correlation between the fraction of sputtered atoms inside molecules and the ion energy, even though there is a small dip at 80 eV c D ≈ 0.12 at.fr.. Such negative correlation is also reported for pure Be sputtering [START_REF] Björkas | Molecules can be sputtered also from pure metals: sputtering of beryllium hydride by fusion plasmawall interactions[END_REF][START_REF] Safi | Multiscale modelling to relate beryllium surface temperature, deuterium concentration and erosion in fusion reactor environment[END_REF]. One could assume that the the mechanism driving the production of these molecules is the swift chemical sputtering (SCS) described in [START_REF] Björkas | Molecules can be sputtered also from pure metals: sputtering of beryllium hydride by fusion plasmawall interactions[END_REF][START_REF] Salonen | Swift chemical sputtering of amorphous hydrogenated carbon[END_REF][START_REF] Salonen | Bond-breaking mechanism of sputtering[END_REF][START_REF] Nordlund | Mechanism of swift chemical sputtering: Comparison of be/c/w dimer bond breaking[END_REF]. The efficiency of such a mechanism is given by a balance between the energy transfer to break a bond in the material and the effective average

Y chem Be+O /Y
Be+O nc cD=0 at.fr. nc cD≈0.12 at.fr. nc 0.25≤cD≤0.3 at.fr. time the incident atom spends in the vicinity of the broken bond. Thus, one could indeed expect that, once the threshold for bond breaking is exceeded, the SCS efficiency would decrease with energy as the effective average time decreases with increasing energy.

In more detail, this mechanism has energy thresholds determined as the minimum/maximum incident energy leading to a broken bond. In [START_REF] Björkas | Molecules can be sputtered also from pure metals: sputtering of beryllium hydride by fusion plasmawall interactions[END_REF][START_REF] Nordlund | Mechanism of swift chemical sputtering: Comparison of be/c/w dimer bond breaking[END_REF], these thresholds are determined for C and Be materials doing MD simulations of dimer irradiations by D ions at various energies and recording at which energies bond breaking occurs. In these dimer irradiation simulations, the ions are directed toward the midpoint of the dimer. Doing similar irradiation with the BeO dimer, one can find that the bond breaking occurs for D energy between 12 eV and 74 eV which give the range in which the SCS should be active. Since the lower limit threshold for bond breaking is 12 eV for the dimer, one could expect no sputtering in irradiation of cell with c D = 0 at.fr. On the other hand, since the upper limit for bond breaking is around 80 eV, one could expect a drop of the fraction of sputtered atoms inside molecules around this energy. This is indeed observed in figure 7 and especially by the dip at 80 eV for c D ≈ 0.12 at.fr.. However, after that initial drop at 80 eV, the fraction of sputtered atoms inside molecules stays constant (c D = 0 at.fr. and 0.25 ≤ c D ≤ 0.3 at.fr.) or increases (c D ≈ 0.12 at.fr.) and at a fairly high level (between 30 an 70 %) which is not expected by SCS. To analyze further this effect, the sputtered molecules are sorted depending on the nature of the species inside them. Five groups can be defined: OD z , Be x O y D z , Be x H z , O 2 and Be x O y .

The first three groups define the molecules containing deuterium which could thus be produced during SCS events. The two last groups define the molecules without deuterium. The indexes x, y and z are mostly one especially when c D = 0 at.fr.. However, x and y might be bigger than one in the Be x O y group, especially when c D > 0 at.fr.. The Be 2 group is not represented since Be 2 dimer has a very low cohesive energy compared to any other Be compound [START_REF] Björkas | Interatomic potentials for the Be-C-H system[END_REF]. It follows that no Be 2 have been observed as a result of sputtering events.

Figure 8 gives the evolution with the ion energy of the fraction of these different groups in the total sputtering yields for the different D concentrations. Adding the contribution of all groups for each concentration would lead to the data reported in figure 7. For all concentrations, figure 8 shows the same trends.

Below 80 eV, the sputtered molecules are mostly in the group of D-containing molecules and mostly in the OD z group. These energies correspond to the range in which the SCS is turned on, which explains the high presence of D-containing molecules. To illustrate the creation of these OD z molecules via the SCS mechanism, trajectories of the incident and sputtered atoms are reported in figure 9 for c D = 0 at.fr. at 30 eV (a) and for c D ≈ 0.12 at.fr. at 10 eV (b).

In figure 9 (a), with an incident energy of 30 eV, the incident D atom transfers an initial kick of kinetic energy to the O atom. The D ion is then bounced back toward the surface by the layer of atoms below the surface and transfers a new kick of kinetic energy to the sputtered O atom, which finally leads to the bond breaking with the surrounding Be atoms ending with a OD dimer.

For c D ≈ 0.12 at.fr. and 10 eV (figure 9(b)), the sputtered atoms are the atoms which are already bonded to some D atoms. In the case presented in figure 9, the sputtered O atom is bound to 2 D atoms initially. Due to these bonds, its potential energy before the sputtering events is -3.43 eV. An O atom on a perfect (000 1) has a potential energy of about -5.30 eV which prevents sputtering by 10 eV D ions. These potential energies are calculated by PARCAS as j V Oj with j ranging over all the neighbors atoms of the considered O atom within the cutoff defined in table A1. As already explained, the threshold for the SCS mechanism for perfect BeO is 12 eV. It means that the energy transfer from 10 eV D ions to O atoms on the perfect (000 1) surface does not exceed the potential energy (-5.30 eV) that binds this atom to the surface. However, as the potential energy of the O decreases (in absolute value), it became easier to sputter it: this is why an O atom with a potential energy of only -3.43 eV can be sputtered by 10 eV D ions. Thus, with an incident energy of 10 eV, only the cells with c D > 0 lead to sputtering events, thanks to the decrease (in absolute value) of the potential energy of the atoms In all cases, one can mostly see OD z molecules and almost no Be x D z which could have two explanations. First, the irradiated surface is (000 1) with O atoms as the topmost atom. It means that these atoms are most likely to be sputtered by the incoming ions. However, for c D = 0 at.fr. at 30 eV, some Be atoms are sputtered, but without being bound while all O sputtered atoms are bounded in a OD dimer. One can also suggest that the preference for OD z molecules is due to the greater chemical affinity between O and D: within the used potentials, the BeH dimer have a cohesive energy of -1.30 eV/atom [START_REF] Björkas | Molecules can be sputtered also from pure metals: sputtering of beryllium hydride by fusion plasmawall interactions[END_REF] while it is -2.272 eV/atom for the OH dimer (D 0 /2 for the O-H part of table A1). Thus, while entering the vicinity of a broken Be-O bond, the D incoming atom would rather bond to the O atom. It has to be noticed that, while irradiating a BeO dimer with D ions directed toward the midpoint of the dimer with a kinetic energy of 30 eV, it produces a single Be atoms and a OD dimer. The fact that the SCS mechanism is not active above 80 eV explains why the fraction of O (figure 6), the fraction of sputtered molecules (figure 7) and the fraction of sputtered OD z molecules (figure 8) drop between 30 eV and 80 eV for c D = 0 at.fr. and c D ≈ 0.12 at.fr.. For ion energies higher than 80 eV, the nature of the sputtered molecules changes completely since they are now dominated by Be x O y molecules (which represent above 90% of all sputtered molecules for these energies). As reported by Vörtler et al., physical sputtering can lead to the formation of sputtered dimer in WC [START_REF] Vörtler | The effect of plasma impurities on the sputtering of tungsten carbide[END_REF]. In these simulations, the mechanism is very similar and starts by the physical sputtering of a Be or O atom. From there, different processes can occur (figure 10):

(i) BeO physical sputtering (figure 10 (a)): the physically sputtered Be/O atom drags with it one of its Be/O neighbors. The bond between the first neighbors of this dragged Be/O atoms are thus broken by the energy transfered from the D incident ion to the BeO dimer. One could then say that the BeO dimer is physically sputtered.

(ii) delayed BeO physical sputtering (figure 10 (b)): the physically sputtered Be/O atom starts to be removed from the initial surface. The incident D atom, through all its collision inside the BeO lattice eventually comes back toward the surface from the bulk and sputters physically another O/Be atoms close to the initial one, which leads to the late formation of a BeO dimer. Since it requires specific conditions, this event is rare.

(iii) delayed SCS sputtering (figure 10 (c)) the physically sputtered Be/O atom is directed toward another part of the surface and removes another O/Be atom via a SCS process (not involving D atoms). In figure 10 (c): the sputtered Be atom is initially at 4.07 Å away from the sputtered O atom with which it will form a BeO dimer. As can be see in figure 10 (c), once the atom is bound to the O atom, vibrational features can be seen in its trajectories that were not present in the first stage of the sputtering events. This process is enhanced when some roughness is present on the surface, since this improves the probability that the sputtered atom interacts again with the surface. It is also enhanced with a disordered surface as the BeO bond energy in such surface are lower than in a perfect surface. It is through that process that O 2 molecules can be formed and released from the surface.

(iv) Detachment-induced sputtering (figure 10 (d)): the physically sputtered Be/O atom was initially part of a cluster of Be x O y losely bonded to the rest of the material by the sputtered atom.

In figure 10 (d), the trajectory of this atom is reported in green and shows that after its interaction with the incident ions on the surface, it is pushed in the bulk. Its removal from the cluster will lead to the detachment of the cluster from the material.

It is worth noting that both delayed SCS sputtering (iii) and detachment-induced sputtering (iv) only occur for c D > 0 because the cumulative irradiation used to build the D concentrations have damaged the irradiated surface. In addition, the time elapsed between the interaction of the incident ion with the material and the release of the BeO dimer is much longer for the delayed SCS sputtering (iii) and detachment-induced sputtering (iv). In the simulations with c D = 0 at.fr., the delayed SCS sputtering (iii) and the detachment-induced sputtering (iv) are not observed since the surface is a perfect (000 1) surface. This explains why the fraction of sputtered atoms in molecules is much smaller for this concentration than for c D > 0 for high energies. Furthermore, this explains why the threshold for the creation of Be x O y molecules is 80 eV for c D = 0 at.fr. while it is 30 eV for c D > 0 (only considering the ion energies used in the simulations). Indeed, if the surface is sufficiently disordered, the detachment-induced sputtering (iv) can require only one bond-breaking event. It can thus be activated by only one interaction between the incident D ion and an atom of the target (which happen for an ion energy above 12 eV as mentioned above). On the other hand, five bond-breaking events are required to sputter physically a BeO dimer ((i) and (ii)) from a perfect surface. Indeed, on the (000 1) surface, each surface O atom has 3 Be first neighbors in the first Be layer under the surface, while each Be atom of this layer has 4 O first neighbors. To create a BeO dimer in this configuration, one has first to break two bonds between O and two of its Be first neighbors. Then, the O atom, which has gained kinetic energy from its interaction with the incident ion, has to transfer part of this kinetic energy to its last remaining Be first neighbor. This transfered kinetic energy has to be sufficient to break the 3 bonds between this Be atoms and its other O first neighbor atoms. To trigger such a process, the initial energy of the incident D ion has obviously to be higher than for breaking a single BeO bond.

One can also notice that the Be x O y D z are mainly observed for the highest value of c D with a particularly high efficiency at 30 eV (figure 8 (c)). Despite the presence of D in these molecules, they have not been created by the SCS mechanism, but via the detachment-induced sputtering (iv) described previously, in which the detached cluster contains D atoms. The formation of these molecules through process (iv) depends on the amount of D on the surface. For high energies, the D profile is spread through a wide depth range and the concentration of D on the surface is very low compared to the maximum D concentration (figure 3). Thus, the probability to form a Be x O y D z molecule rather than a Be x O y one is quite low. Inversely, for low energy, the D profile is less spread and is essentially close to the first few atomic layers below the surface. Thus, the surface D concentration is much higher than for high energy, which increases the probability of creating Be x O y D z . It follows that the amount of produced Be x O y D z is much higher for low energy than for high energy, as shown in figure 8(c). Much less Be x O y D z is observed (only at 50 eV) for c D ≈ 0.12 at.fr.. Indeed, the cells with such concentrations have been created in cumulative irradiations that experienced less cumulative impacts. It is thus less disordered than the almost saturated BeO cell which is close to the "unphysical dramatic" sputtering limit. The detachment-induced sputtering (iv) is thus less efficient. One can notice that from 50 eV to 200 eV, summing the fraction of Be x O y D z and Be x O z y gives a about constant amount of sputtered molecules.

Finally, one can note that detachment-induced sputtering (iv) is actually a part of the processes that lead to the "unphysical dramatic" sputtering. Its particular efficiency for the cells with 0.25 ≤ c D ≤ 0.3 at.fr. explains why the total sputtering yields calculated by MD for this concentration overestimate the experimental data, especially in the energy range of 30 eV -140 eV (figure 5).

Thus, once the concentration of 0.25 ≤ c D ≤ 0.3 at.fr. is reached, the dramatic sputtering regime driven by the detachmentinduced sputtering is activated even in non-cumulative irradiations.

All these chemical processes described here can help understand the behavior of the fraction of sputtered O in the total sputtering yield presented in figure 6. The high fraction of sputtered O at low energy is related to the SCS sputtering. As Be is lighter than O, it is easier to physically sputter it. Around 50-80 eV, the physical sputtering of Be is much efficient than the SCS sputtering. Finally, both O and Be equilibrates at high energy as the physical sputtering of O and Be as well as the formation of Be x O y molecules are activated.

Table 1 summarizes for which energy and for which D concentration each sputtering process is observed in the simulations. The physical sputtering of Be or O is ticked if at least one Be or O atom has been sputtered without being bound to other atoms. The total fractions of physical sputtering (Be+O) are the complementary values of the fraction Y chem Be+O /Y Be+0 reported in figure 7. The fractions of the chemical sputtering processes can be found in figure 8. It can be noted that for c D = 0 at.fr., the range of the SCS mechanism is 30 eV -80 eV as expected from the dimer irradiation. One can also note the different thresholds for physical sputtering of Be and O and the thresholds for the formation of Be x O y molecules.

Reflection of D ions and Energy

For the non-cumulative irradiations, the reflection coefficient of deuterium atoms R D is calculated as the average number of D released from the material per impact. For the cumulative irradiation, it is calculated as the number of released D from the material divided by the number of cumulative impacts. In non-cumulative simulations with c D = 0 at.fr., R D is effectively a reflection coefficient, since it reports the number of incident ions that are reflected back from the surface. In non-cumulative simulations with c D > 0 and in cumulative simulations, this is not strictly a reflection coefficient but a recycling coefficient since previously implanted D can be sputtered from the material.

Similarly to the reflection coefficient of deuterium atoms, the reflection coefficient of energy is calculated as the average of X E X k + E X p /E inc per incident atom with X being the Be, O or D sputtered/reflected atoms, E X k the kinetic energy of atom X at the end of the simulation and E X p its potential energy. The potential energy is added in this formula since sputtered molecules are released with both kinetic energy and potential energy (which is negative). If such molecules go into the plasma, the bonds might break, which will cool down the plasma explaining the negative contribution of this energy. If a large amount of very big molecules are sputtered, as during the detachment presented in figure 1,X E X k + E X p can be negative. It means that the total energy required to break all bonds will be higher than the incident energy on the material. Of course, such events are unrealistic Table 1: Summary of the different processes reported in this work for the different energies and concentration of D: • -c D = 0 at.fr., -c D ≈ 0.12 at.fr., -0.25 ≤ c D ≤ 0.3 at.fr.

Physical sputtering

Chemical sputtering Both reflection coefficients are reported on figure 11 for the D atoms (a) and for the energy (b). The results for cumulative and non-cumulative irradiations are given except the reflection coefficient of energy for the cumulative irradiation with 3000 impacts (for the reason explained above).

O Be SCS OD z Be x O y physical delayed Be x O y physical delayed SCS Be x O y detachment- induced Be x O y (D z ) 10 eV --- --- - --- --- --- --- 30 eV --- • - • - - - - 50 eV - • • • • - - 80 eV • • • • • - - 140 eV • • --- • • - - 200 eV • • --- • • - - 0
The reflection coefficient of D atoms in cumulative irradiation where all the 3000 impacts are considered is close to 1 for all energies. Such large reflection coefficient is an artifact of the dramatic sputtering: while the Be and O atoms are sputtered or detached, they carry with them a lot of previously implanted D atoms. As mentioned in the caption, the error bars for this data set are directed toward the bottom since as the impacts cumulate, the reflection coefficient increases. The values obtained after cumulating 3000 impacts give thus an upper limit. Given the large error bar, it should not be trusted to represent what could happen during an ion exposure with experimentally relevant condition (lower flux).

As for the sputtering yields, relevant values of the reflection coefficient can be obtained by limiting the number of considered impact (using the procedure described in section 3.1.1) or by running non-cumulative irradiations. One can see that for the limited cumulative irradiation and for the noncumulative irradiation with different D concentration, the same trend is observed. The reflection coefficient is close to 1, between 0.85 and 0.65, at low energy (10 eV) and it progressively decreases as the energy of the incident ions increases. Between 10 eV and 200 eV, the reflection coefficient of D atoms drop by about one order of magnitude. Such a trend from tens to hundreds of eV is also reported for instance for Fe and Fe 3 C [START_REF] Safi | Atomistic simulations of deuterium irradiation on iron-based alloys in future fusion reactors[END_REF] or for W [START_REF] Ogorodnikova | Interaction of atomic and low-energy deuterium with tungsten pre-irradiated with self-ions[END_REF]. For that energy range, this trend is due to increased kinetic energy that eases the penetration of D ions in the materials. The reflection coefficient of energy follow exactly the same trend with a drop of one order of magnitude between 10 eV and 200 eV. However, it is lower than the reflection coefficient of D atoms since part of the incident energy have been deposited in the materials (eventually creating defects or sputtering events).

Finally, one can note that the reflection coefficient of D ions seems to be positively correlated with the concentration of D in the cell. Indeed, as the concentration of D increases, the probability for the incident D ions to interact with another D, forming D 2 molecules, increases which enhances the re-emission of previously implanted D. In addition, the previous section shows that the sputtering yields of molecules containing D increase with the concentration of D. This is another source of enhancement of the re-emission of previously implanted D. Considering the error bars, no clear correlation can be deduced between the D concentration and the reflection of the energy.

As for the sputtering, the way of reflection or (a) nc cD=0 at.fr. nc cD≈0.12 at.fr. nc 0.25≤cD≤0.3 at.fr. For the latter group of molecules, it contains all the OD z and Be x D y molecules. Thus at low energy, as for the sputtering data in figure 8, the reflected molecules are mostly OD z . Overall, the most part of the reflected D atoms are single D atoms (≥0.7). However, one can still note the effect of both the D concentration and the ion energy on the fraction of D reflected inside molecules.

c 1 D -Be x O y D z c 2 D -D 2 c 2 D -Be x O y D z c 3 D -D 2 c 3 D -Be x O y D z
A clear positive correlation between the fraction of reflected D inside molecules and the concentration is again observed. For c D = 0 at.fr., this is obvious since the only possible way to be reflected inside a molecule is through the chemical sputtering described in section 3.1.3 (mostly OD dimer as shown in figure 8 (a)). Thus, the re-emission as D 2 is not possible as there is no D before the irradiation. Therefore, large part of the increase between c D = 0 at.fr. and c D > 0 comes from the re-emission of D as D 2 .

This re-emission channel is dominating below 80 eV (figure 12 (b)). Indeed, at low energy, as already mentioned in section 3.1.3, the surface D concentration of D is higher than at high energy since the ions are stopped closer to the surface and spread less in the materials. It follows that the probability of an interaction between the incident D and an adsorbed D is much higher at low energy. In addition, as the surface D concentration is also increasing with the total D concentration, a positive correlation is observed between the fraction of re-emitted D 2 and the concentration of D in the material (figure 12 (b)). It can also break some Be-O bond above a previously implanted D atom which opens a release channel for this D atom. The incident D atom can then interact with this previously implanted D, forming a D 2 molecule which can be released.

The D 2 molecules are not thermally desorbed from the surface. They are produced by the interaction of the incident D atom with another D atom via mechanism similar to the Eley-Rideal or the hot-atom recombination [START_REF] Martinazzo | Hot-atom versus eleyrideal dynamics in hydrogen recombination on ni(100). i. the singleadsorbate case[END_REF]. In the Eley-Rideal mechanism, the incident atom is directly backscattered toward the vacuum taking back a adsorbed D. In the hot-atom recombination, the emission of the D 2 molecule is delayed as the incident atom moves at hyperthermal energy across the surface before the recombination. These two processes are, by definition, only surface processes and are not fully suitable for our energy range. However, the processes are very similar as they involve the recombination of a high energy atom (over the thermal energy) with a D atom at the temperature of the material. For instance, the D atom can be implanted and redirected back to the surface taking with it an adsorbed D.

Above 80 eV, the reflected D atoms inside molecules are mostly reflected in clusters containing Be or/and O atoms (figure 12 (b)). While the fraction of reflected D inside molecules for c D = 0 at.fr. and c D ≈ 0.12 at.fr. stay close together around 10 -3 , it is much high for the highest D concentration. This is a result of the enhanced detachment-induced sputtering described in section 3.1.3 and especially the release of Be x O y D z cluster at the highest concentration. Thus, similarly to what has been concluded for the sputtering, a concentration of 0.25 ≤ c D ≤ 0.3 at.fr. will trigger the dramatic sputtering regime (even in non-cumulative irradiations) leading to a large re-emission of D inside Be x O y D z molecules.

Conclusions

MD simulations have been carried out to simulate the irradiation of (000 1) wurtzite BeO surface by D ions at 300 K with energy from 10 eV to 200 eV. From these simulations, the sputtering yields of Be and O atoms are calculated as well as the reflection coefficient of D atoms from the surface.

Cumulative irradiations show that once a threshold D concentration in the material is exceeded, a 'dramatic' sputtering regime is triggered, leading to detachment of large part of the material. This dramatic sputtering regime is mostly an artifact of the extremely high flux and subsequent large D concentration encountered in cumulative irradiations.

To further analyze the effect of D concentration on the sputtering and reflection properties of BeO, noncumulative irradiations are carried out with three different concentrations. The first one corresponds to a perfect wurtzite BeO cell without D, the second one is the saturation concentration determined experimentally to be 0.12 at.fr. at room temperature. The last one is the threshold concentration to trigger the dramatic sputtering regime around 0.3 at.fr.. The noncumulative irradiations show that the MD sputtering yields increase with the D concentration for a given ion energy: the MD sputtering yields obtained with a D concentration of 0.12 at.fr. are the closest to the experimental values.

Analysis of the chemically bound clusters emitted during the sputtering/reflection events shows that the swift chemical process of O is dominant at low energy up to 80 eV. Above 80 eV, molecules are formed such as BeO dimer through process that does not involve D-induced swift chemical sputtering. The efficiency of the processes leading to these molecules is analyzed. It is found that even in non-cumulative irradiation, some features of the 'dramatic' sputtering are found. Especially, the detachment-induced sputtering causes the rise up of the production of sputtered molecules when the concentration is around 0.3 at.fr.. Table A1: Parameters for the Be-O-H ABOP.

Be-Be [START_REF] Björkas | Interatomic potentials for the Be-C-H system[END_REF] Be-H [START_REF] Björkas | Interatomic potentials for the Be-C-H system[END_REF]

E f j = E BeO+jH -E BeO - j 2 E H2 1 j , (A.12)
in which E BeO+jH is the energy of a BeO cell with j H atoms, E BeO is the energy of a perfect BeO cell of the same size and E H2 is the energy of the H 2 molecule. The formation energies and configurations of j (=1 or 2) H in BeO given by the ABOP were identified by creating more than 100 BeO systems with j randomly positioned H atoms in the lattice. The systems are first relaxed at 300 K for 10 ps and finally quenched in 60 ps at 0K so that the H atoms can find the nearby most stable position. For j = 1, the two configurations which are the most present (151 over 188 initial randomly generated position) involve the H to reside inside an octahedral (Oh) site of the BeO wurtzite lattice formed by 6 O or 6 Be atoms. The description of this site is presented in [START_REF] Hodille | Hydrogen in beryllium oxide investigated by DFT: on the relative stability of charged-state atomic versus molecular hydrogen[END_REF]. In these two configurations, the H is not in the center of the Oh but whether form a O-H bond (called H 1 O-H ) or a Be-H bond (called H Be-H ) (see table A2). Another configuration of a single H, involving a O-H, is also obtained for an H siting in a tetrahedral site but with a higher solution energy compared to the two others.

For the j = 2 case, the most of the final positions show a dissociation of the H 2 molecule and the two H site in two neighboring Oh sites forming O-H bond. However, 54 initial randomly generated position leads to the stabilization of the H 2 molecule in an Oh site oriented in the c crystallographic axis. Most of the final configuration have a solution energy of 1.69 eV/H. However, the formation energy can be further minimized to 1.43 eV/H by distorting a bit the Oh site, but such a configuration represents only 8 of these 54 final positions.

The results of the ABOP are compared to DFT data available in the literature [START_REF] Hodille | Hydrogen in beryllium oxide investigated by DFT: on the relative stability of charged-state atomic versus molecular hydrogen[END_REF][START_REF] Marinopoulos | Defect levels and hyperfine constants of hydrogen in beryllium oxide from hybrid-functional calculations and muonium spectroscopy[END_REF]. BeO is a wide band gap insulator (10.6 eV [START_REF] Weber | CRC Handbook of Laser Science and Technology[END_REF]). In such material, the interstitial hydrogen can act as a donor (H + ), an acceptor (H -) or be in the lattice as a neutral element. Its formation energy depends then on the position of the Fermi level (the chemical potential of the electrons) in the band gap [START_REF] Freysoldt | First-principles calculations for point defects in solids[END_REF]. If both donor and acceptor can exist with the lowest formation energy (for different position of the Fermi level), H has an amphoteric character meaning that it cannot act as a source of carriers. It instead counteracts the existing type of conductivity [START_REF] Marinopoulos | Defect levels and hyperfine constants of hydrogen in beryllium oxide from hybrid-functional calculations and muonium spectroscopy[END_REF]. In BeO, single H has this amphoteric character [START_REF] Hodille | Hydrogen in beryllium oxide investigated by DFT: on the relative stability of charged-state atomic versus molecular hydrogen[END_REF][START_REF] Marinopoulos | Defect levels and hyperfine constants of hydrogen in beryllium oxide from hybrid-functional calculations and muonium spectroscopy[END_REF]. Thus, in the absence of any other defects that would provide carriers, the chemical potential of the electrons lies precisely in the mid-gap [START_REF] Ashcroft | Solid State Physics[END_REF]. We can then take the value of the DFT formation energy of H in BeO for the Fermi energy in the middle of the gap and compare them to our ABOP data. In the DFT calculations [START_REF] Marinopoulos | Defect levels and hyperfine constants of hydrogen in beryllium oxide from hybrid-functional calculations and muonium spectroscopy[END_REF][START_REF] Hodille | Hydrogen in beryllium oxide investigated by DFT: on the relative stability of charged-state atomic versus molecular hydrogen[END_REF], the H + is characterized by the formation of a O-H bond in the material while the H -is characterized by the attraction by H of the Be atoms close to the interstitial H. Thus, the ABOP configuration with a O-H bond is compared to formation energy of H + while the configuration with a Be-H bond is compared to the formation energy of H -solution energy. These values are reported in table A2. The ABOP formation energies are lower than the DFT data. However, the relative error between the DFT data [START_REF] Hodille | Hydrogen in beryllium oxide investigated by DFT: on the relative stability of charged-state atomic versus molecular hydrogen[END_REF] and the ABOP data ranges between 15% (H 2 ) and 40% (H Be-H ). Especially, the most obtained configurations [START_REF] Allouche | First-Principles Study of hydrogen retention and diffusion in beryllium oxide[END_REF] give barriers below 0.9 eV which is inconsistent with the experimental data of Fowler et al. [START_REF] Fowler | Tritium Diffusion in A12O3 and BeO[END_REF] who reported a diffusion barrier of 2.0 eV. This inconsistency has been overcome by calculation of the migration barriers of a H 2 molecule [START_REF] Hodille | Hydrogen in beryllium oxide investigated by DFT: on the relative stability of charged-state atomic versus molecular hydrogen[END_REF] which is calculated to be 1.81 eV in the c direction. This value quantitatively agrees with the experimental data. Thus, we calculate the mobility activation of H 2 molecule by following the mean square displacement R 2 during 10 ns molecular dynamics simulations. Then, the diffusion coefficient D of the H 2 is calculated using the atomistic definition R 2 = 6Dt. An Arrhenius fit is made to the diffusion coefficient data to obtain the diffusion coefficient as D(T ) = D 0 exp(-E D kT ). The simulations are done with a BeO cell containing 1024 atoms which is first relaxed at the desired temperatures for 50 ps. Due to the high energy barrier, the temperature range is 1250-2500K. The diffusion energy and diffusion preexponential factors are E D = 1.19 ± 0.06 eV and D 0 = (7.7 ± 2.9) × 10 -7 m 2 s -1 . The migration energy obtained with the ABOP is lower than the ones calculated by DFT. However, it remains high enough so that the H 2 molecule can be considered immobile in the temperature/time range considered in this study (below 1000 K/10 ns). At these temperature and time ranges, one could have R 2 ≈ 2 Å with the ABOP and R 2 ≈ 0.08 Å with the DFT. These distances stays below the lattice constant which implies that the H 2 will not move at the considered time/temperature ranges whether taking the ABOP or the DFT diffusion coefficient.

Figure 1 :

 1 Figure 1: Snapshots of the BeO cell with a (000 1) surface bombarded at 300 K by 200 eV D ions after, from the left to the right, 350, 750 and 1150 cumulative bombardments. The black circles are the Be atoms, the grey ones are the O atoms and the red ones are the D atoms.

Figure 2 :

 2 Figure 2: Total MD sputtering yields (Be+O) at 300 K obtained after 3000 cumulative impacts and after a limited amount of impacts defined by the condition 0.25 ≤ c max D ≤ 0.30 at.fr.. The MD data are compared to experimental litterature data: exp a : [38] exp b : [40] exp c :[START_REF] Roth | Beryllium selfsputtering: An interpolation of data for D, He, Ne and Ar[END_REF]

2 Figure 3 :

 23 Figure 3: Profiles of the number of Be, O and D atoms after 350 (dashed lines), 750 (solid lines), and 1150 (dashed lines) bombardment of 200 eV D ions. The dotted line at 52 Å represents the position of the initial surface. For sake of clarity, the Be and O atoms are considered together ((N Be + N O )/2) since initially, and all along the simulation, the number of Be stays alsmot the same as the number of O.

  (a) where the slope of the growth of the D concentration decreases as the energy increases.The total sputtering yield obtained for the ≤ c max D ≤ 0.30 at.fr c max D ≈ 0.12 at.fr.

Figure 4 :

 4 Figure 4: (a) Evolution of the maximum D concentration (c D in at.fr.) with the number of cumulative impacts for the different ion energies. (b) Number of cumulative impacts needed to have c D ≈ 0.12 at.fr. (square) and 0.25 ≤ c D ≤ 0.30 at.fr. (circle).cumulative irradiations limited to the number of impacts defined in figure4(b) is plotted on figure2. Compared to the data including all 3000 impacts, the "limited" cumulative data are much closer to the experimental data, especially the high energy ones.However, the low energy data (below 80

Figure 5 :

 5 Figure 5: Total MD sputtering yields (Be+O) at 300K for the non-cumulative irradiations (nc) with three different D concentrations. The sputtering yields for the limited cumulative irradiation (c-limited) are also shown. The MD data are compared to experimental data: exp a : [38] exp b : [40] exp c :[START_REF] Roth | Beryllium selfsputtering: An interpolation of data for D, He, Ne and Ar[END_REF]

Figure 6 :

 6 Figure 6: Ratio Y O /Y Be+O as a function of the ion energy for the different D concentrations in the noncumulative irradiations at 300 K. The black dashed line separates the Be dominated sputtering regime (<0.5) from the O dominated sputtering regime (>0.5).

Figure 7 :

 7 Figure 7: Fration of Be and O sputtered inside chemically bounded clusters Y chem Be+O /Y Be+O as a function of the ion energy for the non-cumulative irradiation at 300 K with different concentration of D.

  nc cD=0 at.f r.

  nc 0.25≤cD≤0.3 at.f r.

Figure 8 :

 8 Figure 8: Fraction of Be and O sputtered inside chemically bounded clusters Y chem Be+O /Y Be+O as a function of the ion energy for the non-cumulative at 300 K with c D = 0 at.fr. (a), c D ≈ 0.12 at.fr. (b) and 0.25 ≤ c D ≤ 0.3 at.fr. (c). The fraction are split over different groups of molecules: OD z , Be x O y D z , Be x H z , O 2 and Be x O y . For sake of clarity, the data point are slightly shifted to prevent an overlapping.

(a) 30

 30 eV -c D = 0 at.fr. (b) 10 eV -c D ≈ 0.12 at.fr.

Figure 9 :

 9 Figure 9: Typical trajectories of incident and sputtered atoms during a SCS event at low energy for c D = 0 at.fr. (a) and c D ≈ 0.12 at.fr. (b). The grey circles are the O atoms, the black ones the Be atoms and the red ones the D atoms.The red lines represents the trajectories of the D atoms and the white lines those of the O atoms. In (b), the trajectory of the incident atom is in orange to distinguish it from other D atoms. The OD 2 molecules at the left of (b) is a result of the cumulative irradiations used to build this cell. None of the sputtering events has led to the sputtering of this molecule.

( a )

 a BeO physical sputteringc D = 0 at.fr. (b) delayed BeO physicalc D = 0 at.fr. (c) delayed SCS -c D ≈ 0.12 at.fr. (d) Detachment-induced sputtering -c D ≈ 0.12 at.fr.

Figure 10 :

 10 Figure 10: Typical trajectories of incident and sputtered atoms during the formation of Be x O y molecules by the BeO physical sputtering process (a), the delayed BeO physical sputtering (b), the delayed SCS sputtering (c) and the detachmentinduced sputtering (d).The grey circles are the O atoms, the black ones the Be atoms and the red ones the D atoms. In each case, the ion energy is 140 eV. For the detachment-induced sputtering (iv), the green trajectory goes from the surface to the bulk and represents the Be atom that attached the sputtered BeO dimer to the rest of the cell.

  cD=0 at.fr. nc cD≈0.12 at.fr. nc 0.25≤cD≤0.3 at.fr.

Figure 11 :

 11 Figure 11: MD reflection coefficient of D atoms, R D (a), and energy, R E (b), at 300 K for different cumulative (labeled c) and non-cumulative (labeled nc) irradiations. The error bars of R D for the 3000impacts cumulative irradiations are directed toward the bottom since the reflection coefficient is increasing as the impacts cumulate making the reported values an upper limit.

Figure 12 : 3 D

 123 Figure 12: Fraction of reflected D atoms inside chemically bound clusters R chem D /R D as a function of the ion energy for non-cumulative irradiations at 300 K with different D concentrations. (a) presents the total fraction of reflected D atoms inside molecules while (b) separates the contribution of D 2 molecules from other Be x O y D z molecules (with x + y > 0). In (b), c 1 D refers to c D = 0 at.fr., c 2 D refers to c D ≈ 0.12 at.fr. and c 3 D

Table A2 :

 A2 Formation energy of H (eV/H) in BeO in different sites. For the ABOP, the H O-H configuration are compared to H + and the H Be-H configuration is compared to H -. For H 2 , the formation energy between brackets is the most stable configuration but the one without brackets is the most obtained one. and single H (H O-H ) agree quantitatively well with the DFT data. The other tested property of the Be-O-H ABOP is the migration energy of the H 2 molecule in the BeO crystal. Calculations of the migration barriers of a single neutral H by Allouche et al.

			DFT
		ABOP	[30]	[29]
	H 2 H 1 O-H H Be-H	1.69-(1.43) 2.04 2.30	1.98 H + 2.82 2.96 -H -3.81 4.54
	of H 2

Acknowledgments This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. Grants of computer capacity from the Finnish Grid and Cloud Infrastructure (persistent identifier urn:nbn:fi:researchinfras-2016072533) are gratefully acknowledged.

Appendix A. O-H potential for Be-O-H applications

We present in this appendix the construction of an interatomic potential for the Be-O-H system. The functional form of the analytical bond order potential (ABOP) has already been extensively discussed [START_REF] Albe | Modeling the metal-semiconductor interaction: Analytical bond-order potential for platinum-carbon[END_REF] and it is only briefly presented here. The total potential energy is obtained by the sum over all atomic bonds:

where the interatomic potential V ij is divided into a repulsive part V R and an attractive part V A . They form together Morse-like pair potentials [START_REF] Morse | Diatomic Molecules According to the Wave Mechanics[END_REF] V

D 0 and r 0 are the dimer bond energy and length. The interactions between the atoms are usually restricted to the nearest-neighbors by the cutoff function:

R is the cutoff distance and D defines the width of the cutoff region. The maximum interaction range is then R+D. The three-body interactions are included in the bond-order term, b ij which includes the angles between the bond

where

with

θ ijk is the angle between atoms i, j and k. α ijk and ω ijk and are fitting parameters. The angular function g ik is define as

) where γ ik , c ik , d ik and h ik are adjustable parameters. g ik (θ ijk ) minimum is γ > 0 and it is reached for h ik =cos θ ijk . In the simulations involving high kinetic energies, the repulsive part of the potential in this formalism is unphysically weak. Thus, it has to be modified to describe short-range interactions more accurately. This is done by joining the universal repulsive Ziegler-Biersack-Littmark potential [START_REF] Ziegler | The Stopping and Range of Ions in Matter[END_REF], V ZBL (r ij ), with the original potential V ij as

which ensures a smooth transition between the original potential and the repulsive ZBL potential. The transition range and distance are defined by the parameter b f and r f respectively. They are chosen such the equilibrium properties are unchanged.

Appendix A.1. Dimer O-H and molecules H n O

The presented ABOP is mainly devoted to simulations of BeO irradiation by plasma ions in fusion devices. We are especially interested in the erosion of this material by hydrogen isotops ions coming from the plasma. Obviously, no liquid water is expected on top of any plasma-facing materials. Thus, the complex interaction of H with O in the liquid water through hydrogen bond is not needed in this version of the ABOP, and we can only stick to the reproduction of the O-H dimer bond as well as reasonably good energetics of H n O molecules, knowing that H 2 O is the most stable.

Ni et al. [START_REF] Ni | A reactive empirical bond order (REBO) potential for hydrocarbonoxygen interactions[END_REF] used density functional theory (DFT) method with the B3LYP exchange-correlation potential to parametrize a reactive empirical bond order potential for hydrocarbon-oxygen interactions.

They calculated a bond length for the O-H dimer of 0.9752 Å and a bond energy of -4.544 eV which are both in good correspondence with the experimental ones (0.971 Å and -4.432 eV respectively [START_REF] Sutton | Tables of Inter-Atomic Distances[END_REF]). We use these DFT values for the ABOP value r 0 and D 0 . Ni et al. [START_REF] Ni | A reactive empirical bond order (REBO) potential for hydrocarbonoxygen interactions[END_REF] also calculated the potential energy surface for the O-H dimer by single-point DFT calculations. We used this potential energy surface to obtain both S and β for the O-H dimer. These parameters are reported in table A1. With these parameters, we are confident that the potential reproduce correctly the O-H dimer which is our most important concern since it is the most probable product of a chemical sputtering of oxide by hydrogen isotope ions. It can be noted that there is a different set of parameters for O-H and for H-O. The difference is only on the parameters of the angular function g ik (θ). When i = H and k = O, for any third atom in the interaction range of that H atom, g HO (θ) will be constant at a high value preventing the formation of unphysical molecules such as HO n>1 . This is the opposite when i = O as molecules of the form OH n>1 exist. The angular function g(θ) sets the most stable configuration of a H n O molecules for n ≥ 2. It has to be pointed out here that with the current formalism, for any H n O molecules, the energy of the O-H bonds are higher than the one for the Dimer. Indeed, the lowest value of g(θ) is γ > 0 which implies that the bond order term b ij cannot be higher than 1. However, as calculated by Ni et al. [START_REF] Ni | A reactive empirical bond order (REBO) potential for hydrocarbonoxygen interactions[END_REF], the O-H bond in the H 2 O molecule is stronger than in the O-H dimer. We thus choose a value of γ so that the energy of the O-H bond in the H 2 O molecule is approximatively the same as in the O-H dimer with a difference of 0.001 eV. The most stable configuration is the H 2 O molecule in which the angle between the two O-H bonds is 104.5 • . Thus, we could have used a value of h so that this angle gives the most stable configuration for H 2 O. However, doing so, the H 3 O configuration would also have angles of 104.5 • between each of its O-H bonds leading to a higher cohesive energy (potential energy per atom) than the H 2 O, which is not realistic. Thus, we decided to use h = 1.0 so it maximizes the angle between any O-H bonds. Doing so, the H 2 O is the most stable molecule among the H n O Appendix A. [START_REF] Brezinsek | Plasma-surface interaction in the be/w environment: Conclusions drawn from the jet-ilw for iter[END_REF]

. H interstitial in BeO

As test properties of the Be-O-H ABOP, we used the formation energies of interstitial hydrogen in different sites inside the BeO wurtzite structure. Recently, DFT calculation of H in BeO have shown that H 2 could form in the BeO lattice [START_REF] Hodille | Hydrogen in beryllium oxide investigated by DFT: on the relative stability of charged-state atomic versus molecular hydrogen[END_REF]. Thus, we also investigated the formation energies of such molecules in the presented ABOP. In the ABOP, the formation energy E f j (in eV/H) of j hydrogen atoms in the lattice is given by