<figure><img src="image_1.png" /></figure>

HAL
open science

# Development of a solute and defect concentration dependant Ising model for the study of transmutation induced segregation in neutron irradiated W-(Re,Os) systems

Matthew J. Lloyd, Enrique Martinez, Luca Messina, Duc Nguyen-Manh

► To cite this version:

Matthew J. Lloyd, Enrique Martinez, Luca Messina, Duc Nguyen-Manh. Development of a solute and defect concentration dependant Ising model for the study of transmutation induced segregation in neutron irradiated W-(Re,Os) systems. Journal of Physics: Condensed Matter, 2021, 33 (47), pp.475902. 10.1088/1361-648X/ac1ec4. cea-03563248

HAL Id: cea-03563248

https://cea.hal.science/cea-03563248

Submitted on 9 Feb 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Development of a Solute and Defect Concentration Dependant Ising Model for the Study of Transmutation Induced Segregation in Neutron Irradiated W-(Re,Os) Systems

M. J. Lloyd¹,²‡, E. Martinez³,⁴, L. Messina⁵ and D. Nguyen-Manh¹,²

¹ Department of Materials, University of Oxford, Parks Road, Oxford, OX13PH, UK

² United Kingdom Atomic Energy Authority, Culham Science Centre, Abingdon, Oxfordshire, OX143DB, UK

³ Department of Mechanical Engineering, Clemson University, Clemson, SC 29634, USA

⁴ Department of Materials Science and Engineering, Clemson University, Clemson, SC 29634, USA

⁵ CEA, DEs, IRESNE, DEC-Service d'Études et de Simulation du Comportement des Combustibles, Cadarache F-13108 Saint-Paul-Lez-Durance, France

May 2021

Keywords: Ising Model; Tungsten Rhenium Osmium Alloys; Density Functional Theory; Atomistic Monte Carlo; Transmutation Induced Segregation

**Abstract.** In this study, radiation-induced precipitation of transmutation products is addressed via the development of a new solute and vacancy concentration dependant Ising model for the W-Re-Os system. This new model includes interactions between both Os and Re atoms with vacancies, thus facilitating more representative simulations of transmutation in fusion reactor components. Local solute concentration dependencies are introduced for the W-Re, W-Os and Re-Os pair interactions. The model correctly accounts for the repulsion between small clusters of vacancies and the attraction between larger clusters/voids, via the introduction of local vacancy concentration dependant interaction coefficients between pairs of atoms and vacancies. To parameterise the pair interactions between atoms and/or vacancies, the enthalpy of mixing, $\Delta H_{mix}$, for various configurations and solute/defect concentrations, was calculated using density functional theory, within 6 binary systems: W-Re, W-Os, Re-Os, W-Vacancy, Re-Vacancy and Os-Vacancy. The new energy model was implemented into the SPPARKS Monte Carlo code, and successfully used to predict the formation of voids decorated with Re and Os solute atoms. Analysis suggests that there is a strong thermodynamic tendency for Os to bind to these voids with a comparatively weaker binding from Re atoms. The binding energies of various solute/vacancy clusters were calculated and showed that Re and Os solute atoms tend to stabilise small clusters of vacancies, increasing the attractive binding energy between the constituents.
A Solute and Defect Concentration Dependant Ising Model for W-Re-Os

2

# 1. Introduction

In the divertor and plasma facing components of a DEMOnstration fusion power station, high heat fluxes and intense neutron loading necessitates the development of novel material solutions. Restrictions based on acceptable activation levels following component removal from a reactor severely limits materials selection options for these components. Tungsten is now considered to be the leading candidate material for these applications after other options such as carbon and beryllium were found to be unsuitable. Tungsten has a high melting temperature (3965 K [1, 2]), favourable activation properties [3, 4], a good thermal conductivity (165 Wm⁻¹K⁻¹ [5]) and a high sputtering threshold under D,T and He plasma exposure (Eth=100 eV [6]). Transmutation effects under intense neutron irradiation leads to the production of Re, Os and Ta. Many studies have therefore focused on the W-Re binary system, as Re is the primary transmutation product produced under neutron irradiation under a DEMO neutron spectrum [4, 7–9], however, recent results indicate that the addition of small quantities of Os (the secondary transmutation product) can lead to major changes in the observed microstructure [10, 11]. In unirradiated W, alloying with small quantities of Re (0-10 at.%) has been shown to introduce a small solution softening effect, but increasing the Re content further results in a large increase in hardness associated with the formation of brittle intermetallic σ phase [12]. Under irradiation, precipitation of Re has been linked to a major degradation of thermal and mechanical properties [10, 13–16].

Transmutation in W is a result of a sequence of neutron absorbtion (n,γ) and neutron loss reactions (n,2n), followed by β-decay of the resultant isotopes [4]. The cross sections for these reaction depends strongly on the isotopes in the material and on the incident neutron energy spectrum. The primary transmutation product, Re, is produced through the β-decay of isotopes produced as a result of neutron absorbtion reactions in W, whereas Ta is produced in smaller quantities via the β-decay of isotopes produced by neutron loss reactions. Os on the other hand is produced via the absorbtion and β-decay of isotopes of Re [4, 7–9, 17]. Nuclear inventory modelling of transmutation in a concept DEMO reactor indicate that as much as 3.8 at.%Re and 1.38 at.%Os could be produced under a spatially homogeneous neutron spectrum, after 5 full-power years and a dose of 25 dpa [4]. By contrast, irradiation under fission spectrum neutrons in a reactor such as the High Flux Reactor (HFR) leads to post-irradiation composition an order of magnitude higher after the same irradiation time [4, 18].

Due to the current absence of a source of fusion spectrum neutrons, experimental studies of transmutation and radiation damage effects in W must rely on simulated neutron damage, either produced using fission spectrum neutrons or using charged particle irradiation (such as with heavy-ions or protons). Neutron irradiation in a number of test reactors have found radiation induced hardening [15, 19–21] and a decrease in the ductile to brittle transition temperature [14, 16]. At low doses the hardening response is primarily a result of an accumulation of dislocation loops and
small vacancy clusters, whereas at higher doses the response is dominated by voids and brittle intermetallic precipitates [15].

Precipitation of non-equilibrium phases in irradiated material is a commonly observed phenomena and is a result of radiation induced precipitation (RIP), whereas radiation enhanced diffusion (RED) can lead to a faster rate of formation of equilibrium phases. Due to their low solubility limits, Re and Os precipitation is therefore expected under irradiation, but is observed to occur at concentrations well below said solubility limits. Neutron irradiation studies have found a combination of σ and χ phase intermetallics, nanoscale Re and Os rich clusters, as well as voids decorated with Re and Os [11, 21–27]. Similarly, heavy-ion irradiation studies have found Re and Os precipitation, and have shown that Os can act to suppress Re clustering [10, 28].

The unique behaviour of W-(Re,Os) alloys suggests that increased solute transport may be a result of strong Re and Os binding to point defects. The most stable self-interstitial configuration in pure W is the 〈 111 〉 dumbbell, which undergoes fast 1D migration [29, 30]. Strong binding between Re [31] and Os atoms with interstitials leads to the formation of thermodynamically stable W-solute mixed interstitial dumbbells. Both Re and Os mixed interstitial dumbbells result in enhanced solute transport thanks to their combined high stability and mobility; the latter through effective 3D migration without dissociation via a series of translation and rotation processes [32–34]. Kinetic Monte Carlo (KMC) models in which this transport process has been implemented have predicted the formation of Re rich precipitates, facilitated primarily by W-Re mixed interstitial migration [35, 36]. However, the W-Re energy model does not capture the expected vacancy clustering, because the pair interactions are independent of the local vacancy concentration.

The di-vacancy interaction in W is repulsive, whereas a 3-vacancy cluster and larger is attractive [37], which leads to an under prediction of void formation in a model in which the pair interactions are independent of local concentration. Both void formation and decoration with Re has been predicted using a cluster expansion Hamiltonian [37]. But as has been outlined in experimental studies, the W-Re binary system alone cannot be used to model transmutation effects in a fusion reactor, because of the strong effect of Os on the resultant microstructure. Modelling radiation damage in W therefore relies on the development of a new energy model for W, Re and Os atoms, the derivation of which is one of the key objectives of this study.

This study features a new broken-bond energy model for KMC simulations in the W-Re-Os alloy, where both vacancy and interstitial dumbbell type defects are considered. The inclusion of Os interactions allows for more representative simulations of W component evolution in a fusion reactor. This new model captures the attractive binding of clusters of 4 or more vacancies, as well as the repulsion between smaller clusters, including the di-vacancy. Density functional theory (DFT) calculations were performed over a wide range of configurations for the W-Re-Os-Vacancy quaternary system and were used to parameterise the various pair interactions present within the new model. The attractive nature of the vacancies in clusters containing 4 or
A Solute and Defect Concentration Dependant Ising Model for W-Re-Os 4

