Mechanism of Re precipitation in irradiated W-Re alloys from kinetic Monte Carlo simulations

TL;DR

This study uses kinetic Monte Carlo simulations to uncover the mechanism of Re cluster formation in irradiated tungsten alloys, revealing how mixed-interstitial transport influences nucleation and growth of Re-rich clusters, which impact material properties.

Contribution

It introduces a mechanistic model incorporating mixed-interstitial transport for Re clustering in irradiated W-Re alloys, validated by simulations aligned with experimental data.

Findings

Re clusters form after ~20 s at 1800 K

Clusters grow to about 4 nm radius with 50% Re

Simulation results agree with atom probe experiments

Abstract

High-temperature, high-dose, neutron irradiation of W results in the formation of Re-rich clusters at concentrations one order of magnitude lower than the thermodynamic solubility limit. These clusters may eventually transform into brittle W-Re intermetallic phases, which can lead to high levels of hardening and thermal conductivity losses. Standard theories of radiation enhanced diffusion and precipitation cannot explain the formation of these precipitates and so understanding the mechanism by which nonequilibrium clusters form under irradiation is crucial to predict materials degradation and devise mitigation strategies. Here we carry out a thermodynamic study of W-Re alloys and conduct kinetic Monte Carlo simulations of Re cluster formation in irradiated W-2Re alloys using a generalized Hamiltonian for crystals containing point defects parameterized entirely with electronic structure calculations. Our model incorporates recently-gained mechanistic information of mixed-interstitial solute transport, which is seen to control cluster nucleation and growth by forming quasi-spherical nuclei after an average incubation time of 20 s at 1800 K. These nuclei are seen to grow by attracting more mixed interstitials bringing solute atoms, which in turns attracts vacancies leading to recombination and solute agglomeration. The clusters grow to a maximum size of approximately 4-nm radius, and are not fully dense with Re, containing 50% or less near the center. Our simulations are in reasonable agreement with recent atom probe examinations of ion irradiated W-2Re systems at 773 K.

Tables

Table 1. Table 1: Energetics of W-Re systems calculated with DFT.

Quantity	Value	Source
$E_{coh}^{\text{A}}$	$8.3276$	This work
$E_{coh}^{\text{B}}$	$7.4070$	This work
$\Delta H^{\mathrm{mix}}$	$-0.1571-0.2311X$	Ref. Gharaee et al. (2016)
$E_f^{\text{V}}$	$3.1690$	This work
$E_b^{\left(a\right)}$	$-0.2096$	This work
$E_b^{\left(b\right)}$	$-0.1520$	This work
$E_b^{\left(c\right)}$	$-0.3079$	This work
$E_b^{\left(d\right)}$	$-0.2992$	This work
$E_{b,\text{1nn}}^{\text{V-V}}$	-0.0146	This work${}^{\left(\text{<a><span>3</span></a>}\right)}$
$E_{b,\text{2nn}}^{\text{V-V}}$	0.3028	This work${}^{\left(\text{<a><span>3</span></a>}\right)}$
$E_f^{\text{AA}}$	$10.16$	Ref. Gharaee and Erhart (2015)
$E_f^{\text{AB}}$	$9.49$	Ref. Gharaee and Erhart (2015)
$E_{b,\text{1nn}}^{\text{AA-B}}$	$-0.52$	Ref. Gharaee and Erhart (2015)
$E_{b,\text{1nn}}^{\text{AB-B}}$	$-0.53$	Ref. Setyawan et al. (2016)
$E_{b,\text{1nn}}^{\text{AA-AA}}$	$-2.12$	Ref. Becquart et al. (2010)
$E_{b,\text{1nn}}^{\text{AA-AB}}$	$-2.12$	Assumed ${}^{\left(\text{<a><span>4</span></a>}\right)}$
$E_{b,\text{1nn}}^{\text{AB-AB}}$	$-3.2$	Ref. Gharaee et al. (2016)
$E_m^{\text{V}\to \text{A}}$ (A)	$1.623$	This work
$E_m^{\text{V}\to \text{B}}$ (A)	$1.651$	This work
$E_m^{\text{V}\to \text{A(1)}}$ (Fig. 2(c))	$1.7151$	This work
$E_m^{\text{V}\to \text{A(2)}}$ (Fig. 2(c))	$1.6378$	This work
$E_m^{\text{V}\to \text{B(3)}}$ (Fig. 2(c))	$1.577$	This work
$E_m^{\text{V}\to \text{A}}$ (V)	$1.623$	This work
$E_m^{\text{V}\to \text{B}}$ (V)	$1.651$	This work

Table 2. Table 2: Bond energy coefficients with the equation used for their calculation, and the literature source. x 𝑥 x is the local solute concentration

${ϵ}_{\text{A-A}}^{\left(1\right)}$	$-1.5815$	cohesive energy, eq. (9)	This work
${ϵ}_{\text{A-A}}^{\left(2\right)}$	$-0.6672$	cohesive energy, eq. (9)	This work
${ϵ}_{\text{B-B}}^{\left(1\right)}$	$-1.4067$	cohesive energy, eq. (9)	This work
${ϵ}_{\text{B-B}}^{\left(2\right)}$	$-0.5935$	cohesive energy, eq. (9)	This work
${ϵ}_{\text{A-B}}^{\left(1\right)}$	$-1.5090-0.0219x$	mixing energy	Ref. Gharaee et al. (2016)
${ϵ}_{\text{A-B}}^{\left(2\right)}$	$-0.6366-0.0092x$	eq. (9)	Ref. Gharaee et al. (2016)
${ϵ}_{\text{A-V}}^{\left(1\right)}$	$-0.4898$	formation energy, eq. (9)	This work
${ϵ}_{\text{A-V}}^{\left(2\right)}$	$-0.2067$	formation energy, eq. (9)	This work
${ϵ}_{\text{B-V}}^{\left(1\right)}$	$-0.4898-0.009432/x$	formation energy fitted to ${ϵ}_{B-V}^{\left(1\right)}=a+b/x$	This work
${ϵ}_{\text{B-V}}^{\left(2\right)}$	$-0.3311+0.036x$	formation energy fitted to ${ϵ}_{B-V}^{\left(1\right)}=a+bx$	This work
${ϵ}_{\text{V-V}}^{\left(1\right)}$	$0.5873$	1nn binding energy	This work
${ϵ}_{\text{V-V}}^{\left(2\right)}$	$0.5566$	2nn binding energy	This work
${ϵ}_{\text{AA-A}}^{\left(1\right)}$	$0.1740$	formation energy, eq. (9)	Ref. Gharaee and Erhart (2015)
${ϵ}_{\text{AA-A}}^{\left(2\right)}$	$0.0734$	formation energy, eq. (9)	Ref. Gharaee and Erhart (2015)
${ϵ}_{\text{AB-A}}^{\left(1\right)}$	$0.1104$	formation energy, eq. (9)	Ref. Gharaee and Erhart (2015)
${ϵ}_{\text{AB-A}}^{\left(2\right)}$	$0.0466$	formation energy, eq. (9)	Ref. Gharaee and Erhart (2015)
${ϵ}_{\text{AA-B}}^{\left(1\right)}$	$-0.2750$	binding energy	Ref. Gharaee and Erhart (2015)
${ϵ}_{\text{AA-B}}^{\left(2\right)}$	$-0.1160$	eq. (9)	Ref. Gharaee and Erhart (2015)
${ϵ}_{\text{AB-B}}^{\left(1\right)}$	$-0.3486$	binding energy	Ref. Setyawan et al. (2016)
${ϵ}_{\text{AB-B}}^{\left(2\right)}$	$-0.1470$	eq. (9)	Ref. Setyawan et al. (2016)
${ϵ}_{\text{AA-AA}}^{\left(1\right)}$	$-0.1905$	binding energy	Ref. Becquart et al. (2010)
${ϵ}_{\text{AA-AA}}^{\left(2\right)}$	$-0.0804$	eq. (9)	Ref. Becquart et al. (2010)
${ϵ}_{\text{AA-AB}}^{\left(1\right)}$	$-0.2505$	binding energy	Assumed ${}^{\left(\text{<a><span>4</span></a>}\right)}$
${ϵ}_{\text{AA-AB}}^{\left(2\right)}$	$-0.1057$	eq. (9)	Assumed ${}^{\left(\text{<a><span>4</span></a>}\right)}$
${ϵ}_{\text{AB-AB}}^{\left(1\right)}$	$-1.3977$	binding energy	Ref. Gharaee et al. (2016)
${ϵ}_{\text{AB-AB}}^{\left(2\right)}$	$-0.5897$	eq. (9)	Ref. Gharaee et al. (2016)

Table 3. Table 3: Saddle-point bond energy coefficients for vacancy jumps (in eV).

