Understanding hydrogen retention in damaged tungsten using experimentally-guided models of complex multispecies evolution

Here we use the extension of the stochastic cluster dynamics method (SCD) [63] to perform all simulations. SCD is a stochastic variant of the mean-field rate theory technique, alternative to the standard implementations based on ordinary differential equation (ODE) systems, that eliminates the need to solve exceedingly large sets of equations and relies instead on sparse stochastic sampling from the underlying kinetic master equation. Rather than dealing with continuously varying defect concentrations C_i, SCD evolves integer-valued defect populations N_i in a finite material volume Ω, thus limiting the number of 'active' ODEs at any given moment. Mathematically, SCD recasts the standard ODE system:

$$\begin{equation}\small \frac{dC_i}{dt} \!=\! g_i+\sum_q\left( \!\!\sum_j s_{jq}C_j-s_{iq}C_i\!+\!\right)\!+\!\sum_{j}\!\left[\left(\!\sum_{k}\mathcal{k}_{jk}C_k-\mathcal{k}_{ij}C_i\!\right)C_j\!\!\right] \end{equation} \tag{ 1 }$$

into stochastic equations of the form:

$$\begin{equation}\small \frac{dN_i}{dt}\! =\! \tilde{g}_i+\!\sum_q \!\left(\!\!\sum_j \tilde{s}_{jq}N_j-\tilde{s}_{iq}N_i+\!\!\right)+\!\!\sum_{j}\!\left[\!\left(\!\sum_{k}\tilde{\mathcal{k}}_{jk}N_k-\tilde{\mathcal{k}}_{ij}N_i\!\!\right)N_j\!\!\right] \end{equation} \tag{ 2 }$$

where the subindices i, j, and k refer to a given defect species, and q runs over all potential defect sinks. The set $\{\tilde{g},\tilde{s},\smash{\tilde{\mathcal{k}}}\}$ represents the reaction rates of 0th (insertion), 1st (thermal dissociation, annihilation at sinks), and 2nd (binary reactions) order kinetic processes taking place inside Ω, and is obtained directly from the standard coefficients $\{g,s,\mathcal{k}\}$ as:

$$\begin{equation*} \tilde{g}\equiv g\Omega,~\tilde{s}\equiv s,~\tilde{\mathcal{k}}\equiv k\Omega^{-1} \end{equation*}$$

The value of Ω chosen must satisfy

$$\begin{equation*} \Omega^{\frac{1}{3}}>\ell \end{equation*}$$

where $\ell$ is the maximum diffusion length l_i of any species i in the system, defined as:

$$\begin{align} \ell = \text{max}_i\{l_i\} \end{align} \tag{ 3 }$$

$$\begin{align} l_i = \sqrt{\frac{D_i}{R_i}} \end{align} \tag{ 4 }$$

Here, D_i and R_i are the diffusivity and the lifetime of a mobile species within Ω. The above expression is akin to the stability criterion in explicit finite difference models. From equation (2), $R_i = \tilde{s}+\sum_j\tilde{\mathcal{k}}_{ij}N_j$. The system of equations (2) is then solved using the kinetic Monte Carlo (residence-time) algorithm by sampling from the set $\{\tilde{g},\tilde{s},\tilde{\mathcal{k}}\}$ with the correct probability and executing the selected events. Details of the microstructure such as dislocation densities and grain size are captured within SCD in the mean-field sense, i.e. in the form of effective sink strengths for hydrogen atoms. For example, the value of s_iq in equation (1) for dislocation and grain boundary defect sinks is:

$$\begin{equation} s_{iq} = s_{\rm d}+s_{\rm gb} = \left(\rho_{\rm d}+\frac{6\sqrt{\rho_{\rm d}}}{d}\right)D_i \end{equation} \tag{ 5 }$$

where $\rho_{\rm d}$ is the dislocation density and d is the grain size (in this case, q = 1 ≡ 'd' and q = 2 ≡ 'gb').

SCD has been applied in a variety of scenarios not involving concentration gradients [63, 64]. However, equation (1) must be expanded into a transport equation (i.e. a partial differential equation, or PDE) by adding a Fickian term of the type:

$$\begin{equation} \frac{dC_i}{dt} = \nabla\left(D_i\cdot\nabla C_i\right)+f(t; C_1,C_2,C_3\ldots) \end{equation} \tag{ 6 }$$

where D_i is the diffusivity of species i, and $f(t; C_1,C_2,C_3\ldots)$ is used for simplicity to represent all of the terms in the r.h.s. of equation (1). To cast equation (6) into a stochastic form, the transport term must be converted to a reaction rate in the finite volume Ω. As several authors have shown, this can be readily done by applying the divergence theorem and approximating the gradient term in terms of the numbers of species in neighboring elements [52, 65]. In the most general 3D case:

$$\begin{equation*} \nabla\left(D_i\cdot\nabla C_i\right)\rightarrow D_i\sum_{\beta} A_{\alpha\beta}\frac{N_i^{\beta}-N_i^{\alpha}}{\Omega_{\alpha}l_{\alpha\beta}} \end{equation*}$$