more vacancies is accounted for by considering the local vacancy concentration in the calculation of the atom-vacancy pair interaction coefficients. Similarly, we consider the local solute concentration when calculating the interactions between Re and or Os solute and W atoms. Metropolis Monte Carlo simulations were performed to investigate the equilibrium configurations of simulation cells containing vacancies. Canonical simulations were also performed to study the interaction between clusters of vacancies and Re/Os solute atoms.

## 2. Fitting Enthalpies of Mixing

### 2.1. Pair Interaction Energy Model

The model presented in this study is based on a pair interaction Hamiltonian, in which the interactions between pairs of individual atoms and point defects are considered. The Hamiltonian for a basic pair interaction model can be expressed as a sum of interaction coefficients between pairs of atoms, $\epsilon$, according to equation 1, where (m) is the nearest neighbour shell, $\alpha$ and $\beta$ are the atom types of the two sites being considered and n is the number of bonds between $\alpha$ and $\beta$.

$$ \mathcal{H} = \sum_{m} \sum_{\alpha, \beta} n_{\alpha-\beta}^{(m)} \epsilon_{\alpha-\beta}^{(m)} \quad (1) $$

In this study, the description given by equation 1 is expanded such that the pair interaction coefficients, $\epsilon$, depend on their local solute or vacancy concentration. The general Hamiltonian in the model presented in this study is therefore given by equation 2,

$$ \mathcal{H} = \sum_{i,j,\alpha,\beta} n_{\alpha,\beta}^{i,j}(x_{\alpha}^{i}, x_{\beta}^{j}) \epsilon_{\alpha,\beta}^{i,j}(x_{\alpha}^{i}, x_{\beta}^{j}) \quad (2) $$

where $n_{\alpha,\beta}^{i,j}$, is the number of bonds between sites i and j, between atoms of type $\alpha$ and $\beta$; $\epsilon_{\alpha,\beta}^{i,j}$, is the bond energy between atoms of type $\alpha$ and $\beta$; and $x_{\alpha}$ is the local concentration of $\alpha$. Interactions are considered between the first nearest neighbour (1NN) and second nearest neighbour atoms (2NN), and the value of the pair interaction, $\epsilon_{\alpha,\beta}^{i,j}$, depends on the distance between the sites, i and j being considered. For a bcc system such as W, the total number of atoms in the first and second nearest neighbour shells are $z_{T}^{(1)} = 8$, and $z_{T}^{(2)} = 6$. The local concentration of $\alpha$, $x_{\alpha}$, is calculated to the 2NN shell of the two sites being considered.

### 2.2. DFT Calculations

The DFT calculations used for fitting enthalpy of mixing [37, 38] in this study were performed with the Vienna Ab-initio Simulation Package (VASP) [39–41]. The interaction among ions and electrons is described by using the projector augmented waves (PAW) method [42]. Exchange and correlation were treated in the generalised
A Solute and Defect Concentration Dependant Ising Model for W-Re-Os

5

gradient approximation GGA-PBE [43], with PAW potentials for W, Re and Os containing semi-core p electron contributions. To treat clusters containing from 2 to 47 sites with solute atoms and vacancies, orthogonal bcc super-cells containing 128 and 250 sites were used in calculations. Total energies were calculated using the Monkhorst-Pack mesh [44] of k-points in the Brillouin zone, with the k-mesh spacing of 0.15 Å⁻¹. This corresponds to 4×4×4 or 3×3×3 k-point meshes for a bcc super-cell of 4×4×4 or 5×5×5 bcc structural units, respectively. The plane wave cut-off energy was 400 eV. The total energy convergence criterion was set to 10⁻⁶ eV/cell, and force components were relaxed to 10⁻³ eV/Å. Super-cell calculations were performed considering vacancy clusters interacting with solute atoms (Re, Os) in bcc W lattice under constant pressure conditions, with structures optimised by relaxing both atomic positions as well as the shape and volume of the super-cell. Besides the data obtained from the large supercell structures for random configurations, the DFT calculations were also performed for 2×2×2 supercell structures generated from decoration of one vacancy and other W, Re, Os atoms from 58 bcc-like ordered structures originally used for binary system [45, 46]. In addition, special quasi-random structures (SQS), created from Alloy Theoretic Automated Toolkit (ATAT) code [37, 38, 46, 47], for configurations with one vacancy and different W and solute atoms were also employed. More detailed information about the DFT energy data used in figures 1, 2, 3, 4 and 6 of the present work can be found from recent publications [13, 37, 38, 48–50].

DFT calculations also enable determining the binding energy of defect clusters. The binding energy of each defect configuration is computed using the following definition:

$$E_b(A_1, A_2, \dots, A_N) = \sum_{i}^{N} E(A_i) - [E(\sum_{i}^{N} A_i) - (N-1)E_{ref}] \quad (3)$$

where $A_i$ marks each component of the N-component cluster. $E(\sum_{i}^{N} A_i)$, $E(A_i)$, and $E_{ref}$ are respectively the energy of the supercell containing the cluster, that containing the isolated component $A_i$, and the perfect W atom super-cell. In this convention, positive energies stand for attractive interactions. For validation sake, an independent VASP study of binding energies for selected defective configurations of vacancy-(Re)-(Os) clusters in W has also been performed [51] and these DFT results are shown in figure 8. Here, each configuration with 1 to 4 vacancies and a variable amount of Re and Os atoms (up to 4 components in total) was embedded in the centre of a 5×5×5 super-cell with a bcc crystal structure and periodic boundary conditions on all sides. For each cluster composition, the configurations were picked semi-randomly to favour the most compact ones, and the maximum distance between one cluster component and the others could not exceed $\sqrt{3}a_0$. With respect to the previous DFT datasets, the latter calculations were performed with a plane-wave cutoff of 350 eV and a 3x3x3 k-point grid. More importantly, in place of zero-pressure conditions, the cell volume and shape were kept fixed to the obtained equilibrium lattice parameter of $a_0 = 3.1725$ Å.
<table><caption>**Table 1.** Local solute and defect dependencies of the pair interactions parametrised in this study, for each of the fitted binary systems.</caption><thead><tr><th>Interaction</th><th>Local Concentration Dependence</th></tr></thead><tbody><tr><td>W-Re</td><td>Re</td></tr><tr><td>W-Os</td><td>Os</td></tr><tr><td>Re-Os</td><td>Os</td></tr><tr><td>W-Vacancy</td><td>Vacancy</td></tr><tr><td>Re-Vacancy</td><td>Vacancy</td></tr><tr><td>Os-Vacancy</td><td>Vacancy</td></tr></tbody></table>

## 2.3. Overview of Fitting Methodology

In this study, we use the enthalpy of mixing, $\Delta H_{\text{mix}}$, calculated using DFT, to parameterise the various pair interactions. The enthalpy of mixing for a 4 component alloy containing constituents A, B, C and D, is given by equation 4, where $\Omega$ is the heat of solution and $x_A$ is the concentration of element A.

$$ \Delta H_{\text{mix}} = \sum_{\text{pairs}} \Omega^p x_A^p x_B^p + \sum_{\text{ternary}} \Omega^t x_A^t x_B^t x_C^t + \sum_{\text{quaternary}} \Omega^q x_A^q x_B^q x_C^q x_D^q + \dots $$

To parameterise the interactions, we assume that the contribution to $\Delta H_{\text{mix}}$ from the ternary and quaternary interactions given in equation 4, are negligible and can therefore be ignored. Therefore, $\Delta H_{\text{mix}}$ can be approximated as the sum of the 6 binary systems: W-Re, W-Os, Re-Os, W-Vacancy, Re-Vacancy and Os-Vacancy, as shown in equation 4.

$$ \Delta H_{\text{mix}} \approx \sum_{A,B} \Omega_s^{AB} x_A x_B \quad (4) $$

where $\Omega_s$ is a function of the pair interaction coefficients. After parameterising the model under equation 4, the validity of this assumption is assessed by comparing with the full set of DFT calculations for the quaternary and ternary systems.

For each of the binary systems, the convex hull of the $\Delta H_{\text{mix}}$ data from DFT is fit with a polynomial function. By comparing the fitted function for $\Delta H_{\text{mix}}$ with equation 4, the pair interactions can be obtained as a function of the solute atom concentration. In the following sections, the details of this fitting method are outlined, first for the interatomic interactions in section 2.4, and then for the vacancy-atom interactions in section 2.5.

## 2.4. Interatomic Interactions

