Modelling of deuterium diffusion and thermal desorption in tungsten exposed to high-flux deuterium plasma

In nuclear fusion, hydrogen isotopes transportation in tungsten has been identified as a crucial process which could affect various key processes in the operation of a fusion device, such as tritium self-sustaining and operational safety. To enhance the understanding of deuterium diffusion and trapping in millimeter-thick tungsten under high-flux plasma conditions, we carried out a systematic simulation investigation on the deuterium plasma exposure and subsequent thermal desorption using the TMAP7 code based on experimental results. The simulation reveals that the defect concentration with a low de-trapping energy , as estimated solely from simulating thermal desorption, is lower than its actual density due to the temperature history of plasma exposure and subsequent cooling processes. Moreover, the growth rate of deuterium mobile concentration exhibits nonmonotonic behavior over time during plasma exposure, which strongly depends on the ratio of defect concentration with a high de-trapping energy to that with a low de-trapping energy. Additionally, we observed that the maximum deuterium diffusion depth does not follow a proportional relationship with the square root of exposure time, which is attributed to the growing number of plasma-induced defects and the reducing injection coefficient with increasing fluences. This work provides valuable insights into the understanding of deuterium transportation in tungsten during high-flux plasma exposure which could contribute to understanding and modelling of hydrogen isotopes recycling in future fusion devices.


Introduction
Transportation of hydrogen isotopes (HIs) in materials has been attracted much attention in fields like nuclear energy and hydrogen energy in recent years.In the context of nuclear fusion, understanding HIs transportation is an important issue from the viewpoint of economy and security [1,2].It could affect many key processes in the operation of a fusion device: the balance of HIs into and out of the plasma-facing materials (PFMs) determines the stability of the fuel recycling process and therefore affects the stable operation of the core plasma, the retention of HIs in the bulk reduces the self-sustaining ability of tritium, and the permeation of tritium into the coolant leads to the formation of tritiated water which is extremely difficult to be dealt with because of its acidity.Consequently, comprehending HIs transportation is important to the modelling of the related processes and the proposing of potential solutions to HIs-related problems.
Tungsten (W) is regarded as one of the most promising PFMs for future tokamak devices due to its advantageous properties of high melting point, high thermal conductivity, low sputtering yield and low hydrogen isotopes retention [3,4].During the operation of devices, the edge plasma freed from magnetic field bondage would impinge on the surface of PFMs causing HIs transportation in tungsten.Due to the large area of tungsten surface in the divertor and potential usage of tungsten as a first wall material, HIs transportation in tungsten is of significance to the operation of a commercial Tokamak device.HIs transportation in tungsten can be briefly summarized as follows.At beginning, some HIs bombarding the tungsten surface will be directly backscattered into the edge plasma, and others will be implanted as interstitial atoms in a shallow subsurface layer which is usually in the order of few nanometers.A reflection coefficient is therefore defined as the ratio of the reflected number of HIs to the total impinging number, which depends on factors such as particle species [5], particle incident energies [6] and surface roughness [7,8].Then, the implanted HIs would diffuse in two directions, one towards the exposure surface where the atoms recombine into molecules and release, the other towards the bulk.Due to the low HIs solubility in tungsten [9], if the implanted flux is greater than the rate at which HIs diffuse out of the implanted zone, a local supersaturation of HIs will appear [10].During HIs diffusion in the bulk, HIs would encounter and be trapped by irradiation-induced as well as intrinsic defects [11], such as vacancies, dislocations, grain boundaries, bubbles, and voids.These trapped HIs contribute to the HIs inventory in tungsten.Interstitial HIs not lost to traps would diffuse deeper and eventually reach the backside of the material when permeation happens.Factors affecting the transportation include the material temperature, the defect concentration, and high-diffusivity paths such as grain boundaries and so on [12][13][14][15].HI transportation is usually described using the macroscopic rate theory equation model.This model is early discussed in [16] and applied to a number of metals like steel [17,18] and tungsten [19][20][21], has led to the development of codes such as TMAP7 [22], DIFFUSE [23], and PIDAT [24].Complex physical and chemical processes are simplified by adjusting the governing equations, boundary conditions and fundamental parameters in these codes.They have been extensively used to make predictions for future nuclear devices and simulate experimental phenomena including D/He exposure as well as thermal desorption process [25][26][27].
In recent years, HIs solution and diffusion in tungsten [28][29][30], copper [31] and other metals [32,33] have been investigated by plasma-driven permeation and gas-driven permeation experiments.Due to the low HIs diffusion coefficient in tungsten, the thickness of the samples used in these experiments is very thin as tens of microns.Measuring millimeter-thick tungsten samples using these methods is challenging.Furthermore, the plasma flux in the plasma-driven permeation experiment is usually a few magnitudes lower than that expected in the first wall or divertor, which limits the understanding of HIs diffusion under high-flux plasma conditions.In our previous work, thermal desorption spectroscopy (TDS) of deuterium (D) in millimeter-thick tungsten under high-flux plasma exposure with fluences up to 1.0 × 10 28 m −2 was measured, and depth profiles of D and defects were estimated by the TDS simulation using TMAP7 [34].In this work, both D diffusion during plasma exposure and thermal desorption after exposure are simulated using TMAP7 based on the previous TDS data, which demonstrates D transportation characteristics in tungsten under high-flux plasma conditions.We find that the defect concentration with a low de-trapping energy, as revealed by the experimental TDS, is underestimated compared to its actual density due to the temperature history of plasma exposure and subsequent cooling process.Besides, the growth rate of D mobile concentration exhibits nonmonotonic behavior over time during plasma exposure, which strongly depends on the defect concentration with a high detrapping energy.At last, the temporal evolution of the maximum D diffusion depth is determined, which shows both the enhanced concentration of plasma-induced defects and the reduced injection coefficient could slow down the D macroscopic diffusion.