${ϵ}_{\text{A-A}}^{sp\left(1\right)}$	$-2.5975$	${ϵ}_{\text{A-A}}^{sp\left(2\right)}$	$-0.5041$
${ϵ}_{\text{A-B}}^{sp\left(1\right)}$	$-2.6451$	${ϵ}_{\text{A-B}}^{sp\left(2\right)}$	$-0.5532$
${ϵ}_{\text{A-V}}^{sp\left(1\right)}$	$0.5465$	${ϵ}_{\text{A-V}}^{sp\left(2\right)}$	$0.1060$
${ϵ}_{\text{B-A}}^{sp\left(1\right)}$	$-2.5188$	${ϵ}_{\text{B-A}}^{sp\left(2\right)}$	$-0.4888$
${ϵ}_{\text{B-B}}^{sp\left(1\right)}$	$-2.5417$	${ϵ}_{\text{B-B}}^{sp\left(2\right)}$	$-0.4943$
${ϵ}_{\text{B-V}}^{sp\left(1\right)}$	$0.2902$	${ϵ}_{\text{B-V}}^{sp\left(2\right)}$	$0.0563$

Table 4. Table 4: Self-interstitial migration parameters. The jump distance for SIA migration is δ = a 0 ​ 3 / 2 𝛿 subscript 𝑎 0 3 2 \delta=a_{0}\sqrt{3}/2 .

$E_m^{\text{AA}}$	0.003	Ref. Suzudo et al. (2014)
$E_r^{\text{AA}}$	0.43	Ref. Suzudo et al. (2014)
$E_m^{\text{AB}}$	0.12	Ref. Gharaee et al. (2016)
$f$	$2.93-0.00055T$	Ref. Zhou et al. (2013)

Table 5. Table 5: Diffusion parameters for vacancy and solute diffusion as a function of solute concentration.

$X$ [at. %]	$D_0$ [m2$\cdot$s-1]	$E_m$ [eV]
Vacancy diffusion
0.0 (eq. (29))	$4.84\times {10}^{-7}$	1.62
0.5	$6.86\times {10}^{-6}$	1.73
1.0	$6.92\times {10}^{-5}$	1.87
2.0	$1.26\times {10}^{-3}$	2.08
5.0	$2.57\times {10}^{-3}$	2.16
Solute diffusion
0.0 (eq. (29))	$3.87\times {10}^{-6}$	$1.62+3.17=4.79$
0.5	$7.56\times {10}^{-7}$	4.67
1.0	$7.80\times {10}^{-7}$	4.67
2.0	$7.89\times {10}^{-7}$	4.66
5.0	$6.75\times {10}^{-7}$	4.59

Equations

$$\mathcal{H}=\sum _i\sum _{\alpha ,\beta }{n}_{\alpha \text{-}\beta }^{(i)}{ϵ}_{\alpha \text{-}\beta }^{(i)}$$

$$\mathcal{H}=\sum _{n,m}\sum _{\alpha ,\beta }{\sigma }_{\alpha }^n{\sigma }_{\beta }^m$$

$$p_{ij}=\exp \left(-\frac{\Delta {\mathcal{H}}_{ij}-N_{\mathrm{B}}\Delta \mu }{k_{B}T}\right)$$

$$p_{ij}=\exp \left(-\frac{\Delta {\mathcal{H}}_{ij}+2\Delta \mu }{k_{B}T}\right)$$

$$p_{ij}=\exp \left(-\frac{\Delta {\mathcal{H}}_{ij}}{k_{B}T}\right)$$

$$r_{ij}=\nu \exp \left(-\frac{\Delta E_{ij}}{k_{B}T}\right)$$

$$\Delta E_{ij}=\sum _p{ϵ}_{\alpha \text{-}p}^{sp}-\sum _q{ϵ}_{\alpha \text{-}q}^{(i)}-\sum _{r\ne \alpha }{ϵ}_{\text{V-}r}^{(i)}+\sum \Delta E_{ij}^{\mathrm{non}\text{-}\mathrm{broken}}$$