The interactions between two identical elements (ie AA, BB, or CC) can be obtained from the cohesive energies for each system, which were calculated using DFT and are
<table><caption>**Table 2.** List of DFT energies used to fit the atom-atom interactions for the W-Re and W-Os. Also listed are the binding energies (positive energies indicate an attraction), $E_{V-V}^b$ for the di-vacancy from DFT calculations at first and second nearest neighbour positions.</caption><thead><tr><th>Configuration</th><th>Description</th><th>DFT Value (eV)</th><th>Ref</th></tr></thead><tbody><tr><td>E<sub>W</sub><sup>coh</sup></td><td>Cohesive Energy</td><td>-8.4117</td><td>[37]</td></tr><tr><td>E<sub>Re</sub><sup>coh</sup></td><td>Cohesive Energy</td><td>-7.4935</td><td>[37]</td></tr><tr><td>E<sub>Os</sub><sup>coh</sup></td><td>Cohesive Energy</td><td>-7.4283</td><td>This Study</td></tr><tr><td>E<sub>V</sub><sup>f</sup></td><td>Formation Energy</td><td>-3.6190</td><td>[37]</td></tr><tr><td>E<sub>V-V</sub><sup>b(1)</sup></td><td>Binding Energy</td><td>-0.048</td><td>[49]</td></tr><tr><td>E<sub>V-V</sub><sup>b(2)</sup></td><td>Binding Energy</td><td>0.2860</td><td>[49]</td></tr></tbody></table>

summarised in table 2. The cohesive energy of element A, $E_{A,coh}$, can be approximated by the sum of the AA pair interactions, $\epsilon_{AA}^{1,2}$ over $z_1 = 8$ 1NN and $z_2 = 6$ 2NN as shown in equation 5.

$$E_{A,coh} = -\frac{z_1}{2}\epsilon_{AA}^{1} - \frac{z_2}{2}\epsilon_{AA}^{2} \qquad (5)$$

To obtain the values of the pair interactions from equation 5, a scaling relation between the first and second nearest neighbour interactions is used, as shown in equation 6. For the bcc system, the scaling factor has a value of $\theta = 0.4219$ [35].

$$\frac{\epsilon_{AA}^2}{\epsilon_{AA}^1} = \frac{\epsilon_{AB}^2}{\epsilon_{AB}^1} = \theta \qquad (6)$$

The $\Delta H_{\text{mix}}$ data from DFT calculations is used to parameterise the cross interactions between pairs of dissimilar atoms, following a similar method to that used by Huang et. al [35]. Figure 1 shows the enthalpy of mixing data for the three binary atomic systems considered in this study: W-Re, W-Os and Re-Os. The data shown in figure 2 consists of configurations containing both random and ordered structures [45]. In the current approach, the minima of these data are fit using a mean field approximation, which may result in overestimation of the stability of some of these ordered configurations.

Fitting of the data was performed in Python using the ConvexHull and curve_fit tools built into the open source scipy library [52, 53]. The DFT data is shown in figure 1 as hollow data points, whereas the extracted convex hull is plotted using filled markers.

As shown in equation 4, $\Delta H_{\text{mix}}$ for a binary system can be expressed in terms of the heat of solution, $\Omega_s$, and the concentration $x_A$ and $x_B$ of A and B respectively. For a binary system, $x_A + x_B = 1$, therefore equation 4 can be written in terms of the concentration of solute only, as shown in equation 7,

$$\Delta H_{\text{mix}} = x_B(1 - x_B)\Omega_s \qquad (7)$$

where $x$ is the concentration of solute atoms. The heat of solution can be expressed in terms of the 1NN and 2NN contributions according to equation 8, where $z_1$ and $z_2$
A Solute and Defect Concentration Dependant Ising Model for W-Re-Os

8

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Plots showing the binary DFT data for the (a) W-Re, (b) W-Os and (c) Re-Os systems, with the DFT data shown as hollow points and the convex hull with filled square markers. Also shown is the function fitted to the convex hull of the data, with the corresponding fitting parameters (see equations 12 and 13) given in the legend.</figcaption></figure>

are the number of 1NN and 2NN and $\Omega_s^{(1)}$ and $\Omega_s^{(2)}$ are the respective heats of mixing

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
A Solute and Defect Concentration Dependant Ising Model for W-Re-Os
9
for the first and second nearest neighbours.

$$
\Omega_s = z_1 \Omega_s^{(1)} + z_2 \Omega_s^{(2)} \tag{8}
$$

They can be written as a function of the pair interactions, $\epsilon$, as shown in equations 9 and 10, where $\epsilon_{A-B}^{(1)}$ is the cross-atom pair interaction between two atoms of type A and B. The pair interactions, $\epsilon_{A-A}^{(1)}$ and $\epsilon_{B-B}^{(1)}$ are obtained from the cohesive energy for each system using equations 5 and 6.

$$
\Omega_s^{(1)} = \epsilon_{A-B}^{(1)} - \frac{1}{2} (\epsilon_{A-A}^{(1)} + \epsilon_{B-B}^{(1)}) \quad (9)
$$

$$
\Omega_s^{(2)} = \epsilon_{A-B}^{(2)} - \frac{1}{2} (\epsilon_{A-A}^{(2)} + \epsilon_{B-B}^{(2)}) \quad (10)
$$

Combining equation 8 with equations 9 and 10 and using the relation between the first and second nearest neighbours given in equation 6, $\Omega_s$ can be written as a function of the 1NN interaction between the cross A-B interaction, as shown in equation 11.

$$
\Omega_s = \epsilon_{A-B}^{(1)}(z_1 + \theta z_2) - \frac{z_1}{2} (\epsilon_{A-A}^{(1)} + \epsilon_{B-B}^{(1)}) - \frac{z_2}{2} (\epsilon_{A-A}^{(2)} + \epsilon_{B-B}^{(2)}) \quad (11)
$$

The convex hull of the DFT data was fit with a custom polynomial function, $f(x) = x(1-x)(\gamma x + \delta)$, using a least squares method employing the Levenberg-Marquardt algorithm. The optimised fits for the interatomic systems are shown in figure 1 along with the fitting constants which are provided in the legend. The lowest order polynomial function that provided an acceptable visual fit to the data shown in figure 1 was chosen. As the $\epsilon_{A-A}^{(1)}$ interactions are constant with respect to local solute and defect concentrations, we assume that the $\epsilon_{A-B}^{(1)}$ interactions have the same form as $\Omega_s$, such that $\epsilon_{A-B}^{(1)} = cx + d$, where c and d are unknown quantities. Comparing the optimised fitting variables, $\delta$ and $\gamma$, with the $\Omega_s$ given by equation 11, the following definitions for the unknown variables c and d can be made:

$$
d = \frac{\delta + \frac{z_1}{2} (\varepsilon_{A-A}^{(1)} + \varepsilon_{B-B}^{(1)}) + \frac{z_2}{2} (\varepsilon_{A-A}^{(2)} + \varepsilon_{B-B}^{(2)})}{(z_1 + \theta z_2)} \quad (12)
$$

$$
c = \frac{\gamma}{(z_1 + \theta z_2)} \qquad (13)
$$

Using equations 12 and 13 and the fitting variables $\gamma$ and $\delta$ obtained from the fitting of the data shown in figure 1 for the W-Re, W-Os and Re-Os systems, the values of the pair interaction coefficients $\epsilon_{W-Re}^{(1,2)}$, $\epsilon_{W-Os}^{(1,2)}$ and $\epsilon_{Re-Os}^{(1,2)}$ were calculated as a function of solute concentration and are given in table 3. The solute concentration dependencies of these interactions are summarised in table 1.
A Solute and Defect Concentration Dependant Ising Model for W-Re-Os

10

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Plots showing the binary DFT data for the (a) W-Vac, (b) Re-Vac and (c) Os-Vac systems, with the DFT data shown as hollow points and the convex hull with filled square markers. Also shown is the function fitted to the convex hull of the data, with the corresponding parameters given in the legend. In (a) where more limited data was available, the convex hull used for fitting consists of all of the hollow DFT points.</figcaption></figure>

## 2.5. Defect Interactions

A similar approach to that of section 2.4 was used to obtain the defect-atom interaction coefficients. The enthalpy of mixing for the W-vacancy, Re-vacancy and
A Solute and Defect Concentration Dependant Ising Model for W-Re-Os 11

Os-vacancy systems was calculated using DFT simulations for the full range of vacancy concentrations, constrained to a bcc lattice, as shown in figure 2. Compared to the interatomic interactions discussed in the previous section, higher order polynomials were required to generate an acceptable fit for the data for each of the atom-vacancy binary systems. Similarly to equation 11, the heat of mixing, $\Omega_s$, for a binary atom-vacancy system is given by equation 14, where A is the atom type and $\epsilon_{V-V}^{(i)}$ is the pair interaction between two vacancies at the ith nearest neighbour distance.

$$ \Omega(x) = \epsilon(x)_{A-V}^{(1)}(z_1 + \theta z_2) - \frac{z_1}{2} (\epsilon_{A-A}^{(1)} + \epsilon_{V-V}^{(1)}) - \frac{z_2}{2} (\epsilon_{A-A}^{(2)} + \epsilon_{V-V}^{(2)}) \quad (14) $$

The vacancy-vacancy pair interaction coefficient, $\epsilon_{V-V}^{(i)}$, is set to zero as there is no vacancy-vacancy interaction within the voids. Interactions occur between the atoms on the periphery of the void with the surrounding atoms in the matrix. The reference energy for the vacancies in the DFT simulations is also taken as zero, which is the equivalent of the cohesive energy for the pure atomic systems. By setting $\epsilon_{V-V}^{(i)} = 0$, equation 14 can be reduced to the following:

$$ \Omega(x) = \epsilon(x)_{A-V}^{(1)}(z_1 + \theta z_2) - \frac{\epsilon_{A-A}^{(1)}}{2}(z_1 + \theta z_2) \quad (15) $$

The function, $g(x)$ fitted to the atom-vacancy interactions is given by equation 16, where $\delta$, $\gamma$, $\beta$ and $\alpha$ are the variables used to optimise the fit.

$$ g(x) = x(1-x)(\delta + \gamma x + \beta x^2 + \alpha x^3) \quad (16) $$

As for the interatomic interactions, we assume that the atom-defect pair interactions, $\epsilon_{A-V}^{(1)}$, follow the same form as the fitted function in equation 16: ($\epsilon_{A-V}^{(1)} = d + cx + bx^2 + ax^3$), where a, b, c and d are the pair interaction components that need to be calculated.

The results of the fitting are shown in figure 2(a), (b) and (c) for the W-vacancy, Re-vacancy and Os-vacancy systems respectively. The DFT data is shown using hollow markers whereas the points used for fitting are shown as filled square markers. For the W-vacancy case, all of the DFT points plotted in figure 2(a) were used to optimise the fit. The convex hull of the data was calculated using the ConvexHull function from the scipy library for vacancy concentrations lower than 0.1. Data points at higher vacancy concentrations were chosen manually to constrain the fitted function and optimise the fit in the data below a vacancy concentration of 0.1.

Comparing the variables optimising during fitting with the heat of solution given by equation 15 allows for the following expressions to be defined for the unknown variables in the atom-vacancy pair interaction coefficient:

$$ a = \frac{\alpha}{(z_1 + \theta z_2)} \quad (17) $$

$$ b = \frac{\beta}{(z_1 + \theta z_2)} \quad (18) $$
<table><caption>**Table 3.** Pair interactions calculated in this study where the pair interaction between atoms of type A and B is given by $\epsilon_{A-B}^{(1)} = ax_B^3 + bx_B^2 + cx_B + d$.</caption><thead><tr><th>Quantity</th><th>Fit based on:</th><th>a</th><th>b</th><th>c</th><th>d</th></tr></thead><tbody><tr><td>ε<sub>W-W</sub><sup>(1)</sup></td><td>Cohesive Energy</td><td>0</td><td>0</td><td>0</td><td>-1.5975</td></tr><tr><td>ε<sub>W-W</sub><sup>(2)</sup></td><td>Cohesive Energy</td><td>0</td><td>0</td><td>0</td><td>-0.6672</td></tr><tr><td>ε<sub>Re-Re</sub><sup>(1)</sup></td><td>Cohesive Energy</td><td>0</td><td>0</td><td>0</td><td>-1.4231</td></tr><tr><td>ε<sub>Re-Re</sub><sup>(2)</sup></td><td>Cohesive Energy</td><td>0</td><td>0</td><td>0</td><td>-0.6004</td></tr><tr><td>ε<sub>Os-Os</sub><sup>(1)</sup></td><td>Cohesive Energy</td><td>0</td><td>0</td><td>0</td><td>-1.4107</td></tr><tr><td>ε<sub>Os-Os</sub><sup>(2)</sup></td><td>Cohesive Energy</td><td>0</td><td>0</td><td>0</td><td>-0.5935</td></tr><tr><td>ε<sub>W-Re</sub><sup>(1)</sup></td><td>ΔH<sub>mix</sub>(W-Re)</td><td>0</td><td>0</td><td>-0.0112</td><td>-1.5266</td></tr><tr><td>ε<sub>W-Re</sub><sup>(2)</sup></td><td>ΔH<sub>mix</sub>(W-Re)</td><td>0</td><td>0</td><td>-0.0047</td><td>-0.6440</td></tr><tr><td>ε<sub>W-Os</sub><sup>(1)</sup></td><td>ΔH<sub>mix</sub>(W-Os)</td><td>0</td><td>0</td><td>-0.0662</td><td>-1.5221</td></tr><tr><td>ε<sub>W-Os</sub><sup>(2)</sup></td><td>ΔH<sub>mix</sub>(W-Os)</td><td>0</td><td>0</td><td>-0.0279</td><td>-0.6422</td></tr><tr><td>ε<sub>Re-Os</sub><sup>(1)</sup></td><td>ΔH<sub>mix</sub>(Re-Os)</td><td>0</td><td>0</td><td>0.0007</td><td>-1.4359</td></tr><tr><td>ε<sub>Re-Os</sub><sup>(2)</sup></td><td>ΔH<sub>mix</sub>(Re-Os)</td><td>0</td><td>0</td><td>0.0003</td><td>-0.6058</td></tr><tr><td>ε<sub>W-V</sub><sup>(1)</sup></td><td>ΔH<sub>mix</sub>(W-Vac)</td><td>-0.0143</td><td>-0.9539</td><td>0.9255</td><td>-0.5472</td></tr><tr><td>ε<sub>W-V</sub><sup>(2)</sup></td><td>ΔH<sub>mix</sub>(W-Vac)</td><td>-0.0060</td><td>-0.4024</td><td>0.3905</td><td>-0.2309</td></tr><tr><td>ε<sub>Re-V</sub><sup>(1)</sup></td><td>ΔH<sub>mix</sub>(Re-Vac)</td><td>2.0347</td><td>-4.2840</td><td>2.7552</td><td>-0.9242</td></tr><tr><td>ε<sub>Re-V</sub><sup>(2)</sup></td><td>ΔH<sub>mix</sub>(Re-Vac)</td><td>0.8584</td><td>-1.8073</td><td>1.1624</td><td>-0.3899</td></tr><tr><td>ε<sub>Os-V</sub><sup>(1)</sup></td><td>ΔH<sub>mix</sub>(Os-Vac)</td><td>0</td><td>-1.9183</td><td>2.6751</td><td>-1.3951</td></tr><tr><td>ε<sub>Os-V</sub><sup>(2)</sup></td><td>ΔH<sub>mix</sub>(Os-Vac)</td><td>0</td><td>-0.8093</td><td>1.1286</td><td>-0.5886</td></tr><tr><td>ε<sub>V-V</sub><sup>(1)</sup></td><td>Reference Energy</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>ε<sub>V-V</sub><sup>(2)</sup></td><td>Reference Energy</td><td>0</td><td>0</td><td>0</td><td>0</td></tr></tbody></table>

$$c = \frac{\gamma}{(z_1 + \theta z_2)} \qquad (19)$$

$$d = \frac{\delta + \frac{z_1}{2} \epsilon_{A-A}^{(1)} + \frac{z_2}{2} \epsilon_{A-A}^{(2)}}{(z_1 + \theta z_2)} \qquad (20)$$

Using equations 17-20 and the optimised fitting variables $\alpha$, $\beta$, $\gamma$ and $\delta$, summarised in figure 2, the values of the pair interaction coefficients $\epsilon_{W-Vac}^{(1,2)}$, $\epsilon_{Re-Vac}^{(1,2)}$ and $\epsilon_{Os-Vac}^{(1,2)}$ were calculated and are summarised in table 3. For the Os-vacancy system a 2nd order polynomial is sufficient; therefore, $a = 0$ for this system.

## 2.6. Evaluating Fit Accuracy

The accuracy of the fitted polynomial functions shown in figures 1. and 2 were assessed by calculating the loss function $\bar{l}$ for each of the binary data sets, as well as for the corresponding ternary systems shown in figure 3 and the quaternary W-Re-Os-Vacancy system. The value of $\bar{l}$ is calculated for each data set by considering the variation between the $\Delta H_{mix}$ predicted by the fitted function, $\Delta H_{mix}^{pred}$ and the minimum value of the DFT data at that point. The DFT data was first binned into $i$ bins, ($i = 30$ for
A Solute and Defect Concentration Dependant Ising Model for W-Re-Os 13

the data shown in figures 1. and 2) and the minimum $\Delta H_{\text{mix}}$, $\Delta H_{\text{mix},i}^{\text{DFT}}$, from that bin, $i$, was obtained. The corresponding value of $\Delta H_{\text{mix}}$ predicted by the pair interaction model for bin $i$, is also calculated and is given by $\Delta H_{\text{mix},i}^{\text{pred}}$. The value of $\ell(i)$ for each of the $i$ bins is calculated. $\bar{\ell}$ is calculated according to equation 21,

$$ \bar{\ell} = \frac{1}{N_i} \sum_i |\Delta H_{\text{mix},i}^{\text{pred}} - \Delta H_{\text{mix},i}^{\text{DFT}}| \quad (21) $$

where $N_i$ is the number of bins for which the difference between DFT and the predicted $\Delta H_{\text{mix}}$ was calculated (i.e. the number of bins which contained DFT data). Overall the fitting showed good agreement with the DFT data, as summarised in table 4, with the highest deviation present for the solute-vacancy systems. In this case, the disagreement between DFT and the $\Delta H_{\text{mix}}$ calculated using the pair interaction model comes from fitting the interactions over the whole vacancy concentration range, whereas at lower vacancy concentrations, which are the most physically relevant for this study, the agreement is better.