Experimental details
A series of D plasma exposures, which are taken for modelling, were executed on recrystallized W samples with exposure durations of 2, 8, 16 and 64 h at 500 K in the previous work [34].The experimental parameters are summarized in table 1.At least one month after exposure, the thermal desorption of D atoms retained in the samples was measured by TDS at Beihang University [35].During TDS, the samples were heated linearly up to 1273 K with a rate of 1 K s −1 at 1 10 Pa.

6
´-HD and D 2 signals were monitored by a quadrupole mass spectrometer (QMS) and the absolute sensitivity was calibrated using a standard leak.

TMAP simulation
TMAP7 is a 1-D reaction-diffusion code which is based on the rate theory equation model.It has a great ability to deal with many kinds of problems, such as radioactive decay, heat source and heat diffusion, etc. Besides, it allows for the consideration of up to three types of defects with different de-trapping energies and defined concentration depth profiles.
In this study, the D plasma exposure processes with the four different fluences and the subsequent D 2 thermal desorption are simulated to give information about D diffusion and trapping by defects.In order to fit the experimental TDS, some parameters and conditions are adjusted as summarized in table 2. The entire sample is divided into 269 diffusion segments, and the defect concentration depth profiles are adjusted.Based on the temperature history of the experiments, the simulation can be roughly divided into three stages which are isothermal heating at 500 K during exposure, cooling down to 300 K within tens of minutes after exposure, and heating linearly up to 1273 K in 973 s during thermal desorption.An implantation depth of 2.6 nm is used according to Eckstein [5].The D injection coefficient (one minus the reflection coefficient) is adjusted as discussed in section 3.1 which shows a tendency to decrease with increasing fluence.The diffusion activation energy taken from Frauenfelder [36] is 0.39 eV, and the heat of solubility is 1.04 eV [37].In TMAP7, a defect can only be occupied by a single HI.The concentrations of defect1 and defect2 as intrinsic defects are both set at 1 10 at.fr.

6
´and the de-trapping energies of them are 1.10 eV and 1.35 eV, respectively.The depth profile of defect concentration and the D injection coefficient remain unchanged over time in a single simulation using TMAP7.