where α is a subindex that represents the current element, and β represents the neighboring elements (β = 1, ..., 6 in 3D). $A_{\alpha\beta}$ is the dividing surface between neighboring elements α and β, $l_{\alpha\beta}$ is the distance between their centers, and $\Omega_\alpha$ is the current element's volume. Figure 2 shows a schematic diagram of the geometry of two generic volume elements. In this fashion, each term of the above sum represents the rate of migration of species i from volume element α to β, which can be now added to the r.h.s. of equation (2) and sampled stochastically as any other event using the residence-time algorithm.

Zoom In Zoom Out Reset image size

The extension of equation (2) from an ODE system into a PDE system can be written for species i in volume element m as:

$$\begin{equation} \begin{split} \frac{dN_i^{\alpha}}{dt}\! =\! D_i\!\sum_{\beta}A_{\alpha\beta}\!\frac{N_i^{\beta}-N_i^{\alpha}}{\Omega_{\alpha}l_{\alpha\beta}}\!+\!\tilde{g}_i\!+\!\sum_q\left(\!\! \sum_j \tilde{s}_{jq}N^{\alpha}_j\!-\!\tilde{s}_{iq}N^{\alpha}_i\!\!\right)+\\ +\sum_{j}\left[\left(\!\!\sum_{k}\tilde{\mathcal{k}}_{jk}N^{\alpha}_k-\tilde{\mathcal{k}}_{ij}N^{\alpha}_i\right)N^{\alpha}_j\right] \end{split} \end{equation} \tag{ 7 }$$

For regular cubic meshes, or in 1D, the Fickian term simply reduces to $D_i\frac{N_i^{\beta}-N_i^{\alpha}}{l^2}$, where l is the element size. This is the approach employed here.

As described by Marian and Bulatov [63], SCD is ideally suited to deal with the probabilistic aspects of irradiation damage. Once a suitable source function is constructed, one can accurately sample it to create defect distributions. Establishing the source term requires the calculation of recoil energy distributions as a function of depth to be used in SR-SCD calculations. Primary damage is generated from a W-recoil distribution obtained using the SRIM software package [66] for approximately 1000 Cu ions with the same incident energy of 3.4 MeV. According to the calculations, only about 43.5% or 1.48 MeV of the total incident energy is expended on lattice damage, with the rest lost on ionization via electronic stopping. We call this energy E_D. The recoil spectra are given as cumulative distribution functions $C(E_{\rm PKA})$, seen in figure 3 as a function of target depth. The average PKA energy from this distribution comes out to be $\langle E_{\rm PKA}\rangle = 2.18$ keV. These functions are then used to obtain random samples of the primary knock-on atom (PKA) energies $E_{\rm PKA}$ by solving $E_{\rm PKA} = C^{-1}(\xi)$, where ξ is a random number uniformly distributed in (0, 1]. Note that $C(E) = \int_0^E P(E)dE$, where P(E) is the normalized recoil energy spectrum.

Zoom In Zoom Out Reset image size

The rate of Cu-ion insertion is calculated as:

$$\begin{equation} r_{\rm ion} = \rho_a^{-1}\phi \gamma \end{equation} \tag{ 8 }$$

where φ is the ion flux, ρ_a is the atomic density of W, and γ is a damage function calculated within SRIM and measured in units of displacements per ion per unit length, and corresponds to the dashed line shown in figure 3. Once an ion impact event is selected in the simulation volume, a sequence of PKA energies is obtained following the procedure just described until the sum of all sampled recoil energies reaches E_D. The insertion coefficients g are obtained by sampling from discrete distributions parameterized to reproduce sub-cascade statistics collected from hundreds of MD cascade simulations at different temperatures [67, 68]. Specifically, the number of Frenkel pairs produced as a function of PKA energy are:

$$\begin{equation} N_F = \Bigg\{ \begin{matrix} & 1.15\times10^{-2} \left(E_{\rm PKA} \right)^{0.74}, & \text{if}~ E_{\rm PKA} \lt 43~\rm{keV}\\ & 1.89\times10^{-5} \left( E_{\rm PKA}\right)^{1.34}, & \text{otherwise} \end{matrix} \end{equation} \tag{ 9 }$$

where $E_{\rm PKA}$ is given in eV. The dependence of N_F with temperature was seen to be weak and is not considered in this study. We set the minimum value of $E_{\rm PKA}$ to produce a stable Frenkel pair inside the material to be 620 eV [64]. For their part, the fractions of clustered SIAs and vacancies $f_{c}^{SIA}$, $f_{c}^{V}$ depend strongly on, respectively, PKA energy and temperature [69]:

$$\begin{align} f_{c}^{SIA} = 0.018\,5 \left(E_{\rm PKA} \right)^{0.326} \end{align} \tag{ 10 }$$

$$\begin{align} f_{c}^{V} = 0.625 - 1.750\times10^{-4} T_{\rm irr} \end{align} \tag{ 11 }$$

where $T_{\rm irr}$ is the absolute irradiation temperature.

The s_ij coefficients used in equation (7) dictate the propensity of a species i to be absorbed at a sink of type j. Generally, these sinks are intrinsic microstructural heterogeneities, such as grain boundaries (GBs), precipitates, free surfaces, and/or dislocations (refer to equation (5)). In the present study, we consider pure W single crystal surfaces, which eliminates the first two types of sinks. Also, since the surface of the material is explicitly resolved spatially, one need not to consider free surfaces as independent sinks (the simulations themselves capture the migration of defects to the surface layer). This leaves dislocations as the only microstructural feature that can absorb defects. Sinks are characterized by their sink strength, which for dislocations is directly the dislocation density ρ_d. The absorption rate of mobile defects i by dislocations is then:

$$\begin{equation} r_i^{\rm abs} = z_i\rho_{\rm d}D_i \end{equation} \tag{ 12 }$$

where z_i is an appropriate bias factor that captures the increased propensity of some defects over others to be absorbed by dislocations (e.g. SIAs vs. vacancies). Consequently, the number of defects absorbed is simply:

$$\begin{equation*} \frac{dN^d_i}{dt} = r_i^{\rm abs}N_i \end{equation*}$$

where the superscript 'd' is used to indicate that particular fraction of defects is trapped at the dislocation.

While the absorption of mobile defects and H atoms by dislocations is spontaneous, the inverse process—i.e. emission—is thermally activated with a dissociation energy of E_d (for H atoms, it is ≈ 0.6 eV according to several estimates [70])). After Friedel [71], this emission rate can be expressed as:

$$\begin{equation} r^{\rm em}_i = \left(\frac{2\pi b\rho_{\rm d}\Omega}{a_0^2}\right)\nu_0\text{exp}\left(-\frac{E_d}{kT}\right)\left(N^d_i-N^0_i\right) \end{equation} \tag{ 13 }$$

The first term (in parentheses) on the r.h.s. of the above expression is a geometric factor that gives the number of possible emission sites from a cylindrical 'tube' around dislocation segments, with a₀ the lattice parameter and $b = a_0\sqrt{3}/2$. ν₀ is an attempt frequency (e.g. the Debye frequency of the material), and $N^0_i$ is the equilibrium concentration of species i. Assuming thermal equilibrium,

$$\begin{equation*} N^0_i = \left(\rho_a\Omega\right)\exp\left(-\frac{E^f_i}{kT}\right) \end{equation*}$$

where $E^f_i$ is the formation energy (for H atoms, this is equivalent to the heat of solution) of species i. In practice, emission from dislocations is only enabled for H atoms, such that—after inserting the expression for $N^0_{\rm H}$ into equation (13)– the emission rate becomes:

$$\begin{equation} r^{\rm em}_{\rm H}\! =\! \left(\!\!\frac{2\pi b\rho_{\rm d}\Omega}{a_0^2}\!\!\right)\nu_0\exp\left(\!\!-\frac{E_d}{kT}\!\!\right)\left(\!N^d_{\rm H}-\rho_a\Omega\exp\left(\!\!-\frac{\Delta E^f_{\rm H}}{kT}\!\!\right)\!\!\right) \end{equation} \tag{ 14 }$$

For their part, the dissociation coefficients of V_mH_n clusters are calculated assuming thermal emission of monomers (H or V) as:

$$\begin{equation} r^{\rm diss}_i = \frac{b}{b+r_i} \frac{4\pi r_i^2}{a_0^2} {\nu}_0 \exp\left(-\frac{E^b_{\rm \{H,V\}}+E^m_{\rm \{H,V\}}}{kT}\right)N_{i} \end{equation} \tag{ 15 }$$

where $E^b_{\rm \{H,V\}}$ are appropriate binding energies of monovacancies or H atoms to a cluster i, and $E^m_{\rm \{H,V\}}$ are the corresponding migration barriers of vacancies or H atoms in the bulk. r_i represents the size of the cluster, typically calculated by assuming a spherical vacancy cluster shape:

$$\begin{equation*} r_i = \left( \frac{m3\Omega_v}{4\pi}\right)^{1/3} \end{equation*}$$

with m the number of vacancies in the V_mH_n cluster and Ω_v the vacancy formation volume.

Broadly speaking, the super abundant vacancy (SAV) mechanism is a process by which hydrogen assists in the generation of vacancies in a material. Numerous examples of the SAV effect exist in metals, e.g. Ni, Cr, Pd, Al, Mo, or Nb [72–75]. At the atomistic level, the fundamental idea is that as vacancies and vacancy clusters capture hydrogen atoms their internal pressure grows up to an unstable point, after which a vacancy-SIA pair is produced with the extra vacancy helping to relieve the pressure (i.e. decrease the H/V ratio) of the cluster, and a SIA emitted into the bulk. In such fashion, the system is in principle able to absorb an indefinite amount of hydrogen, provided that the H/V ratio (n/m) is kept within the stability limits of the equation of state of the clusters. This is akin to the well-known 'trap-mutation' mechanism observed in W exposed to He plasmas [76]. For hydrogen in tungsten, several authors have introduced this mechanism and discussed the energetics of the process to determine the vacancy-to-hydrogen ratios that trigger the reaction [77–79]. When the energetics of the reaction are favorable, the sequence goes as follows:

$$\begin{equation} {\rm V}_m{\rm H}_n + {\rm H}_1 \to {\rm V}_m{\rm H}_{n+1} \to {\rm V}_{m+1}{\rm H}_{n+1} + {\rm I}_1 \end{equation} \tag{ 16 }$$

such that the ratio changes from $\frac{n+1}{m}$ to $\frac{n+1}{m+1}$ and a SIA is inserted into the lattice. The stability of the reaction can be assessed in terms of the energies of each of the reactants and products:

$$\begin{equation} E^f\left({\rm V}_{m+1}{\rm H}_{n+1}\right) + E^f({\rm I}_1)\leq E^f\left({\rm V}_{m}{\rm H}_{n+1}\right) \end{equation} \tag{ 17 }$$