As discussed in section 2.3, we assume that $\Delta H_{\text{mix}}$ could be fitted by neglecting the ternary and quaternary interactions. Therefore, in sections 2.4 and 2.5, we fit $\Delta H_{\text{mix}}$ for each of the 6 binary systems, assuming that $\Delta H_{\text{mix}}$ can be approximated by equation 4. The validity of this assumption is assessed by comparing the ternary and quaternary $\Delta H_{\text{mix}}$ produced by equation 4, to the ternary and quaternary DFT data for the equivalent concentration.

For a ternary system of A, B and C atoms types, $\Delta H_{\text{mix}}$ can be written as a function of the concentrations, $x_A$, $x_B$ and $x_C$ and the heat of solution, $\Omega_s$, for each binary system, as shown in equation 22.

$$ \Delta H_{\text{mix}} \approx \Omega_s^{AB}(1-x_B)x_B + \Omega_s^{AC}(1-x_C)x_C + \Omega_s^{BC}x_Bx_C \quad (22) $$

Figure 3 shows plots of $\Delta H_{\text{mix}}$, given by equation 22, using the calculated pair interactions summarised in table 3 for the W-Re-Os, W-Re-Vacancy, W-Os-Vacancy and Re-Os-Vacancy systems. Multiple stacked points show the locations of DFT simulations for a given solute concentration. The DFT data was binned into concentration bins, 0.01 in width, from which the minimum $\Delta H_{\text{mix}}$ value was obtained. The difference between $\Delta H_{\text{mix}}$ produced by the pair interaction model, $\Delta H_{\text{mix}}^{\text{pred}}$, and the value calculated using DFT, $\Delta H_{\text{mix}}^{\text{DFT}}$, was derived using the loss function $\bar{\ell}$ given by equation 21. Figure 3 shows the full concentration range of both solute atoms and vacancies, with a vacancy concentration of 0.1 marked by the red line in (b), (c) and (d). Higher vacancy concentrations are not observed experimentally, and the vacancy concentration in irradiated materials is typically much lower than this. Therefore, the data shows at higher vacancy concentrations has limited physical meaning, but is included here to illustrate the function that was fit to the DFT data.

The calculated values of $\bar{\ell}$ for the ternary systems shown in figure 3 are summarised in table 4, as well as the corresponding $\bar{\ell}$ calculated for the binary and quaternary
A Solute and Defect Concentration Dependant Ising Model for W-Re-Os 14

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Enthalpy of mixing, $\Delta H_{\text{mix}}$, of the ternary (a) W-Re-Os, (b) W-Re-Vac, (c) W-Os-Vac and (d) Re-Os-Vac systems, described by equation 4 as a function of concentration. Overlaid in black are the enthalpy of mixing points for various concentrations. The opacity of the DFT points is a function of the number of configurations at that concentration point. The red line in the systems containing vacancies shows a vacancy concentration of 0.1, beyond which there is limited physical meaning.</figcaption></figure>

systems. For the W-Re-Os system, the deviation from the estimated surface was small ($\bar{l} = 0.0250$eV). In the W-Re-Vacancy and W-Os-Vacancy systems the calculated $\bar{l}$ was larger but still relatively small at 0.0357eV and 0.0914eV respectively. The value of $\bar{l}$ was larger for the Re-Os-Vacancy system at $\bar{l} = 0.2000$eV and is likely a result of the imperfect fitting shown in figure 2. Overall, these results show that the $\Delta H_{\text{mix}}$ for the 4 ternary systems considered can be appropriately reconstructed using the calculated pair interaction model, without the need for additional concentration dependencies.
<table><caption>**Table 4.** Loss function, $\bar{l}$, calculated using equation 21, for each of the binary, ternary and quaternary systems studied in this project. $\bar{l}$ is calculated for the quaternary system, where $\Delta H_{\text{mix}}^{\text{pred}}$ was calculated using equation 4</caption><thead><tr><th>System</th><th>l̄ (eV)</th></tr></thead><tbody><tr><td>W-Re</td><td>0.0071</td></tr><tr><td>W-Os</td><td>0.0262</td></tr><tr><td>Re-Os</td><td>0.0146</td></tr><tr><td>W-Vacancy</td><td>0.0580</td></tr><tr><td>Re-Vacancy</td><td>0.0967</td></tr><tr><td>Os-Vacancy</td><td>0.2090</td></tr><tr><td>W-Re-Os</td><td>0.0250</td></tr><tr><td>W-Re-Vacancy</td><td>0.0357</td></tr><tr><td>W-Os-Vacancy</td><td>0.0914</td></tr><tr><td>Re-Os-Vacancy</td><td>0.2000</td></tr><tr><td>W-Re-Os-Vacancy</td><td>0.1376</td></tr></tbody></table>

### 3. Monte Carlo Simulations

The parameterisation presented in the previous section was implemented into a modified version of the open source Stochastic Parallel PARticle Kinetic Monte Carlo Simulator (SPPARKS) code, discussed in previous publications [36]. In this model, process rates are calculated using the change in energy, $\Delta E_{ij}$, between the initial and final states, $i$ and $j$, respectively. Each is calculated using the sum of pair interaction coefficients between atoms mapped to a rigid bcc lattice.

Metropolis Monte-Carlo (MMC) simulations were performed in which vacancy type defects exchange positions with a neighbouring atom. The acceptance probability for a metropolis exchange is calculated according to equation 23

$$ \prod_{i \to j} = \min \left\{ 1, \exp \left( -\frac{\Delta \mathcal{H}_{i \to j}}{k_B T} \right) \right\} \qquad (23) $$

where $\mathcal{H}_{i \to j}$ is the change that occurs in an exchange event between sites $i$ and $j$, at a temperature $T$. Metropolis simulations can be performed within this model using either interstitial or vacancy type defects, but in the present study only vacancy type defects are considered.

Canonical simulations were also performed in which atoms can exchange with any neighbouring site within the 1NN shell, defect or atom. This allows for simulations with a fixed solute and defect concentration to be carried out, and for the study of the equilibrium positions of solute atoms and equilibrium configurations of solute-defect clusters. Exchange events in the canonical simulations are also accepted with a probability given by equation 23, where $\Delta \mathcal{H}_{i \to j}$ is the change in energy between the initial and final states in an exchange.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Comparison of binding energies of different vacancy clusters up to 15 vacancies.</figcaption></figure>

Finally, Semi-Grand canonical Monte-Carlo (SGMC) simulations were carried out within the current model. In SGMC simulations, atoms are selected at random from the simulation cell and exchanged with a reservoir in which the atoms do not interact. Both solvent and solute atoms can undergo these exchange events if selected from the simulation cell, which results in a fluctuating solute concentration. The probability for an exchange event, $\prod_{i \to j}$, is calculated according to equation 24

$$ \prod_{i \to j} = \min \left\{ 1, \exp \left( \frac{-\Delta E_{i \to j} - \Delta \mu \Delta N_B}{k_B T} \right) \right\} \quad (24) $$

where $\Delta E_{i \to j}$ is the change in energy between initial state i and final state j; $\Delta \mu$ is the change in chemical potential, $\mu$, given by $\Delta \mu = \mu_A - \mu_B$ for elements A and B; $\Delta N_B$ is the change in the number of solute atoms, $N_B$; and T is the temperature.

In a SGMC simulation using a fixed temperature and chemical potential, the concentration of solute within the box fluctuates through swaps with the reservoir until the system reaches equilibrium. Inside a miscibility gap where a mixture of phases are present, the change in chemical potential does not determine the composition of the system. The system is therefore unstable within these regions and jumps to single phase regions, where the concentration is prescribed by the chemical potential. By plotting the change in chemical potential, $\Delta \mu$, against the equilibrium concentration of solute atoms, the miscibility gap can be identified.
A Solute and Defect Concentration Dependant Ising Model for W-Re-Os

17

<figure><img src="image_7.png" /><figcaption>Figure 5. Void nucleation under the new energy parametersiation, in a 64x64x64 W atom box, simulated using a metropolis algorithm, with number of MC steps and vacancy concentration highlighted. Individual vacancies are shown in blue and are introduced with a dose rate of 1 × 10⁻² per MC step.</figcaption></figure>

4. Results and Discussion

4.1. Binding Energies

The binding energies of various vacancy and solute-vacancy clusters were calculated and compared to DFT results, to assess the capabilities of the new parameterisation. The formation energy, $E_f$, of a vacancy-solute cluster is given by equation 25

$$E_f = E_{\text{box}} - n_A E_{\text{coh}}(A) - n_B E_{\text{coh}}(B) - n_C E_{\text{coh}}(C) \quad (25)$$

where $E_{\text{box}}$ is the total energy of a simulation cell of N atoms containing $n_A$ atoms of type A, $n_B$ B type solute atoms and $n_C$ atom of type C. $E_{\text{coh}}(A)$, $E_{\text{coh}}(B)$ and $E_{\text{coh}}(C)$ are the cohesive energies of elements A, B and C respectively.