$$\Delta E_{ij}=\{\begin{array}{ll}{E}_m+\Delta {\mathcal{H}}_{ij},& \text{if}\Delta {\mathcal{H}}_{ij}>0\\ E_m,& \text{if}\Delta {\mathcal{H}}_{ij}<0\end{array}$$

$${ϵ}^{(2)}={ϵ}^{(1)}{\left(\frac{r_{2nn}}{r_{1nn}}\right)}^{-6}$$

$$E_{coh}^{\text{A}}=-\frac{z_1}{2}{ϵ}_{\text{A-A}}^{(1)}-\frac{z_2}{2}{ϵ}_{\text{A-A}}^{(2)}$$

$$E_{coh}^{\text{B}}=-\frac{z_1}{2}{ϵ}_{\text{B-B}}^{(1)}-\frac{z_2}{2}{ϵ}_{\text{B-B}}^{(2)}$$

$$\Delta H^{\mathrm{mix}}=$$

$$+$$

$${\Omega }_s^{(1)}={ϵ}_{\text{A-B}}^{(1)}-\frac{1}{2}\left({ϵ}_{\text{A-A}}^{(1)}+{ϵ}_{\text{B-B}}^{(1)}\right)$$

$${\Omega }_s^{(2)}={ϵ}_{\text{A-B}}^{(2)}-\frac{1}{2}\left({ϵ}_{\text{A-A}}^{(2)}+{ϵ}_{\text{B-B}}^{(2)}\right)$$

$$X(1-X){\Omega }_s^{\ast }=\Delta H^{\mathrm{mix}}+(1-X)E_{coh}^{\text{A}}+X{E}_{coh}^{\text{B}}$$

$${\Omega }_s^{\ast }=w_0+w_{1}X$$

$$E_f^{\text{V}}=E_{coh}^{\text{A}}-z_1{ϵ}_{\text{A-V}}^{(1)}-z_2{ϵ}_{\text{A-V}}^{(2)}$$

$$E_b^{{\text{B}}_m{\text{V}}_n}=E_f^{{\text{B}}_m{\text{V}}_n}-m{E}_f^{\text{B}}-n{E}_f^{\text{V}}$$

$$E_b^{(a)}={ϵ}_{\text{B-V}}^{(1)}+{ϵ}_{\text{A-A}}^{(1)}-{ϵ}_{\text{A-B}}^{(1)}-{ϵ}_{\text{A-V}}^{(1)}$$

$$E_b^{(b)}={ϵ}_{\text{B-V}}^{(2)}+{ϵ}_{\text{A-A}}^{(2)}-{ϵ}_{\text{A-B}}^{(2)}-{ϵ}_{\text{A-V}}^{(2)}$$

$$\begin{array}{rl}{E}_b^{(c)}=& 2{ϵ}_{\text{B-V}}^{(1)}+{ϵ}_{\text{B-B}}^{(2)}+2{ϵ}_{\text{A-A}}^{(1)}+{ϵ}_{\text{A-A}}^{(2)}-2{ϵ}_{\text{A-V}}^{(1)}-2{ϵ}_{\text{A-B}}^{(1)}\\ & -2{ϵ}_{\text{A-B}}^{(2)}+14\Delta {ϵ}_{\text{A-B}}^{(1)}+10\Delta {ϵ}_{\text{A-B}}^{(2)}\end{array}$$

$$E_b^{(d)}=2{ϵ}_{\text{B-V}}^{(2)}+2{ϵ}_{\text{A-A}}^{(2)}-2{ϵ}_{\text{A-B}}^{(2)}-2{ϵ}_{\text{A-V}}^{(2)}$$

$$E_{b,\text{1nn}}^{\text{V-V}}={ϵ}_{\text{A-A}}^{(1)}+{ϵ}_{\text{V-V}}^{(1)}-2{ϵ}_{\text{A-V}}^{(1)}$$

$$E_{b,\text{2nn}}^{\text{V-V}}={ϵ}_{\text{A-A}}^{(2)}+{ϵ}_{\text{V-V}}^{(2)}-2{ϵ}_{\text{A-V}}^{(2)}$$

$$E_f^{\text{AA}}=-4{ϵ}_{\text{A-A}}^{(1)}-3{ϵ}_{\text{A-A}}^{(2)}+8{ϵ}_{\text{AA-A}}^{(1)}+6{ϵ}_{\text{AA-A}}^{(2)}$$

$$E_f^{\text{AB}}=-4{ϵ}_{\text{A-A}}^{(1)}-3{ϵ}_{\text{A-A}}^{(2)}+8{ϵ}_{\text{AB-A}}^{(1)}+6{ϵ}_{\text{AB-A}}^{(2)}$$

$$E_{b,\text{1nn}}^{\text{AA-B}}$$

$$E_{b,\text{1nn}}^{\text{AB-B}}$$

$$E_{b,\text{1nn}}^{\text{AA-AA}}$$

$$E_{b,\text{1nn}}^{\text{AA-AB}}$$

$$E_{b,\text{1nn}}^{\text{AB-AB}}$$

$$z_1^{sp}{ϵ}_{\alpha \text{-A}}^{sp(1)}+z_2^{sp}{ϵ}_{\alpha \text{-A}}^{sp(2)}=E_m+\sum _{n,q}{ϵ}_{X\text{-}q}^{(n)}+\sum _{n,r\ne X}{ϵ}_{V\text{-}r}^{(n)}$$

$$\eta =N_B^{-1}\sum _i^{N_B}\left(1-\frac{x_i(\mathrm{A})}{1-X}\right)$$

$$D(T)=\nu f{\delta }^2\exp \left(-\frac{E_a}{kT}\right)$$

$$D=\frac{⟨\Delta r^2⟩}{6\Delta t}$$

Full text

W-Re alloys; Solute precipitation; Neutron irradiation; kinetic Monte Carlo

I Introduction

Tungsten is the prime candidate material in magnetic fusion energy devices due to its high strength and excellent high temperature properties Zinkle and Ghoniem (2000); Rieth et al. (2013); Fitzgerald and Nguyen-Manh (2008); Becquart and Domain (2009). Upon fast neutron irradiation in the 600-1000∘C temperature range, W transmutes into Re by the way of beta decay reactions at a rate that depends on the neutron spectrum and the position in the reactor. For the DEMO (demonstration fusion power plant) reactor concept, calculations show that the transmutation rate is 2000 and 7000 atomic parts per million (appm) per displacements per atom (dpa) in the divertor and the equatorial plane of the first wall, respectively (where damage, in each case, accumulates at rates of 3.4 and 4.4 dpa/year) Gilbert and Sublet (2011); Gilbert et al. (2012). The irradiated microstructure initially evolves by accumulating a high density of prismatic dislocation loops and vacancy clusters, approximately up to 0.15 dpa He et al. (2006); Tanno et al. (2009); Hasegawa et al. (2011); Hu et al. (2016). Subsequently, a void lattice emerges and fully develops at fluences of around 1 dpa. After a critical dose that ranges between 5 dpa for fast ($>$1 MeV) neutron irradiation Hasegawa et al. (2011) and 2.2 dpa in modified target rabbits in the HFIR Hu et al. (2015, 2016), W and W-Re alloys develop a high density of nanometric precipitates with acicular shape at Re concentrations well below the solubility limit Hasegawa et al. (2011); Hu et al. (2016). The structure of these precipitates is consistent with $\sigma$ (W7Re6) and $\chi$ (WRe3) intermetallic phases, which under equilibrium conditions only occur at temperatures and Re concentrations substantially higher than those found in neutron irradiation studies Cottrell (2004). A principal signature of the formation of these intermetallic structures in body-centered cubic (bcc) W is the sharp increase in hardness and embrittlement Tanno et al. (2009); Hasegawa et al. (2011); Hu et al. (2016). Qualitatively similar observations have been recently made in W-2Re and W-1Re-1Os alloys subjected to heavy ion irradiation Xu et al. (2015, 2017), clearly establishing a link between primary damage production and Re precipitation.

Precipitation of nonequlibrium phases in irradiated materials is commonplace. The standard theory of irradiation damage includes radiation enhanced diffusion (RED) and radiation induced precipitation (RIP) as mechanisms that can drive the system out of equilibrium due to the onset of point defect cluster fluxes towards defect sinks Was (2007); Dienes and Damask (1958); Cauvin and Martin (1979). Within this picture, Re precipitation in W or W-Re alloys under irradiation would then, in principle, be unsurprising were it not for the fact that Re clustering is seen to occur at concentrations still below the solubility limit even after RED has taken place. In spite of this, recent work using energy models based on the cluster expansion formalism for the W-Re system, and fitted to density functional theory (DFT) calculations, have revealed a direct relationship between excess vacancy concentrations and the formation of Re solute-rich clusters Wrobel et al. (2016). These calculations are substantiated by recent neutron irradiation experiments of pure W at 900*∘C up to 1.6 dpa in the HFR in Petten Klimenkov et al. (2016). Post-irradiation examination of the irradiated specimens reveals the formation of a fine distribution of voids with average 5 nm size surrounded by Re-rich clouds. However, the relative concentration of Re around the voids is still on the order of 12-18% (from a nominal overall concentration of 1.4% from transmutation), well below the precipitation limit of Re in W at 900∘C. However, in the ion beam irradiation experiments of W-2Re alloys by Xu et al. at 300 and 500∘*C, Re-rich clusters with bcc structure are seen to form with concentrations between 12 and 30% Re with no indication of vacancies forming part of the clusters Xu et al. (2015, 2017). Another piece of evidence against a strong association between vacancies and Re atoms comes from irradiation tests of W-Re alloys performed at EBR-II in the 1970s and 80s Matolich et al. (1974); Sikka and Moteff (1974); Williams et al. (1983); Herschitz and Seidman (1984, 1985). In these studies, the presence of Re was seen to suppress swelling, which would seem to suggest a decoupling between vacancy clusters and Re atoms. Clearly, equilibrium thermodynamics involving vacancies alone may not suffice to explain the precipitation tendencies in irradiated W-Re alloys.

All this is suggestive of alternative solute transport mechanisms that may be unique to W-Re systems. Indeed, several recent studies using electronic structure calculations have independently reported a peculiar association between self-interstitial atoms (SIA) and Re solutes that results in very high solute transport efficacy Suzudo et al. (2014, 2015); Gharaee et al. (2016). This mechanism consists of a series of mixed dumbbell rotations and translations such that the mixed nature of the dumbbell is preserved and solutes can be transported over long distances without the need for vacancy exchanges. Furthermore, this mechanism effectively transforms one-dimensional SIA diffusion into a 3D mixed-dumbbell transport process at activation energies considerably lower than that of vacancy diffusion. The objective of this paper is to study the kinetics of Re precipitation in irradiated W accounting for both vacancy and mixed-interstitial solute transport. To this effect, we develop a lattice kinetic Monte Carlo (kMC) model of alloy evolution parameterized solely using first principles calculations. We start in Section II by describing the essential elements of our kinetic model as well as the parameterization effort based on DFT calculations. In Section III we provide our main results, including semi-grand canonical Monte Carlo calculations of ternary W-Re-vacancy and W-Re-SIA systems, and kMC simulations and analysis of the Re-precipitate nucleation and growth. We finish with a discussion of the results and the conclusions in Section IV.

II Theory and methods

II.1 Energy model

The energy model employed throughout this work is a cluster expansion Hamiltonian based on pair interactions truncated to the $2^{\rm nd}$-nearest neighbor (2nn) shell:

[TABLE]

where $(i)$ specifies the type of nearest-neighbor interaction (first or second), $\alpha$ and $\beta$ refer to a pair of lattice sites, separated by a distance specified by the index $i$, $n_{\alpha\text{-}\beta}$ denotes the number of occurrences (bonds) of each $\alpha\text{-}\beta$ pair, and $\epsilon_{\alpha\text{-}\beta}$ is bond energy. In a previous work, we have shown how this Hamiltonian can be reduced to a generalized Ising Hamiltonian involving solvent and solute atoms (A and B), vacancies (V), and pure and mixed interstitials (AA, BB, and AB) Huang and Marian (2016). The Hamiltonian is then expressed as a sum of polynomial terms of various degrees involving spin variables $\sigma_{\alpha}$ and $\sigma_{\beta}$ in the manner of the Ising model:

[TABLE]

One of the advantages of using this notation is that the values assigned to the spin variables conserve the number of atoms $N$ of the system. We refer the reader to ref. Huang and Marian (2016) for more details about this notation. In this paper we focus on the parameterization exercise for irradiated W-Re alloys111With A: W atoms; B: Re atoms; V: vacancies, AA: W-W dumbbell (or self-interstitial atom); BB: Re-Re dumbbell; AB: mixed W-Re dumbbell..

II.2 Semi-Grand Canonical Monte Carlo for AB systems

The thermodynamic phase diagram of the W-Re system can be studied using semi-grand canonical Monte Carlo (SGMC) calculations as a function of temperature and solute concentration Binder et al. (1981); Dunweg and Landau (1993); Pareige et al. (1999); Tavazza et al. (2004); Cannavacciuolo and Landau (2005); Biborski et al. (2010). We seek to minimize the thermodynamic potential of the semi-grand canonical ensemble, characterized by a constant temperature $T$, a constant number of particles $N$, and a constant chemical potential $\mu$. In each SGMC step, an A atom is randomly flipped into a B atom and the new state is accepted with probability:

[TABLE]

where $\Delta{\cal H}_{ij}$ is the energy difference between the initial and final states, $i$ and $j$, $N_{\rm B}=NX$ is the number of solute atoms ($X$: solute concentration), $\Delta\mu$ is the change in chemical potential per atom after the transition, and $k_{B}$ is Boltzmann’s constant. In this work, each transition is defined by changing the chemical nature of one atom chosen at random. In terms of the change in spin variable (in the notation of the generalized Ising Hamiltonian, cf. eq. (2)), this always results in a change of $\delta\sigma=-2$, such that eq. (3) can be simplified to:

[TABLE]

In the calculations, the chemical potential difference $\Delta\mu$ and the temperature $T$ are input variables, while the solute composition $X$ and the equilibrium configurations are obtained when convergence is reached.

II.3 Metropolis Monte Carlo calculations of ABV system configurations

During irradiation, the introduction of large amounts of defects has the potential to impact the thermodynamics of the system. It is therefore of interest to calculate phase diagrams with fixed defect concentrations using equilibrium (Metropolis) Monte Carlo. Defect concentrations are not thermodynamically equilibrated under irradiation –the number of vacancies or interstitials is not controlled by the chemical potential–, and so the AB system must be considered in conjunction with a fixed defect concentration. Take the case of vacancies for example, to properly obtain converged nonequilibrium configurations of ABV systems, we employ a flip and swap approach: (i) initially a system consisting of A atoms and a random distribution of vacancies is considered; (ii) a lattice point is selected at random; (iii) if that lattice point corresponds to an atom, a SGMC step is carried out, resulting in a change in the relative concentrations of A and B; if it, on the contrary, corresponds to a vacant site, then a canonical Monte Carlo step is carried out, leaving $X$ unchanged, and the vacancy is swapped with a randomly-selected atom. This trial swap is then accepted according to the Boltzmann distribution:

[TABLE]

In this fashion, equilibrated AB alloys containing a fixed vacancy concentration are obtained, from which one can determine the changes relative to the thermodynamic equilibrium configurations. Although interstitials are much higher in energy than vacancies (so that only very small concentrations need be explored), the procedure for the ABI system is identical to that of the ABV system.

II.4 Kinetic Monte Carlo simulations of ABVI systems

The kinetic evolution of W-Re alloys under irradiation is studied using standard lattice kMC. The system is evolved by events involving atomic jumps and time is advanced according to the residence-time algorithm Young and Elcock (1966). Jump rates are calculated as:

[TABLE]

where $\nu$ is an attempt frequency and $\Delta E_{ij}$ is the activation energy to jump from state $i$ to state $j$.

II.4.1 Vacancy migration model

Several models have been proposed to describe the activation energy based on different interpretations of the atomic migration process (see, e.g. Soisson et al. (2010) and Huang and Marian (2016) for recent reviews). In this work, the activation energy of vacancy jump is modeled by the saddle-point energy model (or cut-bond model) Soisson et al. (1996); Soisson and Fu (2007); Martínez et al. (2012); Senninger et al. (2016), according to which $E_{ij}$ is given by the energy difference of the configuration when the jumping atom is at saddle point and the initial configuration:

[TABLE]

where $\alpha$ is the jumping atom, V is the vacancy, and $\epsilon^{sp}$ are the bond energies between the atom at the saddle point and the neighboring atoms. The first term on the r.h.s. of eq. (7) reflects the energy of the jumping atom at the saddle point. In this work, we consider interactions up to 2nn distances for this term222In the saddle-point configuration for vacancy migration, there are six 1nn bonds and six 2nn bonds, compared with eight and six for a lattice point configuration.. The second and third terms on the r.h.s. of eq. (7) are the energies of the jumping atom and the vacancy at the initial state $i$. Finally, the fourth term gives the energy difference between state $i$ and $j$ for the non-broken bonds due to local solute concentration changes. The dependence of bond coefficients on local solute concentration will be discussed in Section II.5.

II.4.2 Interstitial defect migration model

Here we consider self-interstitial atoms of the AA type, and mixed-interstitials AB. Due to their rarity, BB interstitials are omitted in our calculations. In bcc metals, AA SIAs are known to migrate athermally in one dimension along $\langle 111\rangle$ directions, with sporadic rotations to other $\langle 111\rangle$ orientations. These processes, however, are treated separately. In contrast to vacancy migration, activation energies of interstitial jumps are calculated using the direct final-initial system energy model Djurabekova et al. (2007); Vincent et al. (2008a, b); Reina et al. (2011):

[TABLE]

In addition, we include a bias due to the well-known phenomenon of correlation, by which a forward jump is slightly more likely to occur than a backward jump. This is reflected in a correlation factor $f$ computed as the ratio of forward to backward jumps Zhou et al. (2013), which in our simulations is temperature dependent. Rotations between $\langle 111\rangle$ directions are simply characterized by an activation energy equal to a rotation energy $E_{r}$.

For their part, as pointed out in Section I, recent DFT studies have revealed a new migration mechanism for mixed dumbbells in W alloys. This mechanism involves an non-dissociative sequence of rotations and translations such that the solute atom is always part of the mixed dumbbell (in contrast with the intersticialcy or ’knock-on’ mechanism commonly displayed by SIAs) Suzudo et al. (2014, 2015); Gharaee and Erhart (2015); Gharaee et al. (2016). This effectively makes AB interstitials move in three dimensions with 2nn jumps along $\langle 100\rangle$ directions. Calculations for the W-Re system have shown that the migration energy in this case is very low, on the order of one tenth of an eV. As we shall see, this has extraordinary implications on the kinetic evolution of irradiated W-Re alloys.

II.4.3 Spontaneous events: recombination and absorption

Any recombination event occurs spontaneously (no sampling involved) when the distance between an interstitial defect and a vacancy is within the $3^{\rm rd}$ nearest neighbor shell. Another reaction considered to be instantaneous is the transition of a SIA into an AB dumbbell when it encounters a solute atom: AA+B$\rightarrow$AB+A. This is because the binding energy between a SIA and a Re solute atom has been calculated to be $-0.8$ eV (negative binding energies represent attraction). The distance for this transformation is set to be equal to the 1nn separation.

Defect absorption represents another type of spontaneous event. Absorption can occur at sinks, such as a plane located in a stationary position within the simulation box Soisson (2006), or a free surface Huang and Marian (2016). Sinks can potentially act also as defect emitters, as in the case of grain boundaries, dislocations, and free surfaces in real microstructures. Details about the implementation of these processes can be found in ref. Huang and Marian (2016).

II.4.4 Frenkel pair generation

In this work, defects are generated as Frenkel pairs at a prescribed rate set by the damage rate. To insert a defect pair, two atomic sites are chosen at random, one is replaced by a vacancy and the other with an interstitial formed by an A atom and the lattice atom.

II.5 Parameters

There are five distinct atomic species used in this work: W atoms (A), Re atoms (B), vacancies (V), SIAs (AA), and mixed-interstitials (AB). As mentioned above, our energy model consists of pairwise interactions up to the 2nn shell. After discounting interstitial-vacancy bonds, this amounts to 26 different types of bonds (13 for each nearest neighbor shell), all of which must be obtained using first-principles calculations. Moreover, as discussed by Martinez et al. and Senninger et al.  Martínez et al. (2012); Senninger et al. (2016), several of these bond energies are sensitive to the local solute concentration and must be computed on the fly in each Monte Carlo step. Following Warczok et al. Warczok et al. (2012), we reduce the number of unknowns from 26 to 13 by partitioning bond energies according the following relation:

[TABLE]

which is used unless both bond energies can be explicitly calculated. For the bcc lattice, this results in $\epsilon_{\alpha\text{-}\beta}^{(1)}/\epsilon_{\alpha\text{-}\beta}^{(2)}=0.421875$ for regular bond coefficients, and $\epsilon_{\alpha\text{-}\beta}^{sp(1)}/\epsilon_{\alpha\text{-}\beta}^{sp(2)}=0.194052$ for saddle-point bond coefficients.

The local solute concentration is always computed up to the 2nn shell. Next we describe the parameterization procedure for each set of bond energies.

II.5.1 W-Re parameters

The W-Re bond coefficients are $\epsilon_{\text{A-A}}$, $\epsilon_{\text{B-B}}$, and $\epsilon_{\text{A-B}}$. They determine the thermodynamic equilibrium phase diagram of the alloy. The bond coefficients of $\epsilon_{\text{A-A}}$ and $\epsilon_{\text{B-B}}$ are obtained from the cohesive energy:

[TABLE]

where $z_{1}$ and $z_{2}$ are coordination numbers for the 1nn and 2nn shells, respectively. The cohesive energies calculated using DFT are given in Table 1.

333With xc-energy correction from Ref. Kato et al. (2011)

The coefficient for the A-B bond is obtained from the enthalpy of mixing of W-Re, $\Delta H^{\rm mix}$, which can be written within the Bragg-Williams approximation Bragg and Williams (1934, 1935); Williams (1935) as:

[TABLE]

where $X$ is the global solute concentration, and $\Omega_{s}$ is the heat of solution, defined as:

[TABLE]

Combining eqs. (10) and (11), $\Delta E^{mix}$ can be expressed as:

[TABLE]

where $\Omega_{s}^{*}=z_{1}\Omega_{s}^{(1)}+z_{2}\Omega_{s}^{(2)}$. To ascertain the dependence on the solute concentration of the heat of solution, we fit the l.h.s. of eq. (14) to the data points for the mixing enthalpies as a function of $X$ calculated in our previous work Gharaee et al. (2016). The best fit, shown in Figure 1, is achieved when $\Omega_{s}^{*}$ is expressed a linear function of the concentration:

[TABLE]

with $w_{0}=-0.1571$ and $w_{1}=-0.2311$. The negative values of $w_{0}$ and $w_{1}$ suggest a strong tendency towards ordering, which becomes larger as the solute concentration increases. Combining eqs. (9), (10), (12), and (14), one can obtain the values of $\Omega_{s}^{(1)}$, $\Omega_{s}^{(2)}$, $\epsilon_{\text{A-B}}^{(1)}$, and $\epsilon_{\text{A-B}}^{(2)}$. A non-constant $\Omega_{s}^{*}$ effectively implies that $\epsilon_{\text{A-B}}$ is also a function of the concentration. Moreover, to reflect local composition variations in the W-Re alloys, we make the assumption that the dependence of $\epsilon_{\text{A-B}}^{(1)}$ and $\epsilon_{\text{A-B}}^{(2)}$ on $X$ can be transferred to the local environment of each atom, such that both bond energy coefficients are functions of the local composition, which we term $x$, and must be computed on the fly for each solute atom in the system.

II.5.2 Vacancy parameters

The vacancy bond coefficients are $\epsilon_{\text{A-V}}$, $\epsilon_{\text{B-V}}$, and $\epsilon_{\text{V-V}}$. $\epsilon_{\text{A-V}}$ can be readily obtained from the value of the vacancy formation energy:

[TABLE]

where $E_{f}^{\text{V}}$ is the vacancy formation energy in pure W (given in Table 1). $\epsilon_{\text{B-V}}$ can be obtained from the binding energies of V-Re configurations. The binding energy of a configuration involving $m$ solute atoms and $n$ vacancies is defined as:

[TABLE]

where the $E_{f}$ are the respective formation energies of each structure. In this work, binding energies for the four vacancy-solute configurations shown in Fig. (2) have been calculated (cf. Table 1). One can now rewrite eq. (16) as a function of the B-V bond coefficients $\epsilon_{\text{B-V}}^{(1)}$ and $\epsilon_{\text{B-V}}^{(2)}$ for each one of the configurations in the figure:

[TABLE]

where $\Delta\epsilon_{\text{A-B}}^{(m)}$ is the change in $\epsilon_{\text{A-B}}^{(m)}$ due to the local solute concentration change resulting from the vacancy jump.

To define the dependence on $x$ of $\epsilon_{\text{B-V}}^{(1)}$, we must consider two factors. First, our DFT calculations show that $\epsilon_{\text{A-V}}^{(1)}>\epsilon_{\text{B-V}}^{(1)}$. Second, $x$ of $\epsilon_{\text{B-V}}^{(1)}$ is seen to increase with local concentration. Both of these conditions are satisfied by assuming a dependence such as $\epsilon_{\text{B-V}}^{(1)}(x)=\epsilon_{\text{A-V}}^{(1)}-ax^{-1}$, where $a$ is a fitting constant. As well, $\epsilon_{\text{B-V}}^{(2)}$ is seen to independently increase with concentration, such that $\epsilon_{\text{B-V}}^{(2)}(x)=bx+c$, where $b$ and $c$ are fitting parameters

$\epsilon_{\text{V-V}}$ can be readily calculated by considering the binding energy of a di-vacancy:

[TABLE]

It is interesting to note that, in accordance with several other studies Becquart and Domain (2007); Ventelon et al. (2012); Muzyk et al. (2011); Oda et al. (2014); Setyawan et al. (2016), $E_{b,\text{2nn}}^{\text{V-V}}$ takes a positive value (cf. Table 1), indicating repulsion between vacancies that are at 2nn distances of each other.

II.5.3 Interstitial defect parameters

The interstitial bond coefficients include $\epsilon_{\text{AA-A}}$, $\epsilon_{\text{AB-A}}$, $\epsilon_{\text{AA-B}}$, $\epsilon_{\text{AB-B}}$, $\epsilon_{\text{AA-AA}}$, $\epsilon_{\text{AA-AB}}$, and $\epsilon_{\text{AB-AB}}$. $\epsilon_{\text{AA-A}}$ and $\epsilon_{\text{AB-A}}$ are calculated directly from the formation energies of SIAs and mixed dumbbells:

[TABLE]

The other bond coefficients are obtained from various binding energies:

[TABLE]

These formation and binding energies are all taken from the literature444The only exception being the binding energy between an AA and an AB interstitial, which is assumed to be equal to the binding energy between two AA..

All the bond energy coefficients, the equation used for their calculation, and the source of the numbers are compiled in Table 2.

II.5.4 Migration parameters

The attempt frequency ($\nu$ in eq. (6)) used for vacancy jumps in this work is set to be equal to Debye frequency of W, or $6.5\times 10^{12}$ Hz Gilbert et al. (2008), while for interstitials we use a value of $1.5\times 10^{12}$ Hz Zhou et al. (2013).

From eq. (7), there are six different saddle-point bond coefficients: $\epsilon_{\text{A-A}}^{sp(m)}$, $\epsilon_{\text{A-B}}^{sp(m)}$, $\epsilon_{\text{A-V}}^{sp(m)}$, $\epsilon_{\text{B-A}}^{sp(m)}$, $\epsilon_{\text{B-B}}^{sp(m)}$, and $\epsilon_{\text{B-V}}^{sp(m)}$, where $m$= 1nn, 2nn. In this notation, $\epsilon_{\alpha\text{-}\beta}^{sp(m)}$ represents the energy of the bond between the atom at the saddle point $\alpha$ and its closest lattice neighbor $\beta$. This means $\epsilon_{\alpha\text{-}\beta}^{sp(m)}\neq\epsilon_{\beta\text{-}\alpha}^{sp(m)}$.

The saddle-point bond coefficients connected to a lattice atom A (W atom), $\epsilon_{\alpha\text{-A}}^{sp(m)}$, can be calculated as:

[TABLE]

where $z_{1}^{sp}$ and $z_{2}^{sp}$ are the numbers of $1^{\rm st}$- and $2^{\rm nd}$ nearest neighbor of an atom at the saddle point, which are both equal to 6 for the bcc lattice, and $E_{m}$ is the migration energy. The term $\Delta E_{ij}^{\rm non\text{-}broken}$ in eq. (7) is zero here since no solute concentration change is involved in an A-atom jump. $\epsilon_{\alpha\text{-A}}^{sp(2)}$ is obtained from $\epsilon_{\alpha\text{-A}}^{sp(1)}$ using eq. (9). Vacancy bonds are calculated in a similar manner.

To calculate the saddle-point bond coefficients pertaining to B (Re) atoms, $\epsilon_{\alpha\text{-B}}^{sp(m)}$, one must consider local solute concentration changes. To this end, we resort to the configurations shown in Fig. 2(c). The A-B saddle-point coefficients $\epsilon_{\text{A-B}}^{sp(m)}$ are obtained from A-atom jumps, labeled ‘1’ and ‘2’ in Fig. 2(c), into the vacant site. The B-B saddle-point coefficient $\epsilon_{\text{B-B}}^{sp(1)}$ is computed assuming a B-atom (labeled ‘3’ in the figure) jump into the vacancy. Equation (9) is then used to obtain the 2nn coefficients. All the necessary DFT calculations to calculate the saddle-point bond coefficients were performed as part of the present work, and are given in Table (3).

The migration energies of SIA and mixed-interstitials, the activation energy for SIA rotation, as well as the correlation factors at different temperatures are taken from the literature, and listed in Table (4).

II.5.5 DFT calculations

Density functional theory calculations were carried out using the projector augmented wave (PAW) method Blöchl (1994); Kresse and Joubert (1999) as implemented in the Vienna ab-initio simulation package Kresse and Hafner (1993, 1994); Kresse and Furthmüller (1996a, b). Since interstitial configurations involve short interatomic distances “hard” PAW setups that include semi-core electron states were employed with a plane wave energy cutoff of 300 eV.

Exchange and correlation effects were described using the generalized gradient approximation Perdew et al. (1996) while the occupation of electronic states was performed using the first order Methfessel-Paxton scheme with a smearing width of 0.2 eV. The Brillouin zone was sampled using $5\times 5\times 5$ $\vec{k}$-point grids. (A detailed discussion of the effect of different computational parameters on the results can be found in Ref. Gharaee and Erhart (2015)). All structures were optimized allowing full relaxation of both ionic positions and cell shape with forces converged to below 10 meV/Å. Migration barriers were computed using $4\times 4\times 4$ supercells and the climbing image-nudged elastic band method with three images Henkelman et al. (2000).

III Results

III.1 Structural phase diagrams

Although our energy model includes thermodynamic information reflective of the phase stability of W-Re alloys, the model consists of a rigid lattice with bcc structure and is thus suitable only for a given, well-defined, concentration range. Our DFT calculations yield bond energies that are consistent with stable solid solutions from zero to approximately 40% at. Re Gharaee et al. (2016). This is confirmed by way of sgMC simulations performed as a function of composition and temperature in $64\times 64\times 64$ computational cells. Figure 3 shows the set of stable compositions obtained as a function of the chemical potential for several temperatures.

The figure shows a clear jump in the Re concentration at a temperature of approximately 100 K. This is typically indicative of a phase boundary, as two distinct phases characterized by widely different solute concentrations seem to coexist at the same temperature and chemical potential. This may be indicative of a miscibility gap between Re concentrations of a few percent and approximately 50 at.%, i.e. beyond the thermodynamic validity of our rigid lattice model. To further characterize the configurations obtained, we calculate the short-range order (SRO) of the configurations obtained according to the Warren-Cowley parameter Cowley (1950):

[TABLE]

which gives the SRO parameter $\eta$ of Re atoms w.r.t. matrix W atoms, with $x_{i}({\rm A})$ being the fraction of A atoms surrounding each solute atom $i$. The sum extends to all B atoms in the system.

According to this definition, $\eta>0$ implies phase separation, $\eta=0$ represents an ideal solid solution, and $\eta<0$ indicates ordering. However, the SRO parameter of a random solution has a range of $-0.003$ to 0.003 regardless of solute composition due to the random occurrence of dimers, trimers, etc. Figure (4) shows the equilibrium SRO as a function of $X$ for several temperatures. As the figure shows, the SRO parameter is near zero for dilute systems, and gradually becomes negative as the concentration increases. Based on the figure we conclude that equilibrium W-Re systems with up to $\approx 40$ at.% solute content are consistent with random solid solutions with a weak tendency to ordering at higher concentrations. The corresponding $T$-$X$ phase diagram is provided in Figure 5.

III.1.1 Effect of vacancies on phase diagram

It is well known that non-equilibrium concentrations of defects can alter the thermodynamic behavior of an alloy. For the W-Re system, Wrobel et al. have studied the ternary W-Re-vacancy system and found that Re clustering occurs in the presence of non-thermodynamic vacancy concentrations Wrobel et al. (2016). Clusters appear as semi-ordered structures of alternating solute and vacancy planes –a necessity given the short-range repulsion between Re atoms on the one hand, and vacancies on the other (cf. Table 1). Next, we carry out a similar study involving various vacancy concentrations, temperatures, and solute concentrations to obtain structural phase diagrams such as that shown in Fig. 5. Each configuration is optimized by combining sgMC steps with energy minimization steps following the process described in Sec. II.3. Figure 6 shows the diagrams for vacancy concentrations of $C_{v}=0.01,~{}0.1,~{}0.2,~{}0.5$ at.% using $64\times 64\times 64$ primitive cells.

As a representative example, Figure 7a shows the equilibrated configuration at 600 K, 1.8 at.% Re (which occurs for $\Delta\mu=0.26$), and $C_{v}=0.5$ at%. The figure shows several Re-vacancy clusters with an ordered structure, consistent with the study by Wrobel et al. Wrobel et al. (2016). Due to their ordered structure, these solute-vacancy clusters form only at Re concentrations that are commensurate with the vacancy concentration in the system, i.e. at values of $X\lesssim 0.04$ in most cases.

III.1.2 Effect of interstitial defects on the phase diagram

Although vacancy concentrations such as those considered in this section are several orders of magnitude larger than the vacancy concentration in thermal equilibrium, one can expect such numbers under far-from-equilibrium conditions such as under high-dose or high-dose rate irradiation. The case is much more difficult to make for SIAs due to their much higher formation energy (3.2 vs. 10.2 eV, to take two representative numbers Gharaee and Erhart (2015)). However, given the inclination of single interstitials to convert into mixed dumbbells in the presence of solute, it is of interest to repeat the same exercise of looking at the clustering propensity of Re in such cases. The results are shown in Figure 8 for a defect concentration of 0.1 at.%. The diagram reveals a stronger clustering tendency when interstitials are present compared to vacancies. Such an effect originates from both more attractive binding energies between mixed-interstitials and solute atoms, and between mixed-interstitials with themselves. A snapshot of the equilibrated atomistic configuration is shown in Fig. (7b), where the precipitates are seen to form platelet-like structures with a mixed dumbbell core surrounded by substitutional solute atoms.

III.2 Kinetic evolution of irradiated W-Re alloys

There are a number of factors that call for performing kMC simulations in W-Re systems.

First, equilibrium Monte Carlo calculations such as those performed in Section III.1 do not provide information about the precipitate nucleation and growth mechanisms, as well as the timescales involved. 2. 2.

Second, there is clear experimental evidence of Re-cluster formation in the absence of vacancies. Hasegawa et al. Tanno et al. (2007); Hasegawa et al. (2011) and Hu et al. Hu et al. (2016) have both reported the formation of W-Re intermetallic precipitates after high-dose, fast neutron irradiation. Moreover, recent irradiation experiments have revealed the formation of Re-rich clusters with bcc structure, i.e. prior to their conversion into $\sigma$ and/or $\chi$ precipitates. For example, Klimenkov et al. note that Re-rich particles not associated with cavities formed in neutron-irradiated single crystal W Klimenkov et al. (2016). As well, using atom-probe tomography Xu et al. have performed detailed analyses of Re-rich atmospheres in bcc W without detecting significant numbers of vacancies Xu et al. (2017). 3. 3.

New understanding regarding interstitial-mediated solute transport in W-Re alloys Suzudo et al. (2015); Gharaee et al. (2016), together with the results in Section III.1.2, call for renewed simulation efforts incorporating these new mechanisms –in particular, the three-dimensional and associative nature of Re transport via mixed-dumbbell diffusion.

These considerations motivate the following detailed study of the Re precipitation kinetics under irradiation conditions. First, however, we proceed to calculate diffusion coefficients and transport coefficients for defect species and solute atoms.

III.2.1 Calculation of diffusion coefficients

Tracer diffusion coefficients (i.e., in the absence of a concentration gradient) for vacancies, interstitials, and solute species in three dimensions are assumed to follow an Arrhenius temperature dependence:

[TABLE]

where $\nu$ is the so-called attempt frequency, $f$ is the correlation factor, $\delta$ is the jump distance, $E_{a}$ is the activation energy, and $D_{0}=\nu f\delta^{2}$ is the so-called diffusion pre-factor. Defect diffusivities can be obtained directly from this equation, with $E_{a}\equiv E_{m}$. For solute diffusion the above expression must be multiplied times the probability of finding a vacancy in one of the 1nn positions, such that $D_{0}=z_{1}\nu f\delta^{2}$ and $E_{a}=E_{m}+E_{f}^{\rm V}$. However, to calculate the diffusivities of solutes and vacancies as a function of the global solute concentration, fluctuations in local chemistry prevent us from using equations for homogeneous systems such as eq. (29), and diffusivities must be obtained by recourse to Einstein’s equation:

[TABLE]

where $\left<\Delta r^{2}\right>$ is the mean squared displacement (msd) and $\Delta t$ is the time interval. This formula assumes equilibrium defect concentrations, which are generally several orders of magnitude smaller than what a typical simulation cell can afford. For this reason, the time in eq. (30) is not directly the time clocked in the kMC simulations, $\Delta t_{\rm kMC}$. Rather, it must be rescaled by a coefficient that accounts for the difference in defect concentration Bouar and Soisson (2002); Nastar and Soisson (2012):

[TABLE]

where $C^{\rm kMC}$ and $C^{\rm eq}$ are the defect concentrations in the kMC simulations and in equilibrium, respectively. For simulations involving only one defect, $C^{\rm kMC}$ is simply equal to the inverse of the number of atoms in the computational cell, $C^{\rm kMC}=N^{-1}$, while $C^{\rm eq}=\exp(-E_{f}/k_{B}T)$, where $E_{f}$ is the instantaneous defect formation energy, i.e. calculated accounting for the local chemical environment. This is the approach used for vacancy mediated diffusion, with $E^{\rm V}_{f}=\sum_{i}\epsilon_{\text{V-}\alpha_{i}}$, where $\alpha_{i}$ symbolizes the neighboring atoms forming a bond with the vacancy. During simulations of solute and vacancy diffusion, $E^{\rm V}_{f}$ is updated in every Monte Carlo time step and time rescaling is performed on the fly. The starting configuration for all calculations involving solute atoms is the equilibrated alloy as obtained in Section III.1 using sgMC simulations. The results for the vacancy and solute diffusivities, $D_{\rm v}$ and $D_{\rm s}$, can be seen in Figure 9, while the parameters resulting from fitting the data points in the above figures to eq. (29) are collected in Table 5. While $D_{\rm v}$ displays a moderate dependence with the solute concentration, $D_{\rm v}$ is quite insensitive to it.

As discussed in Sec. II.5.4, self-interstitial migration occurs by way of fast sequences of $\langle 111\rangle$ transitions punctuated by sporadic rotations, whereas mixed dumbbell diffusion occurs via random $\langle 100\rangle$ hops in three dimensions. Interstitial diffusivities of both types can be calculated straightforwardly by using eq. (29) parameterized with the data in Table 4.

III.2.2 Calculation of transport coefficients

Within linear response theory, mass transport can be related to chemical potential gradients via Onsager’s phenomenological coefficients. The value and sign of these transport coefficients can provide important physical information about the nature of solute and defect fluxes. On a discrete lattice, the transport coefficients $L_{ij}$ coupling two diffusing species can be calculated as Senninger et al. (2016); Allnatt and Lidiard (2003):

[TABLE]

where $V$ is the total volume of the system; $\Delta r_{i}$ is the total displacement of species $i$, and $\Delta t$ is the rescaled time. Here we focus on the relationship between solutes and vacancy and solute atoms, $L_{\rm B\text{-}B}$, and $L_{\rm B\text{-}V}$, as a function of temperature and solute content. Due to the associative transport mechanism of AB interstitials, the corresponding transport coefficient relating interstitials with solute atoms is always positive and we obviate its calculation. Figure 10a shows the results for $L_{\rm B\text{-}B}$, which displays an Arrhenius temperature dependence and is always positive. The dependence with solute concentration is not significant up to 5%, with an average activation energy of 4.7 eV –very similar to the solute diffusion activation energy– and a prefactor of approximately $3.9\times 10^{20}$ m*-1*$\cdot$s*-1*. $L_{\rm B\text{-}B}$ is by definition related to the solute diffusion coefficient presented above.

In Figure 10b we plot the ratio $L_{\rm B\text{-}V}/L_{\rm B\text{-}B}$. Two observations stand out directly from the figure. First, the value of $L_{\rm B\text{-}V}$ is always negative (the exception being at 450 K, when is almost zero). This indicates a reverse coupling between solutes and vacancies, i.e. vacancy fluxes would oppose solute fluxes. The implications of this calculation will become clearer when we study solute precipitation in the next section. Second, $L_{\rm B\text{-}V}$ is on average about an order of magnitude larger than $L_{\rm B\text{-}B}$,which is to be expected for substitutional solutes moving by a vacancy mechanism.

III.2.3 Kinetic Monte Carlo simulations

To narrow down the large parametric space associated with irradiation of W-Re alloys (Re concentration, temperature, dose, dose rate, etc.), we resort to the study performed in Secs. III.1 and III.1.1. It was seen there that segregation occurs most favorably at low solute compositions. For this reason, and to enable comparison with the work by Xu et al. Xu et al. (2015, 2017), we choose a W-2.0% at. Re alloy for our study. By way of reference, this would correspond to the transmutation fraction attained after a dose of 12 dpa or 4 full-power years in DEMO’s first wall according to Gilbert and Sublet Gilbert and Sublet (2011). When relatively high concentrations of defects are present –as one might expect during irradiation– precipitation is also favored at high temperatures, so here we carry out our simulations between 1700 and 2000 K. We use box sizes of $64^{3}$ and $80^{3}$ with a damage insertion rate of $10^{-3}$ dpa per second. As shown in Appendix A, the equivalence relation that exists between both box sizes enables us to compare them directly. The three different defect-sinks discussed in Sec. II.4 are all considered here.

We first investigate the kinetic evolution of a system with no sinks. Eight independent simulations were conducted. It is seen that after an average waiting time of $\approx 21$ seconds (or $\approx 0.02$ dpa) one precipitate starts to grow in all cases. This time can be regarded as the average incubation (nucleation) time for the conditions considered in the study. Figure 11 shows the mean size from all eight cases as a function of growth time, i.e. initializing the clock after the cluster nuclei are formed regardless of the observed incubation time. The dashed line in the figure is the associated spherical growth trend, which the precipitates are seen to follow for approximately 20 s. Subsequently, growth stops at a saturation radius of 4 nm, which is seen to be the stable precipitate size.

A surface reconstruction rendition of one of the precipitates is also provided in the figure as a function of time. This depiction as a compact convex shape is not intended to represent the true diffuse nature of the cluster, and is only shown as an indication of the cluster average size and shape.

The next question we address is the solute concentration inside the precipitate. Xu et al. Xu et al. (2015, 2017) have performed detailed atom probe analyses of radial concentration profiles at 573 and 773 K and find that the precipitates that form might be better characerized as ’solute clouds’, reaching concentrations of around 30% in the center gradually declining as the radius increases. Our analysis is shown in Figure 12, with results averaged over the 8 cases tried here. The figure shows that the concentration at the precipitate core (within the inner 1.5 nanometers) surpasses 50% –the thermodynamic limit for the formation of intermetallic phases–, which could provide the driving force for such a transformation. Because our energy model is not valid above the solid solution regime, we limit the interpretation of such phenomenon however. What is clear is that the precipitates are not fully-dense, even near their center. In fact, the relative solute concentration appears to diminish near the precipitate core once the saturation point has been reached.

Finally, we address the issue of whether it is vacancy or interstitial mediated transport that is primarily responsible for solute agglomeration and the formation of Re-rich clusters. To this end, we track the evolution with time of the incremental SRO change brought about by any given kMC event during the formation of one the precipitates discussed above. The results are given in Figure 13, where contributions from SIA and mixed interstitial jumps, vacancy jumps, and Frenkel pair insertion are plotted. These results conclusively demonstrate that mixed-interstitial transport is dominant among all other events to bring solute together. Vacancies, on the other hand, serve a dual purpose. They first act as a ‘hinge’ between solute atoms that would otherwise repel, much in the manner shown in Fig. 7a. This results in an initial positive contribution to the SRO, as shown in the inset to Fig. 13, by forming dimers, trimers, or other small solute clusters. However, once a critical nucleus forms and starts to grow, vacancies reverse this behavior and act to dissolve the precipitate (differential SRO turns negative in Fig. 13), mostly by making the precipitate/matrix interface more diffuse. As expected, Frenkel pair insertion has practically no effect on the overall precipitate evolution.

The precipitate grows by a sustained capture of mixed interstitials and subsequent attraction of vacancies. This gives rise to localized recombination at the precipitate, which makes the precipitates incorporate solute atoms over time. Figure 14 shows the spatial location of the recombination events during a period of 2.0 s before, during, and after precipitate growth. The figure clearly shows that, once formed, the precipitate becomes a preferential site for recombinations, which results in further growth and eventually in saturation. Because the primary source of solute is via interstitial transport, which also brings W atoms, the precipitates are never fully compact ($x\sim 1$). Instead, maximum concentrations of around 50% are seen near the center when the precipitates reach their saturation size of 4-nm radius. As we will discuss in the next section, this is consistent with experimental measurements and observations of both coherent bcc clusters and incoherent $\sigma$ and $\chi$ phases.

Simulations performed with defect sinks in the same temperature range simply result in solute segregation in the manner described in our previous work Huang and Marian (2016). Radiation induced precipitation results from the onset of defect fluxes to the sinks, providing sufficient competition to delay the formation of bulk precipitates beyond the time scales coverable in our kMC simulations. More information is provided in the Supplementary Information.

IV Discussion and conclusions

IV.1 Mechanism of nucleation and growth

On the basis of our results, the sequence of events that leads to the formation of Re-rich precipitates in irradiated W-2Re (at%) alloys is as follows:

First, a Frenkel pair is inserted in the computational box following the procedure specified in Sec. II.4. As interstitials enter the system, they perform a one-dimensional migration until they encounter a solute atom, after which they become mixed AB dumbbells capable of transporting solute in 3D. As these mixed dumbbells diffuse throughout the lattice, they encounter substitutional solute atoms and become trapped forming a B-AB complex with a binding energy of 0.15 eV (cf. Table 1). 2. 2.

The vacancy in the Frenkel pair migrates throughout the lattice contributing to the formation of small Re complexes (dimers, trimers, tetramers, etc.). Vacancy motion does not necessary imply solute drag, as indicated by the negative value of transport coefficients in Sec. III.2.2. However, as the evidence from the Metropolis Monte Carlo simulations in Fig. 6 suggests, they can form small complexes of stable Re-V clusters. 3. 3.

The vacancy can become locally trapped in the small Re-V complexes mentioned above. However, at the high temperatures considered here, it is likely to de-trap and continue migrating until it finds the immobilized interstitial from (1), as this provides the largest thermodynamic driving force to reduce the energy of the crystal. When the vacancy and the interstitial meet, another small V-Re cluster is formed. Throughout this process, both mixed interstitial and vacancy hops are characterized by an increasing differential SRO parameter (cf. Fig. 13). 4. 4.

Eventually, one of these Re clusters grows larger than the rest due to natural fluctuations. When that happens the likelihood that the V-AB recombination will take place at that larger cluster grows. This signals the onset of the growth process, fueled by continued attraction of AB mixed dumbbells and the subsequent associated recombination. At this stage, vacancies reverse their role as solute-atom ‘hinges’ and begin to contribute to cluster dissolution (negative differential SRO parameter in Fig. 13). This results in the development of a more or less diffuse interface as the precipitate grows, which delays the next recombination event and slows down growth. 5. 5.

Although the precipitate continues to be the main pole of attraction for vacancy-interstitial recombinations (cf. Fig. 14), the system reaches a point where most of the solute is consumed into a diffuse precipitate that halts further growth. Vacancies then have more time to interact with the interface atoms before the next recombination event, which results in a smearing of the precipitate interface. In the absence of sinks, or other precipitates, the existing cluster is the sole focus of solute agglomeration, which allows it grow to its maximum size for the current alloy content of 2% Re. It is to be expected that with competing solute sinks, the precipitates might either be slightly smaller in size or less solute-dense internally.

This qualitative explanation is built on direct evidence and interpretation from our results, described in detail in Sec. III. However, to support some of the above points more explicitly, we provide additional details as Supplementary Information.

Interestingly, the essential features of our mechanism were originally proposed by Herschitz and Seidman Herschitz and Seidman (1984, 1985) on the basis of atom probe observations of neutron-irradiated W-25Re alloys. Remarkably, these authors had the intuition to propose the basic elements needed to lead to Re precipitate formation identified in our work with the more limited understanding available at the time.

IV.2 Brief discussion on the validity of our results

Our simulations are based on a highly-optimized implementation of the standard kMC algorithm. With the computational resources available to us, we can reasonably simulate systems with less than 500,000 atoms into timescales of tens of seconds. This has proven sufficient to study Re clustering at high temperatures, where vacancy mobility is high and comparable to mixed-interstitial mobility. Recall from the previous section that the formation of clusters is predicated on the concerted action of both defect species, with mixed interstitials becoming trapped at small Re clusters followed by a recombination with a vacancy that makes the cluster grow over time. Clustering and precipitation of Re in irradiated W has been seen at temperatures sensibly lower than those explored here, such as 573 and 773 K for ion-irradiated W-Re Xu et al. (2015); Edmondson et al. (2015); Xu et al. (2017), 773 and 1073 K for neutron irradiated W in HFIR Hu et al. (2016), 1173 K in neutron irradiated W in the HFR reactor Klimenkov et al. (2016), and Williams et al. at 973$\sim$1173 K in EBR-II Williams et al. (1983). The work by Hasegawa et al. in JOYO Tanno et al. (2008); Hasegawa et al. (2011) does cover –by contrast– a similar temperature range as ours. The principle is that the mechanism proposed here can be conceivably extended to lower temperatures without changes with just a timescale adjustment due to the significantly slower mobility of vacancies at those temperatures. This would require simulated times that are far too long to cover with kMC.

An intrinsic limitation of our model is that it is based on a rigid bcc lattice and cannot thus capture the transition of precipitates to the intermetallic phase. As such, our model does not necessarily reflect the true microstructural state when the local concentration surpasses 40$\sim$50%, which is when phase coexistence is expected to occur according to the phase diagram Ekman et al. (2000). However, our simulations are useful to determine the kinetic pathway towards the accumulation of Re concentrations in the vicinity of that amount. Neutron irradiation experiments such as those performed at JOYO and HFIR reveal the formation of acicular $\sigma$ and $\chi$ precipitates Williams et al. (1983); Tanno et al. (2008); Hasegawa et al. (2011); Hu et al. (2016), which presumably indicates reaching local values of Re concentration of or higher than 40$\sim$50% at the site of precipitate formation. However, in controlled ion irradiation experiments Xu et al. (2015); Edmondson et al. (2015); Xu et al. (2017) there is clear evidence that the precursor to the formation of these intermetallic precipitates are noncompact Re-rich clusters with bcc structure. We cannot but speculate how the transition from these solute-rich clusters to well-defined line compounds $\sigma$ and $\chi$ takes place (perhaps via a martensitic transformation, as in Fe-Cu systems Erhart et al. (2013)), but it is clear that it is preceded by the nucleation and growth of coherent Re clusters. In our simulations, we find that the clusters have a maximum concentration of $\approx$50% in the center, in contrast with Xu et al., who observe concentrations no larger than 30%. This disparity may simply be a consequence of the different temperatures considered relative to our simulations (773 vs 1800 K), as it is expected that the accumulation of solute by the mechanism proposed here will be accelerated by temperature.

As well, our Re clustering mechanism is predicated on the insertion of Frenkel pairs, when it is well known that fast neutron and heavy-ion irradiation generally result in the formation of clusters of vacancies and interstitials directly in displacement cascades. However, even here tungsten is somewhat of a special case. Recent work Sand et al. (2013); Yi et al. (2015); Setyawan et al. (2015) suggests that most of the defects in high-energy ($>$150 keV) cascades in W appear in the form of isolated vacancies and interstitials. This, together with the fact that most displacement cascades for non-fusion neutrons and heavy ions have energies well below the 150-keV baseline, gives us confidence that our mechanism would be operative even in such scenarios.

IV.3 Implications of our study

Beyond the obvious interest behind understanding the kinetics of Re-cluster formation in irradiated W-Re alloys, our model is useful to interpret other physical phenomena. For example, it is well known that swelling is suppressed in irradiated W-Re alloys compared to pure W Matolich et al. (1974). By providing enhanced avenues for interstitial-vacancy recombination, small Re clusters capture mixed interstitials, allowing sufficient time for vacancies to subsequently find them and suppressing the onset of swelling. Intrinsic 3D mobility of mixed dumbbells is likely to favor recombination as well. Note that this explanation for swelling suppression is different to the one proposed for Fe-Cr alloys, where 1D migration of SIAs is restrained by Cr atoms Terentyev et al. (2005).

Finally, the mechanisms proposed here refer to homogeneous nucleation, i.e. Re clustering occurs without any assistance from RED or RIP, and hence without the need for defect sinks. This is again a remarkable feature of these alloys, confirmed in several studies Herschitz and Seidman (1985); Xu et al. (2017); Klimenkov et al. (2016). As noted by Herschitz and Seidman, “The coherent precipitates were not associated with either linear or planar defects or with any impurity atoms; i.e. a true homogeneous radiation-induced precipitation occurs in this alloy”, or by Klimenkov et al., who point out that “The formation of Re-rich particles with a round shape was detected in the single crystal material. These particles were formed independently of cavities”. We leave out heterogeneous precipitation at voids, as the evidence in the literature is conflicting at this stage: discounted in some works Herschitz and Seidman (1985); Xu et al. (2017) and observed in others Klimenkov et al. (2016).

Acknowledgements

C. H., Y. Z., and J. M.’s work has been supported by the US Department of Energy’s Office of Fusion Energy Sciences, grant DE-SC0012774:0001. Computer time allocations at UCLA’s IDRE Hoffman2 supercomputer are also acknowledged. L. G. and P. E. acknowledge support from the Swedish Research Council in the form of a Young Researcher grant, and the European Research Council via a Marie Curie Career Integration Grant. Computer time allocations by the Swedish National Infrastructure for Computing at NSC (Linköping) and C3SE (Gothenburg) are gratefully acknowledged.

Appendix A Size dependence of physical time in kMC simulations

As explained in Section IV, the mechanism of formation of Re clusters requires the concerted action of both interstitials and vacancies. In order to be able to capture their formation during reasonable computational times, the temperature regime considered must be one where the mobility of both species is comparable ($1700\sim 2000$ K in our case). Then, the rate of arrival of solute atoms to a previously-nucleated Re cluster can be approximated by:

[TABLE]

where $t_{\rm FP}$ and $t_{diff}$ are the average time in between successive Frenkel-pair insertions and a characteristic diffusion time required by a vacancy and an interstitial to recombine with one another. $r_{s}$ is measured in units of atoms per unit time. At the temperatures and dose rates considered here, $t_{\rm FP}\gg t_{diff}$, such that $r_{s}\approx{t_{\rm FP}}^{-1}$. Assuming then that for each Frenkel pair inserted a minimum of one solute atom is transported:

[TABLE]

where $N_{B}$ is the total number of solute atoms in the precipitate. $r_{\rm dpa}$ in the above equation is the damage rate, expressed in units of [dpa$\cdot$s*-1*]. The precipitate volume growth rate is directly equal to the atomic volume times $r_{s}$:

[TABLE]

Assuming that the precipitate is close to spherical:

[TABLE]

And, operating, we arrive at the equation for the evolution of the precipitate radius with time:

[TABLE]

which is the equation used for fitting in Fig. 11.

Then, from eq. (34), for a given constant dpa rate, it is clear that the ratio $r_{s}(V_{1})N_{1}^{-1}=r_{s}(V_{2})N_{2}^{-1}=~{}{\rm constant}$, where $V_{1}$ and $V_{2}$ are two different box sizes. For as long as the approximation in eq. (33) is valid, then:

[TABLE]

which allows us to compare simulations done on box sizes of $64^{3}$ and $80^{3}$ directly. We emphasize that at lower temperatures, and/or high dose rate, where $t_{\rm FP}\approx t_{diff}$, this comparison is no longer valid.