such that reaction (16) will take place when the excess energy

$$\begin{equation} \Delta E_{\rm SAV} = E^f\left({\rm V}_{m+1}{\rm H}_{n+1}\right) + E^f({\rm I}_1)-E^f\left({\rm V}_{m}{\rm H}_{n+1}\right) \end{equation} \tag{ 18 }$$

becomes negative. These formation energies can be calculated as:

$$\begin{align*} E^f\left({\rm V}_{m}{\rm H}_{n}\right) & = nE^f({\rm H}_1)+E^f({\rm V}_m)-\left[E^b({\rm V}_m\text{-}{\rm H}_n)+m\mu_{\rm W}\right]\\ E^f({\rm V}_m) & = E^f({\rm V}_1) +E^f({\rm V}_{m-1}) - E^b({\rm V}_{m-1}\text{-}{\rm V}_1) \end{align*}$$

where $E^b({\rm V}_m\text{-}{\rm H}_n)$ is the binding energy between a pure H aggregate with n atoms and a V_m vacancy cluster, and $\mu_{\rm W}$ is the cohesive energy (eV/atom) of W. $E^f({\rm H}_1)$, $E^f({\rm I}_1)$, and $E^f({\rm V}_1)$ are the heat of solution of H in W, and the vacancy and self-interstitial atom formation energies, respectively. Values for the required binding energies are given in appendix A (tables A1 and A3), while the rest of the parameters are: $\mu_{\rm W} = 8.90$ eV [80], $E^f({\rm H}_1) = 1.04$ eV [81], $E^f({\rm I}_1) = 9.96$ eV [82], and $E^f({\rm V}_1) = 3.23$ eV [83].

Figure 4 shows a color map of $\Delta E_{\rm SAV}$. The white solid line marks the limit of stability of the SAV mechanism. Above it, reaction (16) becomes energetically favorable and the SAV mechanism becomes active, while, below it, V_mH_n clusters are stable. This criterion guides the behavior of the clusters and is an independent source of growth during the simulations. This solid line can be approximated by a linear relation whose slope yields the critical x = n/m ratio. For our model, this is equal to x = 4.0, represented as a white dashed line in the figure. However, while figure 4 indicates the thermodynamic stability region of the SAV mechanism, its occurrence may also be controlled by a kinetic barrier (in the manner of standard chemical reactions). Such an energy, which we term $E^{\rm SAV}_0$, defines a transition rate that is enabled only when $\Delta E_{\rm SAV} \lt 0$:

$$\begin{equation} r_{\rm SAV} = \nu_0 \exp\left(-\frac{E^{\rm SAV}_0}{kT}\right) \end{equation} \tag{ 19 }$$

At present, there are no experimental or numerical estimates of $E^{\rm SAV}_0$, and it thus must be obtained through other means. We will return to this point when we discuss the model results below.

Zoom In Zoom Out Reset image size

Table 2 gives all the numerical parameters and simulation conditions in accordance with the three experimental phases discussed in section 2.

Table 2. Material and simulation parameters employed in the simulations.

W material parameters:	Symbol	Value	Units
Atomic density	ρa	6.31 × 1028	$[\mathrm{m^{-3}}]$
Lattice parameter	a0	3.16	Å
Dislocation density	ρd	1010	$[\mathrm{m^{-2}}]$
Cu irradiation parameters:
Irradiation temperatures	$T_{\rm irr}$	300,	573,	873, 1023,	1243	[K]
Ion flux	φ	1.34 × 1015	$[\mathrm{m^{-2}s^{-1}}]$
Ion energy	$E_{\rm ion}$	3.4	[MeV]
Dose	(φt)	0.2	[dpa]
H plasma exposure parameters:
Temperature	-	383	[K]
Flux	$\mathrm{\phi_H}$	4.0 × 1020	$[\mathrm{m^{-2}s^{-1}}]$
Time duration	Δt	2500	s
Deposition energy	$E_{\rm H}$	113	$[\mathrm{eV}]$
Thermal desorption parameters:
Heating rate	β	0.5	[K· $\mathrm{s^{-1}}$]
Temperature range	ΔT	300 ∼ 1300	[K]
Numerical parameters:
Number of spatial elements	n	101	-
Element volume	Ωi	10−23	[m3]
Top layer thickness	ls	0.54	[nm]
Element thickness	li	20	[nm]

The calculation of the transport energetics of defects and hydrogen atoms in tungsten has attracted a great deal of attention in recent years [43, 80, 84, 84–104]. Based on an extensive literature review, we have selected a set of parameters for the diffusivities and energetics of defects and hydrogen, summarized in tables 3 and 4. The first table gives the values of the diffusivities of all mobile species assuming the standard Arrhenius expression for thermally-activated motion $D(T) = D_0\exp\left(-\frac{E^m}{kT}\right)$. Table 4 gives additional values for the parameters used in figure 1.

Table 3. Diffusion coefficients of the mobile species considered in this work.