Following the definition of Eq.3, for a cluster configuration $i$, in a box of element A, with $n$ vacancies, $m$ type B solute atoms and $k$ type C solute atoms, the cluster binding energy, $E_{\text{bind}}^{(i)}$ is given by equation 26 [37],

$$E_{\text{bind}}^{(i)} = \underbrace{nE_{\text{box}}(1\text{V}) + mE_{\text{box}}(1\text{B}) + kE_{\text{box}}(1\text{C})}_{\text{Single Solute/VacancyBoxEnergy}} - \underbrace{E_{\text{box}}^{(i)}(n\text{V}; m\text{B}; k\text{C})}_{\text{ClusteriBoxEnergy}} - \underbrace{(n + m + k - 1)E_{\text{box}}(A)}_{\text{PureABox}} \quad (26)$$
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Cluster binding energies for a single vacancy surrounded with Re and Os solute calculated using the new W-Re-Os-Vac energy model. Also shown by the hollow markers is the binding energy calculated for the equivalent configuration using DFT.</figcaption></figure>

where $E_{\text{box}}(\text{1B})$ is the total energy of a simulation cell containing a single atom of type B; $E_{\text{box}}(\text{A})$ is the energy of a simulation cell containing no solute atoms, and $E_{\text{box}}^{(i)}(\text{nv}; \text{mB}; \text{kC})$ is the energy of the box containing the cluster, $i$. Equation 26 can be used to calculate the binding energy of both vacancies or a mixture of vacancy and solute atoms. In this form, a positive value of the binding energy indicates an attractive interaction.

## 4.2. Vacancy Clustering

The binding energy of vacancy clusters, ranging from the di-vacancy to clusters of 15 vacancies, was evaluated using equation 26 (see figure 4). Void formation is predicted at clusters of 4 or more vacancies. There is a negative binding between pairs of vacancies in all of the configurations simulated, corresponding to a repulsive interaction.

In the clusters containing 4 vacancies, the majority of the configurations are also negative, indicating a repulsion, but in some configurations the binding energy is positive, indicating an attractive interaction in some circumstances. This trend continues in larger vacancy clusters, which show an increasing attractive binding energy as the size of the cluster grows. This suggests that under the new parameterisation, the formation of large voids is predicted, as has been observed experimentally under neutron irradiation. In previous pair interaction models in which the local vacancy concentration was not considered, this phenomena was not observed.
A Solute and Defect Concentration Dependant Ising Model for W-Re-Os

19

<figure><img src="image_9.png" /></figure>

**Figure 7.** Canonical simulations for a W box containing 0.5at.% vacancies, evolved with metropolis to form voids as shown in (a). Canonical simulations were carried out for (b) W-1.4Re (c) W-1.4Os and (d) W-1.4Re-0.1Os, randomly distributed through the initial simulation cell shown in (a) at a temperature of 1100K.

### 4.3. Void Nucleation

To demonstrate the formation of voids under the new parameterisation, MMC simulations were carried out at 1800 K. Figure 5 shows a 64x64x64 atom simulation cell into which vacancies were inserted at a rate of $1 \times 10^{-2}$ per MC step. The vacancies are shown in blue with W atoms hidden and the MC step and vacancy concentration is labelled under each image. Initially there is a random distribution of vacancies in the simulation cell, but after $2 \times 10^6$ MC steps, at a vacancy concentration of $1.3 \times 10^{-4}$, a small vacancy cluster starts to appear. Over time this cluster begin to grow into a large void, as shown at $3 \times 10^7$ MC steps at a vacancy concentration of $3.6 \times 10^{-3}$, together with the nucleation of other voids which also continue to grow.

### 4.4. Vacancy-Solute Clustering

To study the effect of Re and Os solute atoms on the stability of vacancy clusters, the binding energy of a single vacancy surrounded by an increasing number of solute atoms was calculated using equation 26. A central vacancy was decorated with an increasing number of Re or Os atoms within the 1NN shell, until the vacancy was fully surrounded by solute atoms. The calculated binding energies for the single vacancy-solute clusters are shown in figure 6.

Figure 6 shows an increasingly attractive binding energy as the number of Re and Os atoms surrounding the vacancy increases. The equivalent DFT data is also plotted in figure 6, showing a similarly increasing trend with increasing solute atoms. Os had a greater effect on the cluster stability than the Re, with the binding energy of the 8 atom cluster being more than double for Os than Re. Although trends are similar, the pair interaction model clearly overestimates the binding energy, most significantly for voids decorated with a high concentration of Os atoms.

The binding energies of a randomised set of cluster configurations were calculated for a database of Vacancy - (Re,Os) cluster configurations, according to equation 26.
<caption>A Solute and Defect Concentration Dependant Ising Model for W-Re-Os</caption>

20

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Vacancy solute cluster binding energies, calculated using new W-Re-Os energy parametrisation. Label indicates the composition of each cluster, with each point on the x-axis representing a different cluster composition. The first axis segment is comprised of vacancy clusters only. The second and third are vacancy-Re and vacancy-Os clusters respectively and the final axis segment is for mixed clusters of vacancy-Re-Os. The color represents the number of vacancies in a cluster, as shown in the legend.</figcaption></figure>
A Solute and Defect Concentration Dependant Ising Model for W-Re-Os 21

The calculated binding energies are shown in figure 8, with each point on the x-axis representing a different cluster composition, marked in the label for each x axis bin. The colour of each bin represents the number of vacancies in the cluster, as shown in the legend. The first axis segment shows only vacancy clusters, up to a maximum of 4. The following 2 axis segments are clusters of vacancies and Re and vacancies and Os respectively. The final axis segment shows a mixture of vacancies, Re and Os. All of the clusters considered here contain a maximum of 4 solute atoms and vacancies. Figure 8 shows clearly that both Re and Os stabilise small vacancy clusters including di-vacancies. Os has the strongest effect on stabilisation, with vacancy-Os clusters having the strongest attractive binding energies.

Figure 7 shows the results of 3 simulations in which a starting configuration of voids (with a vacancy concentration of 0.5at.%), formed using MMC, is used as a starting point for canonical simulations containing 1.4 at.%Re, 1.4 at.%Os and 1.4 at.%Re - 0.1 at.%Os (approximately the same composition found in neutron irradiation experiments [11]). The solute atoms were distributed randomly and then allowed to evolve using the canonical algorithm. The results of these simulations after a total of 3 × 10⁷ MC steps are shown in figure 7. There is a clear decoration of Re and Os visible in figures 7 (b) and (c), with Os decorating the voids most strongly. This was expected due to the stronger Os-Vacancy binding observed previously under the current parameterisation. When both Re and Os are present, Os segregates most strongly, followed by Re to a lesser extent. This is fully consistent with the stronger Os-vacancy than Re-vacancy binding shown in figure 8. This shows that there is a thermodynamic driving force for the segregation of these elements to voids in W. In order for void decoration to actually take place, however, a kinetic-transport mechanism needs to be active. This is most likely provided by the positive kinetic coupling existing between Re/Os solutes and point defects [51], and will be further addressed in future studies. In addition, further study using this new energy parameterisation within a kinetic MC model is needed to further understand the decoration of voids with Re and Os transmutation products.

4.5. *Phase Analysis*

To investigate the phases predicted in W-Re and W-Os binary systems under the new parameterisation, SGMC simulations were carried out. Solute atoms are exchanged with a reservoir with an acceptance probability determined by equation 24. The chemical potential for solute exchange, $\Delta\mu$, is incremented and the solute concentration after 100 MC steps/atom in the simulation cell, was measured. Plots of $\Delta\mu$ vs solute concentration are shown for the W-Re and W-Os systems in figures 9 (a) and (c) respectively. Figure 9 shows the presence of miscibility gaps in both Re and Os systems, indicating a phase change. The corresponding jump in the solute concentration indicates that a mixture of phases are present, and that the chemical potential for solute exchange does not determine the composition of the system, and is therefore unstable. No transformation is observed at higher temperatures, suggesting that longer time scales are required to reach
A Solute and Defect Concentration Dependant Ising Model for W-Re-Os 22

equilibrium at these temperatures or that the system is over the critical temperature.

For W-Re a jump to approximately 50 at.%Re is observed, particularly steep at lower temperatures, as shown in figure 9. This is consistent with the concentration of the σ phase intermetallic in W-Re. A similar transformation was observed in previous parameterisations of W-Re [35, 36]. For the W-Os system, two transformations are visible, one from 0-50 at.%Os and a second from 50-75 at.% Os.

The structure of the curves at low temperatures shown in figure 9 suggest that under the current parameterisation, several intermetallic phases are stable in both Re and Os. These phases have an ordered structure that is a result of fitting the ordered configurations shown in figure 2 with a mean field approach. Further improvement of the treatment of W-Re and W-Os binary interactions will be the subject of future work.