Results and discussion
3.1.Simulated D thermal desorption by modelling plasma exposure using concentration depth profiles of defects obtained by only fitting experimental TDS The experimental and simulated TDS of the four exposed samples have been reported in the previous work [34] as shown in figure 1.The experimental TDS shows two peaks in all cases with the initial desorption temperature at ∼400 K.As fluence increases, the first main desorption peak remains fixed at 654 K, while the second peak varies from 744 K to 823 K.The simulated TDS using TMAP7 reveals two types of defects with de-trapping energies of 1.10 eV (defetct1) and 1.35 eV (defect2), respectively.Notably, a sharp drop in the concentration depth profiles of defects obtained by only fitting experimental TDS is suggested at ∼40 μm as shown in figure 1(a), which could be attributed to the effect of the first layer of grain boundaries parallel to the exposure surface.
In the simulation of plasma exposure and subsequent thermal desorption, the concentration depth profile of the defect2 in [34] is employed because D trapped in defect2 with a high de-trapping energy is less likely to desorb at the low exposure temperature of 500 K or during cooling [38].Thus, the concentration depth profile of D trapped in defect2 in figure 1(a) could represent the actual concentration profile during exposure.In the simulation, the injection coefficient is adjusted so that the preset defect2 is just completely filled, resulting in injection coefficients of 0.014, 0.0065, 0.006, and 0.0031 for RECW-2h, RECW-8h, RECW-16h, and RECW-64h, respectively.It shows a decreasing trend, which may be related to the surface morphology modifications due to plasma exposure.Scanning electron microscopy (SEM) analysis reveals that the number density of blisters gradually develops with the increasing fluence, and burst blisters are observed in RECW-16h, and RECW-64h [34].It has been reported that the cavities of the burst blisters connected to the sample surface can provide a shortcut for the release of D atoms, which results in a reduced incident flux of D atoms into the bulk [39].In addition, protrusions and blisters can increase the surface area as well as the reflection coefficient of HIs on the surface [39].
After the concentration depth profiles of defect2 and the injection coefficients are determined, the simulation of plasma exposure and subsequent thermal desorption is first performed using the concentration depth profiles of defect1 in [34] and the simulated TDS is shown in figure 1(c).It is obvious that the intensity of the first peak at 654 K is significantly lower than the experimental value, indicating that the preset defect1 concentration obtained by only fitting the experimental TDS is less than the actual one formed during plasma exposure.This discrepancy arises because a large portion of defect1 could not be filled at the exposure temperature of 500 K due to its low de-trapping energy.Besides, there is a cooling time for the temperature to return to RT after terminating the plasma exposure.Consequently, a part of the D atoms trapped in defect1 near the exposure surface could detrap and release outside the bulk.During the cooling period of the simulation, the desorbed D atoms from RECW-2h, RECW-8h, RECW-16h and RECW-64h accounts for 7%, 4%, 3% and 2% of the D retention at the time when terminating the plasma exposure.The defect1 concentration should therefore be modified, as discussed in the following section.

Simulated D thermal desorption by modelling plasma exposure using adjusted concentration depth profiles of defects
In order to better match the experimental TDS, the concentration depth profiles of defect1 are adjusted by considering D occupation in defects during exposure and D desorption during cooling, as shown in figure 2(a).It should be noted that the concentration of defetct1 within the first 20 μm is adjusted, especially the first 8 μm as D atoms retained near the surface are most likely affected.The maximum concentration increases from 1.9 10 4 ´to at fr 6.0 10 ..,