Species	D0 [m2 s−1]	Em [eV]	Source
SIA clusters, Iq:
$\mathrm{I_1}$	$\mathrm{8.74\times 10^{-8}}$	0.009	[103]
$\mathrm{I_2}$	$\mathrm{7.97 \times 10^{-8}}$	0.024	[103]
$\mathrm{I_3}$	$\mathrm{3.92 \times 10^{-8}}$	0.033	[103]
q  > 3	$2.99\times 10^{-7}q^{-0.5}$	0.013	[104]
Vacancy clusters, Vn:
$\mathrm{V_1}$	$\mathrm{177\times 10^{-6}}$	1.29	[103]
$\mathrm{V_2}$	$\mathrm{2.91\times 10^{-9}}$	1.66	[104]
m > 2	$4.01\times 10^{-(5+3\ {\rm m})}$	1.66	[104]
Hydrogen:
$\mathrm{H_1,~H_2}$	$\mathrm{1.58\times10^{-7}}$	0.25	[43]

Table 4. H-atom energetics in W (refer to figure 1 for details). All energies are given in eV (from reference [42]).

Symbol	Value [eV]
$E^{\rm diss}$	0.0
$E^{\rm abs}$	1.44
Er	0.25
Es	0.79

Cu-ion irradiation, described in section 2, to a total dose of 0.2 dpa at 300, 573, 873, 1023, and 1243 K results in the accumulation of point defect and defect clusters in the irradiated specimens. Figure 5(a) shows the concentration of vacancy-type defects⁸ in the entire specimen as a function of dose calculated with SR-SCD. As the figure shows, the defect concentration experiences an incubation period of rapid increase characterized by a power law with scaling exponent 0.5 (which appears to be independent of temperature), followed by convergence to steady state. This steady state is reached more rapidly at elevated temperature, and its final value is indicative of the point at which vacancy and vacancy cluster production is balanced by absorption at sinks. At 300 K the system has not reached this steady state after 0.2 dpa. The accumulation of vacancy defects is substantial, with concentrations approaching 0.02 at.% at 300 and 573 K. The steady state concentrations roughly decrease linearly with temperature. The buildup of SIA clusters is seen to follow identical qualitative trends as vacancies, although with much lower concentrations.

Zoom In Zoom Out Reset image size

The depth distributions at 0.2 dpa for each temperature are plotted in figure 5(b). Consistent with figure 3, defects are distributed uniformly into the material to a depth of approximately 1.4 microns. These initial vacancies and vacancy clusters become the seeds for the formation of larger hydrogen-trapping V_mH_n complexes during the plasma exposure phase. The clusters are generally small. Figure 6 shows the cluster size distribution after the irradiation phase, with the largest one being a V₆ at 573 K.

Zoom In Zoom Out Reset image size

Next, we simulate the next two phases, i.e. hydrogen plasma exposure and thermal desorption of all irradiated samples. First, however, as introduced in section 3.4, a key parameter to be defined before entering a complete set of simulations of the H exposure phase is the energy barrier $E_0^{\rm SAV}$. In the absence of direct measurements or calculations of $E_0^{\rm SAV}$, here we take a backend approach to ascertain its value: we systematically vary it from zero to 1.5 eV and simulate plasma exposure and thermal desorption for each value. We then compare the simulated TDS to the experimental ones and select the value of $E_0^{\rm SAV}$ that produces the best match as evaluated using a least-squares error analysis. As an advance of more results to come, in figure 7 we show results for six different scenarios: 0 (equivalent to a 'spontaneous' SAV mechanism), 0.90, 0.95, 1.00, and 1.20 eV. We also show results for simulations under no SAV mechanism. Our results suggest that the best match is obtained for $E_0^{\rm SAV} = 0.95$ eV, which is the value used hereafter to showcase the hydrogen exposure and thermal desorption phases.

Zoom In Zoom Out Reset image size

In accordance with experimental conditions, the hydrogen exposure stage is simulated at a temperature of 383 K and a duration 2500 s under a constant hydrogen flux of 4.0 × 10²⁰ m⁻²s⁻¹ (total fluence 10²⁴ m⁻²). Figure 8 shows the distribution of hydrogen as a function of depth for each irradiation temperature at the end of the exposure period. Note that the concentrations shown include hydrogen in any form, i.e. free hydrogen monomers as well as H atoms trapped in vacancy clusters. Irradiation temperature is seen to have a noticeable impact on the depth distribution of hydrogen, largely correlated with the distribution of vacancy clusters during the irradiation phase. In general, hydrogen is found deeper, albeit in lower concentrations, the larger $T_{\rm irr}$ is.

Zoom In Zoom Out Reset image size

The analysis of the V_mH_n populations displayed in figure 8 is critical, as it quantitatively establishes the trapping propensity of the material under the conditions explored here. There are two main aspects of this population that must be further studied. One is the relative partition of hydrogen-to-vacancy ratios, x = n/m within the cluster subpopulation. The other is the absolute size of the clusters, which has implications for both the total amount of H stored in the system and the visibility under the microscope of the defects. Next, we provide a detailed analysis of both of these aspects in an attempt to ascertain the fuel footprint in irradiated single-crystal W surfaces.

5.3.1. Stability of hydrogen-vacancy clusters and distribution of hydrogen-to-vacancy ratios.

Figure 9 shows the integrated concentration of V_mH_n clusters as a function of x for each irradiation temperature case. In keeping with the limits established in figure 4, a dashed line is drawn at x = 4 to show the relative stability of clusters based on their hydrogen-to-vacancy ratio. Clusters with $x \lt 4$ are generically deemed to be 'underpressurized', i.e. with a propensity to absorb more hydrogen atoms, while clusters with x > 4 are said to be 'overpressurized', i.e. with propensity to release hydrogen (or, when the SAV mechanism is active, to produce vacancy-SIA pairs). On the basis of this partition, two things can be established about the different $T_{\rm irr}$ cases: (i) a higher irradiation temperature leads to preferentially underpressurized clusters, while lower temperatures result in an equipartition between under- and overpressurized bubbles; (ii) the range of x observed is 0 to 9.