A similar approach was applied to the W-vacancy system to study vacancy precipitation. The equilibrium concentration of vacancies was however too low to be produced using a reasonable simulation box size. Instead, the vacancy precipitation is studied using a mean field approximation. The formation energy of the mono vacancy, $E_f$, can be calculated according to equation 25 for a simulation cell containing a single vacancy. The vacancy formation energy under the new parametrisation is 3.25 eV. The equilibrium vacancy concentration at temperature $T$ is then given by

$$C_{\text{eq}} = e^{-E_f/k_B T} \qquad (27)$$

Under the pair interaction model, the internal energy, $U$, can be expressed in terms of the heat of solution, $\Omega$, and the solute concentration, $x$, resulting the Helmholtz's free shown in equation 28.

$$\Delta F = x(1-x)\Omega(x) - T\Delta S \qquad (28)$$

The entropy of mixing contribution, $\Delta S$ is approximated as $k_B(x_1 \log(x_1) + x_2 \log(x_2))$, where $x_1$ is the concentration of solvent and $x_2$ is the concentration of solute. Substituting into equation 28, the Helmholtz free energy $\Delta F$ can be expressed in terms of temperature, $T$, and the solute concentration $x$, where $x_1 + x_2 = 1$, as shown in equation 29, following the Redlich-Kister expansion [54].

$$\Delta F = x(1-x)\Omega(x) - k_B T(x \log(x) + (1-x) \log(1-x)) \qquad (29)$$

The heat of solution, $\Omega(x)_s$ can be expressed in terms of the pair interactions according to equation 11. A plot of equation 29 against $x$ shows two minima, one at low solute concentrations and one at high solute concentrations. The common tangent between the minima was calculated at multiple temperatures and a phase diagram for the W-vacancy system was constructed, as shown in figure 10.

Figure 10 shows the W-vacancy phase diagram over the lower concentration region, on a normal and log scale respectively. Plotted in red is the equilibrium concentration of vacancies given by equation 27 as a function of temperature. The phase boundary determined by the common tangent method is higher in concentration than the equilibrium concentration across the whole range of temperature. This
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Plot of the change in chemical potential, $\Delta\mu$, for (a) Re atom exchange with a Re solute reservoir as a function of Re concentration, and (b) Os atom exchange with a Os solute reservoir as a function of Os concentration.</figcaption></figure>

implies that a certain concentration of vacancies is possible before void formation is predicted. The difference in the curves is due to the repulsive interaction between pairs of vacancies which prevents di-vacancy formation. The concentration at the boundary is very low, indicating that under typical irradiation conditions the material will be in a supersaturated state and will be characterised by the precipitation of vacancy clusters.

5. Conclusion

In this study, a new local solute and defect concentration dependent Ising model, for the W-Re-Os-Vacancy system was presented. The energy model was used to carry
A Solute and Defect Concentration Dependant Ising Model for W-Re-Os

24

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Low vacancy concentration section of the W-Vacancy phase diagram predicted by new parametrisation plotted on (a) a regular x axis and (b) on a log axis of concentration</figcaption></figure>

out metropolis, canonical and semi-grand canonical simulations within the W-Re-Os-Vacancy system, to assess the formation of vacancy and transmutation solute clusters in irradiated tungsten. The model presented in this study predicts an attractive binding energy in vacancy clusters containing 4 or more vacancies, whereas in smaller cluster a repulsive binding energy is observed. This represents an improvement on previous Ising models in which void formation is not predicted. The nucleation and growth of voids was demonstrated in MMC simulations, and decoration with Re and Os was observed in canonical simulations. Binding energy calculations showed that both Re and Os stabilise small vacancy clusters, with the effect of Os being the strongest. This implies that the presence of Re and Os transmutation products in a fusion reactor could stabilise small
## REFERENCES

25

voids. These findings agreed well with DFT calculations for equivalent configurations, however, the model overestimated the magnitude of the binding, particularly in the case with Os. The overestimation of the binding energies between vacancies and Re/Os solute atoms likely results in an overestimation of the void decoration in the MC simulations, however, this does not impact the relative strengths of the Os and Re interactions. The prediction of void formation within the new model allows for more accurate simulation of neutron damage under fusion spectrum neutrons, in which void formation has been observed in experiments from a range of different reactors. Further study using this new energy model, including KMC with Re, Os and vacancies will provide insight into the formation of decorated voids under neutron irradiation in W.

## Acknowledgments

This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement number 633053. The views expressed herein do not necessarily reflect those of the European Commission. MJL's work is also supported by the UK Engineering and Physical Sciences Research Council [EP/N509711/1] and the Culham Centre for Fusion Energy, United Kingdom Atomic Energy Authority through an Industrial CASE scholarship, [Project Reference Number 1802461]. DNM's work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 and 2019-2020 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. DNM also acknowledged funding from the RCUK [grant number EP/T012250/1]. LM and DNM also acknowledged the support from high-performing computing facility MARCONI (Bologna, Italy) provided by EUROfusion. EM would like to acknowledge start-up funds from Clemson University. MJL would also like to thank David Armstrong, Michael Moody and Paul Bagot for discussions of experimental work related to the study presented in this paper. Figures showing different simulation cells were created using the open-source visualisation software, OVITO (Open Visualisation Tool) [55].

## References

[1] D. E. J. Armstrong, P. D. Edmondson, and S. G. Roberts. "Effects of Sequential Tungsten and Helium Ion Implantation on Nano-Indentation Hardness of Tungsten". In: *Applied Physics Letters* 102.25 (2013), p. 251901.

[2] D. R. Lide. *CRC Handbook of Chemistry and Physics, 85th Edition*. Vol. 85. Boca Raton, FL: CRC press, 2005. ISBN: 0-8493-0485-7.

[3] M. R. Gilbert, M. Fleming, and J.-Ch. Sublet. "Automated Inventory and Material Science Scoping Calculations under Fission and Fusion Conditions". In: *Nuclear Engineering and Technology* 49.6 (2017), pp. 1346-1353.
[4] Mark R. Gilbert and Jean-Christophe Sublet. *Handbook of activation, transmutation, and radiation damage properties of the elements simulated using FISPACT-II and TENDL-2015; Magnetic Fusion Plants*. Version CCFE-R(16)36. UKAEA.

[5] T. Tanabe et al. "Temperature Dependence of Thermal Conductivity in W and W-Re Alloys from 300 to 1000 K". In: *Materials Letters* 57.19 (2003), pp. 2950-2953.

[6] M. J. Baldwin and R. P. Doerner. "Helium Induced Nanoscopic Morphology on Tungsten under Fusion Relevant Plasma Conditions". In: *Nuclear Fusion* 48.3 (2008), p. 035001.

[7] M. R. Gilbert et al. "Activation, Decay Heat, and Waste Classification Studies of the European DEMO Concept". In: *Nuclear Fusion* 57.4 (2017), p. 046015.

[8] MR Gilbert and J-Ch Sublet. "Neutron-induced transmutation effects in W and W-alloys in a fusion environment". In: *Nuclear Fusion* 51.4 (2011), p. 043005.

[9] M.R. Gilbert, J-Ch. Sublet, and S.L. Dudarev. "Spatial Heterogeneity of Tungsten Transmutation in a Fusion Device". In: *Nuclear Fusion* 57.4 (2017), p. 044002.

[10] A. Xu et al. "Ion-Irradiation-Induced Clustering in W-Re and W-Re-Os Alloys: A Comparative Study using Atom Probe Tomography and Nanoindentation Measurements". In: *Acta Materialia* 87 (2015), pp. 121-127.

[11] M. J. Lloyd et al. "Decoration of Voids with Rhenium and Osmium Transmutation Products in Neutron Irradiated Single Crystal Tungsten". In: *Scripta Materialia* 173 (2019), pp. 96-100.

[12] William D Klopp. *Review of ductilizing of group VIA elements by rhenium and other solutes*. Vol. 4955. National Aeronautics and Space Administration, 1968.

[13] F. Hofmann et al. "Non-Contact Measurement of Thermal Diffusivity in Ion-Implanted Nuclear Materials". In: *Scientific Reports* 5 (2015), p. 16042.

[14] R. G. Abernethy et al. "Effects of Neutron Irradiation on the Brittle to Ductile Transition in Single Crystal Tungsten". In: *Journal of Nuclear Materials* 527 (2019), p. 151799.

[15] X. Hu et al. "Defect Evolution in Single Crystalline Tungsten following Low Temperature and Low Dose Neutron Irradiation". In: *Journal of Nuclear Materials* 470 (2016), pp. 278-289.

[16] E. Gaganidze et al. "Fracture-mechanical properties of neutron irradiated ITER specification tungsten". In: *Journal of Nuclear Materials* 547 (2021), p. 152761.
ISSN: 0022-3115.

[17] C.-H. Huang, M. R. Gilbert, and J. Marian. "Simulating Irradiation Hardening in Tungsten under Fast Neutron Irradiation including Re Production by Transmutation". In: *Journal of Nuclear Materials* 499 (2018), pp. 204-215.
REFERENCES 27