3
´-For RECW-2h, RECW-8h, RECW-16h, and RECW-64h, respectively, as compared to [34]. Figure 2(b) shows the simulated TDS, which is more consistent with the experimental TDS compared to  figure 1(c).The adjustments made to the concentration depth profiles of defect1 lead to a better agreement between the simulated and experimental D thermal desorption.
The nature of the two defects cannot be determined by TDS but could be preliminarily discussed.Some researchers suggest that the de-trapping energy for H in dislocations in tungsten varies in the range of 0.85-1.25 eV [37,[40][41][42][43].And the de-trapping energy of 1.10 eV used in this work is very close to the interaction energy between H and edge dislocation calculated by the first-principles method [44].For defect2, 1.35 eV is in reasonable agreement with the de-trapping energy of 1.29-1.6eV for vacancy-type defects [21,[45][46][47][48]. Besides, it has been reported that the de-trapping energy of H in grain boundary is comparable to that of a vacancy [49,50].Therefore, defect1 may be a dislocation-type defect and defect2 could be vacancy-type defects such as vacancies or grain boundaries.Due to the limitation of TMAP7, a defect can only trap a HI.It has been shown that density functional theory calculations that a vacancy can trap multiple HIs with varying binding energies [47,51].The defect concentration set in the code therefore must be higher than that in the reality to compensate the amount of trapped D atoms.And the de-trapping energy used in the code should be the statistical average of de-trapping energies when multiple HIs could detrap from one defect.As reported in [47,52], with the increasing number of desorbed D atoms, the de-trapping energy of the defect gradually increases, leading to a shift in the desorption peak towards higher temperature.The defect distribution depth set in TMAP7 may be deeper than the actual depth to achieve a match with the high-temperature signal in experimental TDS.

Simulation of D transportation during plasma exposure
The 'mobile concentration' in the unit of Dm , 3 -defined as the concentration of the interstitial D atoms as solute which are not trapped at a depth with a certain thickness, could help reveal the characteristics of D transportation during plasma exposure.The temporal evolutions of the mobile concentration at different depths within the first 60 μm are shown in figures 3 and 4, using the concentration depth profiles of defects in sections 3.1 and 3.2, respectively.
Despite imperfect TDS fitting in section 3.1, the mobile concentration as a function of time shows interesting results.In figure 3, the mobile concentration at shallower depths grows rapidly and eventually reaches 'saturation' as time increases.This 'saturation' state is not an equilibrium state, but a stationary state where the net flux ( net in out G = G -G ) of D atoms remains constant over time.The time taken for the mobile concentration at the same depth in the four samples to reach a stationary state increases with the increasing fluence, which is attributed to the increase of the defect concentration with the increasing fluence.Besides, it is found that the growth rate of the mobile concentration as a function of time isn't monotonic but demonstrates many steps.In general, the number of D atoms implanted into W should gradually increase over time, but these steps at certain depths indicate that the net G suddenly slows down in a faster manner at a given time.The question arises: Where does out G go?The answer lies in the D trapping by the defects at the near-neighbour depths.Take the mobile and D trapping concentration at 2 μm in figure 3(a2) as an example.In region I, it is obvious that D atoms are more likely being trapped by defects rather than dissolving as solutes.And the time at which the mobile concentration starts to increase is later than that for the D trapping concentration.In region II, when a step at 2 μm appears for the first time, the mobile concentrations in the next microns are negligible, but D atoms are trapped by defects at 3 and 4 μm.The trapping rate is highest at the depth of 3 μm and remains constant for a period of time until 40%-50% of the defects are filled by D. In region III, as the step at 2 μm ends, the growth rate of the mobile concentration at 2 μm resumes its normal value.The appearance of the steps is therefore attributed to solute D atoms at 2 μm diffuses into 3 μm and being trapped by the defects.Due to the strong attraction of defects to D atoms, out G at 2 μm suddenly increases.
However, as shown in figure 4, the steps become less pronounced after increasing the concentration of defect1 within the first 20 μm.This suggests the sudden increase out G at a certain depth due to D trapping at next depths is prevented.In fact, it is the increased defect concentration of defect1 which traps D atoms at a certain depth that prevents D trapping by defects at next depths.Under extreme conditions, when the concentration of the defect with a low de-trapping energy is extremely large, the steps would disappear, and when the concentration of the defect with a high de-trapping energy is extremely high, the steps will be more apparent where the growth rate of the mobile concentration over time is close to zero.

Diffusion front
As mentioned in the introduction section, D diffusion in millimeter-thick tungsten under high-flux plasma exposure is less investigated due to limitations in the experiment conditions, but can be studied in this work by analyzing the maximum depth that D could reach under the four fluences.
A 'diffusion front' is defined as the depth where the mobile concentration of D atoms reaches a minimum of 1.6 10 at.fr.

9
´or zero at the end of exposure.The following diffusion fronts are obtained according to the simulation of plasma exposure with the modified concentration depth profiles of defects in section 3.2.In case of several factors.First, reducing injection coefficients in section 3.1 indeed reduces the D incident flux, leading to a decrease in the mobile concentration and thereby shortening the final diffusion front.Second, due to the increasing plasma-induced defect concentration with the increasing fluence as revealed by TDS, the movement of solute D atoms towards the bulk slows down.In order to verify the two effects , two additional groups of simulations are carried out as follows.In the first group, the same uniform defect concentration of 1.0 10 at.fr.

4
´is set in the four samples, with injection coefficient varied as 0.0031, 0.006, 0.0065 and 0.014.In the second group, the injection coefficient in the four samples is set to be identical as 0.014, and the defect concentrations are varied from 1.0 10 at.fr. 4 ´to 2.0 10 at.fr. 4 ´-As presented in figure 5(b), the diffusion front is deeper if the injection coefficient is larger.In case of the higher defect concentration, the diffusion front is visibly shallower.Additionally, as the defect concentration and the injection coefficient remain constant, the expected square-root time dependence of diffusion front is observed with the increasing fluence/time.Furthermore, it has been reported that stress and strain in the material could alter HIs diffusion [15,53].The defects formed by D plasma exposure could create a complex strain field in the bulk, which might reduce the diffusion coefficient of D atoms under certain conditions.However, due to the limitation of the code, this factor could not be evaluated in this work.A new code is under development to include the contribution of strain or stress to HIs transporation.

Conclusion
By simulating the D plasma exposure and thermal desorption using TMAP7, characteristics of D diffusion and trapping in tungsten under high-flux plasma exposure are investigated.Two types of defects with de-trapping energies of 1.10 eV (defect1) and 1.35 eV (defect2) are used and the injection coefficients of the four samples are adjusted to be 0.014, 0.0065, 0.006, 0.0031, respectively.Besides, it is found that the defect concentration with a low de-trapping energy, as revealed by only simulating thermal desorption , is lower than its actual density due to the temperature history of plasma exposure and subsequent cooling processes.And defetct1 within the first 20 μm are most likely affected.In addition, the growth rate of D mobile concentration exhibits nonmonotonic behavior over time during plasma exposure, which strongly depends on the ratio of defect concentration with a high de-trapping energy to that with a low de-trapping energy.At last, the temporal evolution of the maximum D diffusion depth is determined by locating the position with a minimum mobile concentration of 1.6 10 at.fr. ´or zero, which shows both the increasing concentration of plasma-induced defects and the reduced injection coefficient could slow down the D macroscopic diffusion.

Figure 1 .
Figure 1.(a) Concentration depth profiles of defects obtained by only fitting experimental TDS [34].(b) Comparison between the simulated TDS and the experimental TDS [34].(c) Simulated TDS using the concentration depth profiles of defects in [34] but considering both plasma exposure and thermal desorption processes.The de-trapping energies of the two defects are 1.10 eV and 1.35 eV, respectively.

Figure 2 .
Figure 2. (a)Modified concentration depth profiles of defects used in the simulation.(b)Comparison between the simulated TDS and the experimental TDS.

Figure 4 .
Figure 4. D mobile concentration at different depths within the first 60 μm in the four samples during the simulation of plasma exposure using the concentration depth profiles of defects in figure 2(a).(a1), (b1), (c1), (d1) are the mobile concentration during the whole exposure period of 2, 8, 16 and 64 h, respectively.(a2), (b2), (c2), (d2) are the detailed images of the mobile concentration and D trapping concentration for the first 7 μm.

Table 1 .
Experimental parameters of samples.

Table 2 .
Summary of parameters used in TMAP7 code.