Zoom In Zoom Out Reset image size

Figure 10 shows the depth distribution of the V_mH_n clusters in terms of their x value as a function of $T_{\rm irr}$ at the end of the H-exposure phase. The colored bands around each curve represent error bars from numerical fluctuations over five independent SR-SCD runs. Generally, little correlation is observed between x and the depth at which each cluster type is found. Consistent with figure 8, the main effect of irradiation temperature is to result into deeper V_mH_n cluster penetration.

Zoom In Zoom Out Reset image size

Finally, it is useful to represent the cluster population not just in terms of their x ratios but also in terms of their absolute m and n numbers. Figure 11 shows three-dimensional plots of V_mH_n clusters as a function of m, n, and the depth at which they are found. Their concentration is given using the color key as specified in the figure. Each subplot corresponds to a specific irradiation temperature. In a similar manner to the distribution of x ratios, key qualitative differences can be appreciated between the lower and higher temperature cases. At lower temperatures, figures 11(a)–(c), most V_mH_n clusters are located within the first micron of depth, with m and n covering a large span between V$_{1\sim 6}$H$_{10\sim 30}$ to approximately V$_{7500\sim 8000}$H$_{2200\sim 2350}$ and relatively high concentrations. By contrast, at 1023 and 1243 K (figure 11(d), the figure for $T_{\rm irr} = 1243$ K is not shown due to its qualitative and quantitative similarity with the 1023-K case) only clusters composed of more than 5500 vacancies are observed at high depths. These large clusters appear in relatively low concentrations. We can synthesize the information contained in these figures into a global size distribution for the end of the exposure phase. This represent essential information to facilitate comparison with experimental distributions, which are limited on the low end by machine resolution. Our results are shown in figure 12 in the form of temperature-smeared Gaussian distributions, with the raw data shown as a histogram in the background. While the figure clearly shows a bi-modal distribution with a peak within one nanometer sizes and another one centered around 6 nm, the standard resolution limit in microscopy of 0.5-to-1.0 nm may prevent revealing the former in experimental size distributions. Techniques such as x-ray diffuse scattering can achieve sub-nanometer resolutions, thus enabling a better comparison with our results [39]. We discuss this further in section 6.2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Thermal desorption processes are thought to be composed of a series of emission events characterized by Gaussian profiles [105, 106]. These Gaussians are centered at temperatures defined by characteristic binding energies representing specific hydrogen release events. It is therefore useful to decompose the desorption profiles shown in the figure into a set of overlapping Gaussian distributions that coincide with the most dominant peaks (or protuberances) in each TD spectrum. Doing this, however, may involve a certain degree of arbitrariness, as these peaks are not always evident (particularly at high temperatures). In any case, we believe that some clear mappings exist and some consistency among the different curves can be found.

Here we focus on the 0.95-eV case in figure 7. Figure 13 shows the decomposed thermal desorption spectra for each irradiation temperature. In performing each decomposition, Gaussians with like colors are used when they are centered at the same desorption temperature (horizontal axis), regardless of the irradiation temperature case. The width of the Gaussians is temperature dependent and has been adjusted so that the contribution to the TD spectrum from the overlap among them is consistent with the full curve. As such, the 300-K and 573-K temperature case (figure 13(a) and 13(b)) is composed of seven distinct Gaussians, each representing a specific desorption mechanism, centered at, respectively, 333, 425, 485, 527, 574, 623, and 836 K. Two of these are shared with the 873-K case, figure 13(c), which is suggestive of identical H-release mechanisms in these two cases. However, the remaining spectra corresponding to the higher irradiation temperatures are characterized by different peaks, namely 310, 842 K (figure 13(d)).

Zoom In Zoom Out Reset image size

Table 5. Operative dissociation mechanisms at each thermal desorption peak temperature for each irradiation temperature. The symbol '⊥' represents 'dislocations'. The simulated dissociation energy, $E^{d}_{\rm sim}$ is obtained as $E^{d}_{\rm sim} = E^b+E^m$. The values for E^b are extracted from table A3 for each specific reaction, while E^m comes from table 3. The values for $E^d_{\rm TDS}$ were obtained from equation (20).

T [K]	Irradiation temperature [K]					$E^{d}_{\rm sim}$ [eV]	$E^{d}_{\rm TDS}$ [eV]
	300	573	873	1023	1243
310	–	–	–	⊥	⊥	0.85	0.70
333	⊥	⊥	⊥	–	–	0.85	0.75
425	V1H8 →V1H7	V1H8 → V1H7	–	–	–	0.88	0.97
485	V1H7 → V1H6	V1H7 → V1H6	–	–	–	1.05	1.11
	V6H23 → V6H22		–	–	–	1.06
527	V1H6 → V1H5	V1H6 → V1H5	–	–	–	1.18	1.21
	V6H21 → V6H20	V6H22 → V6H21	–	–	–	1.18
574	V1H5 → V1H3	V1H5 → V1H3	–	–	–	1.28	1.32
	V6H18 → V6H17	–	–	–	–	1.34