[18] M. R. Gilbert and J-Ch. Sublet. *Handbook of Activation, Transmutation, and Radiation Damage Properties of the Elements Simulated using FISPACT-II and TENDL-2015; Nuclear Fission plants (HFR focus)*. Oxford: United Kingdom Atomic Energy Authority, 2015.

[19] A. Hasegawa et al. "Neutron Irradiation Effects on the Microstructural Development of Tungsten and Tungsten Alloys". In: *Journal of Nuclear Materials* 471 (2016), pp. 175-183.

[20] J. C. He et al. "Microstructural Development and Irradiation Hardening of W and W-(3-26) wt% Re Alloys after High-Temperature Neutron Irradiation to 0.15 dpa". In: *Nuclear fusion* 46.11 (2006), p. 877.

[21] Y Katoh et al. "Response of unalloyed tungsten to mixed spectrum neutrons". In: *Journal of Nuclear Materials* 520 (2019), pp. 193-207.

[22] M. Klimenkov et al. "Effect of Neutron Irradiation on the Microstructure of Tungsten". In: *Nuclear Materials and Energy* 9 (2016), pp. 480-483.

[23] M Dürrschnabel et al. "New insights into microstructure of neutron-irradiated tungsten". In: *Scientific Reports* 11.1 (2021), pp. 1-17.

[24] Philip D Edmondson, Baptiste Gault, and Mark R Gilbert. "An atom probe tomography and inventory calculation examination of second phase precipitates in neutron irradiated single crystal tungsten". In: *Nuclear Fusion* 60.12 (2020), p. 126013.

[25] Taehyun Hwang et al. "Effect of neutron irradiation on rhenium cluster formation in tungsten and tungsten-rhenium alloys". In: *Journal of Nuclear Materials* 507 (2018), pp. 78-86. ISSN: 0022-3115.

[26] Akira Hasegawa et al. "Neutron irradiation effects on the microstructural development of tungsten and tungsten alloys". In: *Journal of Nuclear Materials* 471 (2016), pp. 175-183.

[27] Xunxiang Hu et al. "Transmutation-induced precipitation in tungsten irradiated with a mixed energy neutron spectrum". In: *Acta Materialia* 165 (2019), pp. 51-61.

[28] A. Xu et al. "Ion-Irradiation Induced Clustering in W-Re-Ta, W-Re and W-Ta Alloys: An Atom Probe Tomography and Nanoindentation Study". In: *Acta Materialia* 124 (2017), pp. 71-78.

[29] S. P. Fitzgerald and D. Nguyen-Manh. "Peierls Potential for Crowdions in the bcc Transition Metals". In: *Physical Review Letters* 101.11 (2008), p. 115504.

[30] D. Nguyen-Manh, A. P. Horsfield, and S. L. Dudarev. "Self-Interstitial Atom Defects in bcc Transition Metals: Group-Specific Trends". In: *Physical Review B* 73.2 (2006), p. 020101.

[31] C. S. Becquart et al. "Microstructural Evolution of Irradiated Tungsten: Ab Initio Parameterisation of an OKMC Model". In: *Journal of Nuclear Materials* 403.1-3 (2010), pp. 75-88.
[32] T. Suzudo, M. Yamaguchi, and A. Hasegawa. "Migration of Rhenium and Osmium Interstitials in Tungsten". In: *Journal of Nuclear Materials* 467 (2015), pp. 418-423.

[33] T. Suzudo and A. Hasegawa. "Suppression of Radiation-Induced Point Defects by Rhenium and Osmium Interstitials in Tungsten". In: *Scientific reports* 6 (2016), p. 36738.

[34] N. Castin et al. "The Effect of Rhenium on the Diffusion of Small Interstitial Clusters in Tungsten". In: *Computational Materials Science* 177 (2020), p. 109580.

[35] C.-H. Huang et al. "Mechanism of Nucleation and Incipient Growth of Re Clusters in Irradiated W-Re Alloys from Kinetic Monte Carlo Simulations". In: *Physical Review B* 96.9 (2017), p. 094108.

[36] M. J. Lloyd et al. "Radiation-Induced Segregation in W-Re: from Kinetic Monte Carlo Simulations to Atom Probe Tomography Experiments". In: *The European Physical Journal B* 92.10 (2019), p. 241.

[37] J. S. Wróbel et al. "A First-Principles Model for Anomalous Segregation in Dilute Ternary Tungsten-Rhenium-Vacancy Alloys". In: *Journal of Physics: Condensed Matter* 29.14 (2017), p. 145403.

[38] D. Nguyen-Manh et al. "First-Principles Model for Voids Decorated by Transmutation Solutes: Short-Range Order Effects and Application to Neutron Irradiated Tungsten". In: *Phys. Rev. Mater.* 5 (2021), p. 065401.

[39] G. Kresse and J. Hafner. "Ab initio molecular-dynamics simulation of the liquid-metal-amorphous-semiconductor transition in germanium". In: *Phys. Rev. B* 49.20 (1994), pp. 14251-14269. DOI: 10.1103/PhysRevB.49.14251.

[40] G. Kresse and J. Furthmüller. "Efficient Iterative Schemes for Ab Initio Total-Energy Calculations using a Plane-Wave Basis Set". In: *Phys. Rev. B* 54.16 (1996), pp. 11169-11186. DOI: 10.1103/PhysRevB.54.11169.

[41] G. Kresse and J. Furthmüller. "Efficiency of Ab-Initio Total Energy Calculations for Metals and Semiconductors using a Plane-Wave Basis Set". In: *Comput. Mater. Sci.* 6.1 (1996), pp. 15-50. DOI: 10.1016/0927-0256(96)00008-0.

[42] P. E. Blöchl. "Projector augmented-wave method". In: *Phys. Rev. B* 50.24 (1994), pp. 17953-17979. DOI: 10.1103/PhysRevB.50.17953.

[43] John P. Perdew, Kieron Burke, and Matthias Ernzerhof. "Generalized gradient approximation made simple". In: *Phys. Rev. Lett.* 77.18 (1996), pp. 3865-3868. DOI: 10.1103/PhysRevLett.77.3865.

[44] James D. Pack and Hendrik J. Monkhorst. "Special points for Brillouin-zone integrations". In: *Phys. Rev. B* 13.12 (1976), pp. 5188-5192. DOI: 10.1103/PhysRevB.13.5188.

[45] D. Nguyen-Manh, M. Yu. Lavrentiev, and S. L. Dudarev. "Magnetic origin of nano-clustering and point defect interaction in Fe-Cr alloys: an ab-initio study". In: *J. Computer-Aided Mater. Des.* 14 (2007), pp. 159-169.
REFERENCES

29

[46] M. Muzyk et al. "Phase stability, point defects, and elastic properties of W-V and W-Ta alloys". In: *Phys. Rev. B* 84.10 (2011), p. 104115. ISSN: 10980121. DOI: 10.1103/PhysRevB.84.104115.

[47] A. van de Walle, M. Asta, and G. Ceder. "The alloy theoretic automated toolkit: A user guide". In: *CALPHAD* 26.4 (2002), pp. 539-553. DOI: 10.1016/S0364-5916(02)80006-2.

[48] D. Nguyen-Manh and S. L. Dudarev. "Trapping of He clusters by inert-gas impurities in tungsten: First-principles predictions and experimental validation". In: *Nucl. Instrum. Methods Phys. Res. B* 352 (2015), pp. 86-91.

[49] Daniel R Mason, Duc Nguyen-Manh, and Charlotte S Becquart. "An empirical potential for simulating vacancy clusters in tungsten". In: *Journal of Physics: Condensed Matter* 29.50 (2017), p. 505501.

[50] D. R. Mason et al. "Relaxation of microscopic and mesoscopic irradiation-induced defects in tungsten". In: *J. Appl. Phys.* 126 (2019), p. 075112. DOI: 10.1063.1.5098452.

[51] L. Messina et al. *Transport Properties of Transmutation Products in Concentrated W-Re-Os Alloys by Smart Fitting of Ab Initio Datasets*. Presented at the 19th International Conference on Fusion Reactor Materials (ICFRM19), New Modeling Methods Symposium, San Diego, CA, USA, 2019.

[52] *scipy.optimize.curve_fit*. [online] https://docs.scipy.org/doc/scipy/reference/generated/scipy Mixedudes.html. 2020 (accessed 23 September, 2020).

[53] *scipy.spatial.ConvexHull*. [online] https://docs.scipy.org/doc/scipy/reference/generated/scipy-temporal-convex-hull.html. 2020 (accessed 23 September, 2020).

[54] Otto Redlich and AT Kister. "Thermodynamics of Nonelectrolyte Solutions-xyt Relations in a Binary System". In: *Industrial & Engineering Chemistry* 40.2 (1948), pp. 341-345.

[55] A. Stukowski. "Visualization and Analysis of Atomistic Simulation Data with OVITO-the Open Visualization Tool". In: *Modelling and Simulation in Materials Science and Engineering* 18.1 (Dec. 2009), p. 015012.