623	V1H3 → V1	V1H3 → V1	–	–	–	1.49	1.44
	V6H16 → V6H12	V6H16 → V6H12	–	–	–	1.42
836	V$_{202\sim 818}$H$_{1\sim 67}$ → V$_{202\sim 818}$H$_{0\sim 66}$	V$_{172\sim 3593}$H$_{1\sim 181}$ → V$_{172\sim 3593}$H$_{0\sim 181}$	V$_{145\sim 2100}$H$_{1\sim 109}$ → V$_{145\sim 2100}$H$_{0\sim 108}$	–	–	1.69∼ 1.96	1.96
842	–	–	–	V$_{563\sim 840}$H$_{104\sim 171}$ → V$_{563\sim 840}$H$_{103\sim 170}$	V$_{692\sim 844}$H$_{1\sim 51}$ → V$_{692\sim 844}$H$_{0\sim 50}$	1.94∼ 1.96	1.97

As mentioned above, each one of the Gaussians shown in figure 13 represents a distinct hydrogen emission mechanism. These mechanisms can be directly identified in the simulations from processes defined by equations (14) and (15)), so that we can ascribe a specific dissociation reaction to each emission temperature, T_p. Further, we characterize each Gaussian by a representative dissociation energy $E^d_{\rm TDS}$, determined as [107]:

$$\begin{equation} \frac{E^{d}_{\rm TDS}}{kT_p^2} = \frac{\nu_0}{\beta}\exp\left(-\frac{E^d_{\rm TDS}}{kT_p}\right) \end{equation} \tag{ 20 }$$

where β is the heating rate. Each $E^d_{\rm TDS}$ can subsequently be compared with the corresponding dissociation energies corresponding to the specific reactions taking place in the simulations. Those ultimate originate from the data given in table A3. The results of such analyses are provided in table 5.

The information contained in the table is the culmination of the simulation effort of the three-stage experimental study described in section 2. The first column represents all the Gaussian peaks used in figure 13. Columns 2-5 specify the transitions suffered by the clusters during each emission process. With increasing desorption temperature, the emission events include hydrogen emission from dislocations, small overpressurized clusters gradually releasing hydrogen to reduce their hydrogen-to-vacancy ratio towards x = 4, and large underpressurized clusters emitting H atoms until all the hydrogen is evacuated from the system. Finally, the last column lists the corresponding dissociation energies of each dissociation event as calculated with equation (20). We find that, in general, the discrepancy between the $E^d_{\rm TDS}$ given in table 5 and the actual dissociation energies from table A3 is small. While this simply corroborates that the reactions identified in each temperature interval are consistent with the Gaussian emission peaks, it is an encouraging result that adds confidence to our results.

Hydrogen thermal desorption in damaged crystals is an extremely complex problem involving numerous processes often co-occurring in a synergistic or non-linear manner. Just listing those considered here gives an idea of this complexity: cascade damage, defect diffusion, hydrogen deposition and penetration, hydrogen passivation, internal hydrogen diffusion, hydrogen absorption and trapping at dislocations and grain boundaries, vacancy cluster-hydrogen reaction, production of vacancies from cluster internal overpressurization, and hydrogen dissociation from clusters, all of it captured with spatiotemporal resolution. The amount of physics needed to undertake simulations of such a challenging process can be staggering, both in terms of the kinetics and thermodynamics of the atomic species involved, as well as in terms of the detailed parameterization needed for physical accuracy and validation. As well, carefully-conducted experiments covering the above processes and in such a way as to facilitate comparison with the models are an invaluable tool to validate and refine the simulations. Fortunately, we believe that the fusion materials community has reached sufficient maturity in all of these fronts to enable the type of study presented here. For example, a systematic effort by researchers to map the energetics of the W-H system using electronic structure calculations and semi-empirical potentials, there are now very reliable data sets available. Thanks to advances in algorithms and computer power, there now exist powerful reaction-diffusion PDE models to study strongly-coupled multispecies problems in great detail. This, combined with several decades of 'know-how' in irradiation damage modeling has open the door to simulations such as those presented here which allow us to push the physical accuracy/computational efficiency tradeoff to unprecedented levels.

In undertaking the present simulation effort, we must also be transparent about two key premises adopted here, namely the presumed infallibility of both experimental results and atomistic calculations. It is beyond the scope of this work to comment in depth about the validity of such premises, and the reader is simply referred to the pertinent sources for details (e.g. references [33, 108] with regard to atomistic calculations, references [44, 58] for the irradiation conditions and plasma exposure, and references [109, 110] for uncertainties in thermal desorption spectroscopy) for details. In general, it is probably safe to say that, as has become common practice in the computational materials science community, both experimental characterization and ab initio calculations have reached a degree of reliability as it relates to PMI studies that allow us to confidently use them in our studies as the 'true' baseline against which to compare model predictions.

With this in mind, the main result of this work is presented in figure 14 (subpanel in figure 7), which we use to discuss the main physical insights gained from the current exercise. Analysis of figures such as those points to two principal features that stand out from within the body of results obtained here:

• (a)  
  Our results conclusively show that, vis-à-vis the experiments, the superabundant vacancy mechanism is a key piece of physics needed to understand the accumulation of vacancy-hydrogen clusters in damaged W subjected to hydrogen exposure point to the existence of a transition barrier. Moreover, we find that the SAV mechanism is accompanied by a transition barrier that modulates changes in the cluster structure. Mapping the simulated TDS to the experimental ones points to a value of 0.95 eV for this transition energy.
• (b)  
  The thermal desorption spectra of H-exposed damaged W surfaces can be broadly decomposed into three distinct emission temperature regions:These temperature regimes are schematically depicted in figure 15, which are in strong qualitatively agreement with the experimental TDS obtained in reference [58].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this context, the SAV mechanism plays a fundamental role in the final accumulation and release of hydrogen from damaged tungsten surfaces. Going back to figure 7, the absence of a SAV-type mechanism overemphasizes smaller, overpressurized cluster populations, while an unabated⁹ SAV mechanism leads to an abundance of larger, underpressurized bubbles. In view of these trends, the effect of the SAV activation barrier is clear. $E_0^{\rm SAV}$ modulates the two extremes (no SAV, spontaneous SAV) and bridges the temperature regions where overpressurized and underpressurized clusters dominate. This is clearly seen in figure 7 and in the experiments, where for $E_0^{\rm SAV} = 0.95$ eV our results show remarkable agreement with both the emission peaks and the emission fluxes measured experimentally.

A first attempt at validation of our results can be made by comparing the Sun et al [39] (figure 6 in their paper). At 300 K and 0.2 dpa with 5-MeV Cu ions, the measured vacancy defect concentration was 0.103± 0.02 at. %, while our integrated value from figure 5(b) is 0.023± 0.004 at. %. Part of this 5 × difference is likely attributable to the difference in irradiation energy (5.0 vs 3.4 MeV), although we do not discount other factors related to more fundamental issues such as the model framework. Noteworthy is also the fact that the average cluster sizes in each case (measured vs. simulated) are 10 ±  3 Å and 8.32 ±  0.3 Å, in very close agreement with one another.

Second, related to the experimental effort discussed in section 2, the desorption hydrogen fluence can be measured for each of the irradiation-temperature specimens and compared to the integrated hydrogen release in the simulations. The measurements are shown in figure 16 for two Cu-ion irradiation energies together with the present results. A difference factor of about 2.5 × is observed, which can be partially explained from the analysis of the difference in TDS features that we discuss next.

Zoom In Zoom Out Reset image size

Analyses of the experimental peaks between 400 and 700 K (intermediate temperature range) suggests that our simulations underestimate the concentrations of clusters that give rise to the first set of emission peaks (between 400 and 550 K), while they slightly overestimate those in the range between 600 and 700 K. We can specifically point to the reactions that govern each subregime from among those listed in table 5 to analyze this discrepancy. Between 400 and 550 K the peaks are due to the decomposition of V₁H₇, V₁H₈, V₆H₂₁, and V₆H₂₂ clusters. The implication would thus be that the simulations are underestimating these populations relative to the experimental observations. Conversely, between 600 and 700 K the relevant clusters are V₁H₃, V₁H₅, V₆H₁₆, and V₆H₁₈, whose populations would appear to be overestimated in the simulations. This could be due to a number of causes, although given how similar the clusters are in both subregimes, subtle changes in the binding energies would likely suffice to correct this. This is one way simulations such as these can inform the atomistic calculations of the energetics of the W-H system.

Typically, H emission peaks below room temperature are not captured in tungsten in standard TDS experiments. In the few instances where the starting temperature is 300 K or below, experiments have indeed revealed emission of hydrogen and deuterium [111–113]. This is consistent with electronic structure calculations of H trapping energies at edge and screw dislocations ranging between 0.5 and 0.9 eV [70, 114], which would result in desorption peaks below 400 K. Our model includes this latest computational input and, accordingly, it produces a peak at around 300 K. It is worth adding that grain boundaries are predicted to trap hydrogen with energies similar to those of dislocations [ZHOU2011240], such that a contribution to fuel retention from GBs or other two-dimensional defects, if they exist, should not be discounted.

In any case, only dedicated experiments designed for TDS below room temperature can conclusively determine whether the peak exists or not. Results such as those shown here can provide an opportunity for experimental science to verify such low-temperature behavior.

The standard treatment of SIA-type defects in crystalline metals is to assume rapid one-dimensional migration away from the core of displacement cascades towards sinks. When cascades are sufficiently dense from high-energy PKA, dislocation loops can be formed, which, after a critical fluence, may lead to so-called 'rafts' or loop networks. These networks have been shown to be effective hydrogen traps when exposed to a H atmosphere [115].

In the present irradiations, the combination of superficial damage production (the peak dpa occurs at a depth of 0.7 microns) and relatively low PKA energy (≈ 2 keV) results in little SIA-defect accumulation. Additionally, recent DFT calculations have found trapping energies for small H-SIA complexes of around 0.3 eV in W [116, 117]. These energies result in H release in thermal desorption tests at temperatures well below room temperature. For these reasons, we are confident that the SIA-defect contribution to H trapping in the conditions studied in this work can be neglected without much effect on the final results.

When SIA clusters grow to loop sizes, they can become indistinguishable from the dislocation network and trap hydrogen with 0.5∼ 0.9 eV. This of course is already captured in our model, except that the dislocation density is time-invariant. Instances of higher energy PKA over long doses would probably warrant the inclusion of a dislocation density evolution model coupled to the cluster dynamics model to account for loop network formation.