Modeling refractory high-entropy alloys with efficient machine-learned interatomic potentials: defects and segregation

We first use the Gaussian approximation potential (GAP) framework [16, 17] to train an interatomic potential for the complete Mo–Nb–Ta–V–W system. With five elements and the aim to get a fairly general potential that is at least reasonably transferable to any alloy composition, constructing the training database is a challenge. We also require that the potential can describe not only bulk crystals, but also any form of defects that may be present or form due to irradiation. We found that with a reasonably sized training data set, using the many-body SOAP descriptor [18] that is typically used in GAPs leads to overfitting issues (showing significantly larger test errors than training errors). Converging a SOAP-GAP towards DFT accuracy for independent test data would likely require a prohibitively large training set due to the vast descriptor space with five elements. Our solution is therefore to rely on low-dimensional two- and three-body descriptors, which require less data but as we demonstrate can still achieve remarkably good accuracy for multicomponent alloys. In fact, it outperforms a SOAP-GAP trained to the same training data. Additionally, the use of only low-dimensional descriptors allows for a tabulation scheme that provides a significant computational speedup, as discussed in Sec. II.2.

We use the following expression for the total energy of a system of $N$ atoms

$\displaystyle\begin{split}E_{\mathrm{tot.}}^{\mathrm{GAP}}&=\sum_{i<j}^{N}E_{ij}^{\mathrm{rep.}}(r_{ij})+\delta_{\mathrm{2b}}^{2}\sum_{i<j}^{N}\sum_{p}^{M_{\mathrm{2b}}}\alpha_{ij,p}K_{\mathrm{2b}}(r_{ij},r_{p})\\ &+\delta_{\mathrm{3b}}^{2}\sum_{i,j<k}^{N}\sum_{t}^{M_{\mathrm{3b}}}\alpha_{ijk,t}K_{\mathrm{3b}}(\bm{q}_{ijk},\bm{q}_{t}).\end{split}$

(1)

Here, the first term is a purely repulsive screened Coulomb potential fitted to each element pair using all-electron DFT data [19]. Our methods for fitting and including the repulsive pair potentials are described in detail in Ref. [20]. The second term is a machine-learning pair potential, carried out as Gaussian process regression with the interatomic distance of each pair $ij$ as the descriptor. The final term is the three-body machine-learning term as a sum over all atom triplets $ijk$. $\alpha$ are the optimised regression coefficients, and $K_{\mathrm{2b}}$, $K_{\mathrm{3b}}$ are the kernel functions, for which we use the squared-exponential kernel with 1 Å as the length-scale hyperparameter. $M_{\mathrm{2b}}$ and $M_{\mathrm{3b}}$ are the number of representative pairs $p$ and triplets $t$ from the training structures used to carry out the regression. We use $M_{\mathrm{2b}}=20$ and $M_{\mathrm{3b}}=300$ for all element combinations. Values higher than $M_{\mathrm{3b}}=300$ provided similar accuracy, but showed an increasing (although small) tendency of overfitting. The machine-learning energy predictions are scaled by the parameters $\delta_{\mathrm{2b}}^{2}=10$ eV and $\delta_{\mathrm{3b}}^{2}=1$ eV. The GAP framework and its parameters are described in more detail in e.g. Refs. [17, 21].

The three-body descriptor calculated for each triplet of atoms $ijk$ is the vector [17]

$$\bm{q}_{ijk}=\begin{pmatrix}r_{ij}+r_{ik}\\ (r_{ij}-r_{ik})^{2}\\ r_{jk}\\ \end{pmatrix}f_{\mathrm{cut}}(r_{ij})f_{\mathrm{cut}}(r_{ik}),$$

(2)

which is invariant to permutations of $jk$. Smooth cutoff functions $f_{\mathrm{cut}}(r)$ are imposed on the bonds to the central atom $i$. We use 5 Å cutoff radii for both the two-body and three-body descriptors and all element pairs and triplets. Initially, we experimented with longer cutoffs for the pairwise potentials (to capture possible long-range interactions) and shorter three-body cutoffs, but found that 5 Å cutoffs for both descriptors lead to better accuracy.

The potential is trained to total energies, forces, and virial stresses obtained from DFT calculations using vasp [22, 23, 24, 25] for a variety of Mo–Nb–Ta–V–W structures. The training set includes bulk Mo–Nb–Ta–V–W bcc crystals sampled at all compositions, including the pure elements. The majority of the alloys are randomly ordered, but we also found it necessary to include ordered alloys. Defected structures (up to 5 vacancies and self-interstitial atoms) are only included in the pure elements and in equiatomic random MoNbTaVW alloys. Liquid structures are included for pure elements and all equiatomic binary, ternary, quaternary, and quinary alloys of Mo–Nb–Ta–V–W, but with most structures for the equiatomic MoNbTaVW high-entropy alloy (HEA). The training set additionally includes HEA crystals with one randomly added interstitial atom that is relatively close (but not too close) to its neighbour atom. These structures ensure that the repulsion between all elements in crystals is captured. Surface structures are included for the pure elements along with a few HEA surfaces to ensure some transferability to alloy surfaces, although we do not specifically target surface properties. Our strategy for constructing the training structures and the training process is described in more detail in the Supplemental material online. The potential is trained using the quip code (https://github.com/libAtoms/QUIP).

Carrying out the Gaussian process regression terms of Eq. 1 is computationally expensive due to the loop over the (sparsified) training points $M_{\mathrm{2b}}$, $M_{\mathrm{3b}}$. Because our GAP only contains pure two- and three-body-dependent terms, it is possible to create a tabulated version of the potential, which we will call tabGAP. Creating computationally efficient tabulated machine-learning potentials was initially demonstrated by Glielmo et al. [26] and recently further developed in the flare framework [27, 28], although details vary from our approach. The pairwise terms can be evaluated as 1D interpolations between tabulated pair energies and the three-body terms as 3D interpolations between tabulated triplet energies. We use a grid of $(r_{ij},r_{ik},\cos\theta_{ijk})$ points for the three-body interpolation. With $S^{\mathrm{1D}}$ and $S^{\mathrm{3D}}$ representing 1D and 3D cubic spline interpolations, the total energy of the tabGAP becomes

$$E_{\mathrm{tot.}}^{\mathrm{tabGAP}}=\sum_{i<j}^{N}S_{ij}^{\mathrm{1D}}(r_{ij})+\sum_{i,j<k}^{N}S_{ijk}^{\mathrm{3D}}(r_{ij},r_{ik},\cos\theta_{ijk}),$$

(3)

where the repulsive pair potential and the GAP pair potential in Eq. 1 are merged into one 1D spline per element pair. With five elements, there are 15 unique pairs and 75 triplets (with symmetry on the $jk$ elements), each requiring a tabulated grid of energies. With sufficiently dense grids, the tabGAP is virtually identical to the original GAP. The convergence of the interpolation error with increasing numbers of grid points is demonstrated in the Supplemental material online.

We have implemented this cubic-spline-interpolated potential as a pair_style for lammps [29] (based on the implementation from the open-source flare code [27, 28]). The speed-up compared to the original GAP evaluated as in Eq. 1 is around two orders of magnitude with the current implementation, making it comparable in speed to traditional angular-dependent potentials like Tersoff and MEAM potentials [30, 31, 32]. The potential file is available from [33] and the tabulation and lammps implementation from https://gitlab.com/jezper/tabgap.

All DFT calculations are carried out with vasp [22, 23, 24, 25] with PAW potentials [34, 35] (_pv for Ta and _sv for all other elements), the PBE GGA exchange-correlation functional [36], a 500 eV cutoff energy, 0.15 Å^-1 maximum $k$-point spacing on $\Gamma$-centred Monkhorst-Pack grids [37], and 0.1 eV Methfessel-Paxton smearing [38]. All molecular statics and dynamics simulations are done with lammps [29] with a custom implementation of the tabulated machine-learned potential (tabGAP) as discussed above.

The average lattice constant and mixing energy of bulk equiatomic HEA properties are obtained by relaxing 50 different randomly ordered 2000-atom bcc systems. The simulation cells were cubic during the energy and pressure minimization. Allowing for noncubic relaxation produced cubic systems within the statistical uncertainty of the cubic bcc lattice constant. The elastic constants are computed for the same relaxed 50 HEA systems. The bulk modulus is obtained from volume-energy fits to the Birch-Murnaghan equation of state [39]. The remaining elastic constants are solved from parabolic fits to strain-energy curves of the tetragonal and trigonal deformation modes (see e.g. Ref. [40]), using $\pm 2\%$ strain intervals and optimising the atom positions at every strain.

Nudged elastic band calculations for vacancy migration barriers are performed in boxes of 128 atoms using lammps. The migration energy is obtained from the saddle point of the converged barrier. For relaxing single self-interstitials we use boxes of 251 atoms. For the formation energies of single vacancies we use 54-atom boxes in both vasp and lammps. We optimise both the positions and relax the box size to zero pressure for the vacancy and self-interstitial calculations. The formation energy for an $A$ vacancy (where $A$ is Mo, Nb, Ta, V, or W) is calculated as

$$E_{\mathrm{f}}^{A\mathrm{vac}}=E_{\mathrm{vac}}-E_{\mathrm{bulk}}+\mu_{A},$$

(4)

where $E_{\mathrm{vac}}$ is the total energy of the relaxed vacancy system and $E_{\mathrm{bulk}}$ is the total energy of the relaxed HEA bulk with the vacancy filled by an atom of element $A$. The chemical potential $\mu_{A}$ is approximated simply as the energy per atom of $A$ in its ground state (bcc bulk), which is close to the real chemical potential estimated using other methods [41]. The vacancy relaxation volume is calculated as

$$\Omega_{\mathrm{rel.}}^{A\mathrm{vac}}=\frac{V_{\mathrm{vac}}-V_{\mathrm{bulk}}}{\Omega_{A}},$$

(5)

where $V_{\mathrm{vac}}$ and $V_{\mathrm{bulk}}$ are the total volumes of the relaxed vacancy and bulk HEA systems and $\Omega_{A}$ is the volume per atom of pure bcc $A$.

Relaxing self-interstitial atoms in the HEA often results in stable interstitial configurations far from the initial configuration, as many local atomic environments provide no local minima. This, combined with the fact that many interstitial configurations are mixed dumbbells, makes the choice or construction of the corresponding reference bulk system somewhat ambiguous. Hence, we define and compute the formation energy of an $A$–$B$ interstitial dumbbell configuration as the formation energy of the entire interstitial system with $N=251$ atoms ($E_{\mathrm{f}}^{\mathrm{SIA}}$), and then compare it to an average defect-free HEA bulk reference to get only the formation energy associated with the interstitial as

$$E_{\mathrm{f}}^{A\text{-}B}=E_{\mathrm{f}}^{\mathrm{SIA}}-NE_{\mathrm{mix}}^{\mathrm{HEA}}.$$

(6)

The formation energy and mixing energy (per atom) of a HEA system of $N$ atoms with total energy $E$ is $E_{\mathrm{f}}=E-\sum_{i}^{N}E_{i}$ and $E_{\mathrm{mix}}=E_{\mathrm{f}}/N$, where $E_{i}$ is the energy per atom of the pure ground state of atom $i$. $E_{\mathrm{mix}}^{\mathrm{HEA}}$ is computed as the average mixing energy of randomly ordered HEA systems and is $-41.85$ meV/atom in the tabGAP (as listed in Tab. 2 and discussed later).

Hybrid Monte Carlo–Molecular dynamics (MC+MD) simulations are carried out as implemented in lammps. MD is done in the $NPT$ ensemble at different constant temperatures and zero pressure. For the single-crystal simulations we use 100 trial MC atom swaps every 10 MD steps using boxes of 6750 atoms. The simulations are continued until the potential energy and short-range order parameters have sufficiently converged, although for at least 1 million MC steps and about 8 million MC steps for the lower temperatures. The short-range order parameter for element pairs $AB$ with interatomic bonds in the interval $\Delta r_{ij}$ is computed using the definition

$$S_{\Delta r_{ij}}^{AB}=1-\frac{p_{\Delta r_{ij}}^{AB}}{c_{B}},$$

(7)

where $p_{\Delta r_{ij}}^{AB}$ is the probability of finding a $B$ atom around $A$ in the neighbour shell $\Delta r_{ij}$ and $c_{B}$ is the concentration of $B$ in the alloy. For first-nearest-neighbour pairs, we used the interval $\Delta r_{ij}=[0,3]$ Å and for the second-nearest neighbours $\Delta r_{ij}=[3,3.9]$ Å.

For the void and dislocation loop MC+MD simulations, we allow for more MD relaxation with 10 MC trial swaps every 10 MD steps in boxes of around 20 000 atoms. Here, we run the simulations until the concentrations around the defects have stabilised (disregarding the surrounding bulk), which only required a few hundred thousand MC steps.

The defect annealing simulations are done using boxes of 250 000 atoms with 10 000 randomly created Frenkel pairs, corresponding to a supersaturated (4%) defect concentration. The atomic positions are first optimised, which already annihilates many Frenkel pairs and causes some initial clustering of interstitials and vacancies. This is followed by a 1 ns $NPT$ annealing run at 2000 K and zero pressure. The final frames are then again optimised to allow a more reliable analysis of the final defect structure. The dislocations are identified using the dislocation extraction algorithm (DXA) [42] in ovito [43]. Interstitials and vacancies are found using the Wigner-Seitz method and grouped into clusters with the cutoff radius between the second and third nearest neighbour for vacancies and between the third and fourth for interstitials.

Before using the tabGAP to study segregation and the defect structure of the equiatomic MoNbTaVW high-entropy alloy, we here briefly demonstrate the accuracy of the potential. Table 1 lists root-mean-square errors (RMSE) of the tabGAP compared to DFT for a range of structures used as test sets. All atoms in the bulk crystals have been randomly displaced from the perfect lattice positions to produce significant interatomic forces. The accuracy of the tabGAP is around 3 meV/atom and 0.05–0.2 eV/Å for any given composition of Mo–Nb–Ta–V–W as well as for HEA lattices containing vacancies or self-interstitial atoms. Some further test results are provided in the Supplemental material, where we show that the tabGAP also reproduces the W–Mo alloy training dataset from Ref. [44] with similar accuracy as in Tab. 1, verifying that the potential is accurate also for the binary alloys.

To further verify that the tabGAP can describe the phase stability and elastic response of arbitrary alloy compositions, we compare the 0 K equation of state for a variety of alloys between DFT and the tabGAP in Fig. 1. From the energy-volume data, we obtain the bulk moduli and equilibrium mixing energies from fits to the Birch-Murnaghan equation of state [39]. Fig 1(a) shows the energy-volume data for randomly ordered alloys at all equiatomic alloy compositions (one system per composition) and Fig. 1(c) shows data for 20 different quinary alloys randomly ordered and at randomly sampled compositions. The corresponding bulk moduli and mixing energies are shown in Figs. 1(b) and (d), compared between DFT and the tabGAP. The tabGAP reproduces the DFT data within only a few meV/atom for mixing energies and a few GPa for the bulk moduli. Also noteworthy is that the agreement between DFT and the tabGAP is good for the entire wide range of volumes for all alloys in Figs. 1(a) and (c), where each curve spans an energy difference around 1–2 eV/atom.

In table 2, we list average properties of MoNbTaVW predicted by the tabGAP and computed from 50 relaxed 2000-atom systems. The lattice constant in the tabGAP (3.195 Å) is close to the experimental value 3.1832 Å [45].

Reproducing liquid properties for arbitrary compositions, densities, and temperatures with very high accuracy cannot be expected with a simple three-body-dependent potential. As seen in Tab. 1, the RMSEs for liquids are up to 10 times higher than for any crystalline structures. Nevertheless, we confirmed that the tabGAP still provides a reasonable descripton of the liquid phase and melting. The melting point of the HEA predicted by the tabGAP is $2760\pm 20$ K as determined by $NPT$ molecular dynamics (MD) simulations of a solid-liquid interface. A rough estimate of the experimental melting point can be taken as the average of the pure-element melting temperatures, yielding 2961 K [46]. However, previously we found that even highly accurate many-body GAPs tend to underestimate the melting point slightly [47]. The corresponding average of the pure-element many-body GAPs from Ref. [47] is 2796 K [47], which is very close to the HEA melting point predicted by the tabGAP.

We first apply the tabGAP in a study of ordering and segregation in the single-crystal equiatomic HEA using hybrid MC+MD simulations. Previous studies using cluster expansion models and DFT calculations have established that there is a strong preference for local ordering of Mo–Ta binary crystals due to their relatively low mixing energy [48, 11, 49]. Additionally, it has been previously shown that including lattice relaxation is crucial in order to reproduce the correct phase stabilities at finite temperatures [13]. Relaxation effects and realistic thermal vibrations are here naturally included in the $NPT$ MC+MD simulations with the tabGAP. Since the tabGAP is trained to both random and ordered alloys, including all the binary alloys and MC+MD-optimised ordered systems obtained in an active-learning fashion (see the Supplemental material), we are confident that it can accurately model order and disorder in the HEA.

Fig. 2 shows the results from the MC+MD simulations. The 1NN and 2NN short-range order parameters as functions of temperature in Fig. 2(a) show, consistent with previous studies, that Mo and Ta are locally ordered at low temperatures. A strongly negative 1NN Mo–Ta short-range order (SRO) value and a positive 2NN SRO value indicate the presence of MoTa in the CsCl (B2) order, which is confirmed by visually analysing the lattice. Fig. 2(b) shows snapshots of the systems at two temperatures. In addition to Mo–Ta ordering, Mo–Nb and W–V 1NN pairs are also favoured. W–V pairs stand out by also having negative 2NN SRO value, indicating a more complex structure than the CsCl order, as was also observed in Ref. [11]. All SRO values start approaching zero after 300 K, initiating the transition from the ordered phases to a random solid solution. Still, a relatively strong local ordering of mainly Mo–Ta is retained up to temperatures beyond 1000 K.

The preferential ordering and segregation can largely be understood by the mixing energies of the binary alloys. Fig. 2(c) shows mixing curves at 0 K for all binaries as given by the tabGAP (which are identical within 2.9 meV/atom to DFT data, see figure S2 in the Supplemental material). The connected data points in Fig. 2(c) are mixing energies for randomly ordered alloys (using 1024-atom boxes) and the large solid points are the ordered CsCl phase. The mixing curves show that Mo–Ta alloys have the lowest mixing energies and that the ordered MoTa phase is by far the most preferred binary alloy. The mixing energies also suggest that Nb–V alloys are the least favoured, with clearly positive mixing energies for both random and the ordered CsCl phase. This is reflected in the MC+MD simulations by the fact that Nb and V are strongly segregated from each other, showing positive SRO values for both 1NN and 2NN pairs.

The presence of defects introduce local lattice relaxation and long-range strain fields that may influence the preferential ordering and segregation of the different-sized atoms in the HEA. We use the tabGAP to investigate the energetics of single vacancies and self-interstitials in the HEA. The chemical complexity provides exceedingly many possible configurations for even these simple defects, which calls for a statistical treatment. Here, we only consider randomly ordered HEAs with randomly added single vacancies and interstitials, followed by relaxation. Fig. 3(a)–(b) shows distributions of formation energies and relaxation volumes of single vacancies in the HEA. The vacancies are created in 10 different HEA systems. For each vacancy system, reference bulk systems are created by filling the vacancy with each element separately, so that we in total obtain data for 50 different vacancies. Fig. 3 also shows results from DFT calculations to further validate the predictive accuracy of the tabGAP. The tabGAP slightly underestimates the average formation energy and the relaxation volume compared to DFT, but overall provides similar distributions. The average vacancy formation energy is 3.3 eV according to DFT and 3.1 eV according to the tabGAP. The average relaxation volume is $-0.35$ at. vol. in DFT and $-0.40$ at. vol. in the tabGAP.

The single-vacancy data are separated by element in Fig. 3(c), revealing additional insight. First, there is no correlation in the formation energy of vacancies between the elements, with energies in the range 2.5–4 eV for all elements. There is, however, a clear separation between the elements for relaxation volumes. Filling a vacancy with a V atom, the smallest atom, cause the least relaxation around the vacancy with most relaxation volumes around $-0.2$ to $-0.3$ at. vol. Vice versa, being the largest atoms, Nb and Ta vacancies produce significantly stronger relaxation with relaxation volumes around $-0.4$ to $-0.5$ at. vol.

Fig. 4 shows distributions of migration energies for single vacancies in the HEA. In total, we calculated over 1100 migration barriers in different randomly ordered HEAs with the tabGAP to obtain reliable statistics. Fig. 4 reveals that the migration energies for each element are roughly normally distributed and span a wide energy range of around 1 eV. There is again a clear separation between the elements of the migrating atom. The group 6 elements W and Mo have the highest migration energies. The average migration energy for W is 1.56 eV and for Mo 1.48 eV. The group 5 elements show significantly lower migration energies, with the averages 1.06 eV for Nb, 1.12 eV for Ta, and 1.11 eV for V. The trend is consistent with the vacancy migration energies in the pure elements, where W and Mo have by far the highest migration energies (around 1.7 and 1.3 eV [50]) while the group 5 elements all have migration energies around 0.6–0.7 eV [50]. The result of alloying on the vacancy migration is thus a wide distribution of migration energies where the averages are reduced for W and Mo and increased for Nb, Ta, and V.

To explore the preferential chemical and geometrical configurations of self-interstitial atoms, we relaxed 1000 HEA systems each containing one randomly inserted interstitial atom. Relaxation produces mainly dumbbell configurations of atom pairs. Distributions of the relaxed dumbbell configurations and formation energies are shown in Fig. 5. Fig. 5(a)-(b) reveal that there is a clear preference for V-containing interstitial dumbbells, and in particular pure V–V dumbbells which make up the majority of stable interstitials. V–V dumbbells also have on average the lowest formation energy, as shown in Fig. 5(b). Note that there are not enough non-V dumbbells to provide reliable distributions and average formation energies, but the trend of V-containing dumbbells having the lowest formation energies is clear. This can again be understood by the fact that V is the smallest atom and therefore prefers shorter interatomic bonds than the other elements.

Fig. 5(c) shows the relaxed dumbbell directions. Based on the pure elements, one would expect mainly $\hkl<111>$ directions, which are the lowest-energy configurations in Nb, V, and Ta, and directions between $\hkl<111>$ and $\hkl<221>$, which are preferred in pure Mo and W (i.e. the $\hkl<11\xi>$ interstitial revealed in Ref. [51]). Although these directions make up a large fraction of the relaxed interstitial dumbbells, Fig. 5(b) shows that \hkl¡110¿ dumbbells are observed in almost equal numbers to \hkl¡221¿ dumbbells. By looking at the relaxed systems, we find that the preferred direction is strongly influenced by the chemical environment. Stable configurations are found by a competition between low-energy dumbbell directions and the possibility to form V–V or other V-containing pairs. We found no statistically significant correlation between element pair and dumbbell direction, and hence only show the total distribution of directions in Fig. 5(c).

We also find that the preferential formation of V-containing dumbbell interstitials guides the mobility and migration mechanisms of single interstitials. In contrast to the pure metals, where single interstitials migrate one-dimensionally along the \hkl¡111¿ direction in all elements [51], interstitial migration in the HEA is three-dimensional. From constant-temperature MD simulations, we observe that the interstitials migrate mainly through connections of V atoms, so that the V-containing dumbbells is most likely to migrate to other neighbouring V atoms. If no V atoms are present in the local environment, the V interstitial remains stationary for significantly longer times.

Many of the above-discussed observations are consistent with the results by Zhao in a recent DFT study of point defects in a similar bcc alloy (VTaCrW) [41]. Namely, (1) there is no clear difference in the vacancy formation energies between the different elements, (2) the relaxation volume of vacancies directly correlates with the size of the removed atom, (3) vacancy migration energies vary significantly between the elements, (4) self-interstitial atoms are most stable when containing smaller atoms, like V–V dumbbells, and (5) self-interstitial dumbbells adopt different directions depending on the chemical environment, with $\hkl<110>$ as one of the most probable low-index directions, in clear contrast to the pure elements.

Going beyond single vacancies and interstitials, we now investigate how clusters of defects influence their local chemical ordering. The most stable vacancy- and interstitial-type defect clusters formed during irradiation in bcc materials are small voids and dislocation loops. We carry out MC+MD simulations for HEAs containing (separately) voids of two sizes (15 vacancies and 65 vacancies) and interstitial dislocation loops with the Burgers vectors 1/2\hkl¡111¿ and \hkl¡100¿ and diameter 2 nm.

The simulations show that elemental segregation around the defect clusters is clearly energetically favoured. Fig. 6 shows the concentration profiles as functions of the distance from the void centre for both void sizes. Snapshots of the inner surface atoms of the voids are shown as insets. The equilibrated inner void surface is almost completely covered by Nb. This is understandable as Nb is (along with Ta) the largest atom and also has the lowest surface energy (lower than Ta) [47]. Fig. 6 also shows that there is an excess of V at the interface of the Nb-covered inner surface and the bulk. The snapshots reveal that this is because V atoms tend to decorate the edges of the void surface. The presence of the large Nb atoms and the polyhedral shape of the voids make these sub-surface regions compressed, which explains why they are preferentially occupied by V, the smallest atom.

Interstitial dislocation loops produce strong strain fields with a locally compressed lattice, which our previous observations suggest should be favoured by V. Indeed, Fig. 7 shows that essentially 100% of the atoms within the strain field of the dislocation loops are V after the MC+MD simulations reach a steady state. This is also in line with the previous observation that the single V–V interstitial dumbbell is the most favoured interstitial configuration. The results for both voids and interstitial loops are almost identical in MC+MD simulations at both 300 K and 1000 K, showing that the preferential segregation, driven mainly by the atom size, is strong and occurs also at high temperatures.

In the previous sections, we have described the energetics of point defects as well as segregation and ordering in the bulk and around possible radiation-induced defect clusters in the HEA. It is not straightforward to combine all these results into an understanding of how the HEA behaves during irradiation and the subsequent defect annealing and recovery. To this end, we carry out annealing MD simulations of a supersaturated concentration of defects and observe how the defect structure evolves over time. We prepare three 250,000-atom HEA samples with randomly inserted interstitial atoms and vacancies corresponding to a defect concentration of 4%. After an initial minimisation of positions, the samples are annealed at zero pressure and 2000 K for 1 ns. For comparison, we also anneal a pure W sample in the same way. The results are summarised in Fig. 8. Animations of the defect evolution are provided in the Supplemental material online.

Fig. 8 shows striking differences between W and the HEA during annealing. In W, interstitial atoms are extremely mobile and join to form 1/2\hkl¡111¿ dislocation loops already during the first few picoseconds. This is clear from Fig. 8(a), which shows that the fraction of interstitials in clusters reaches 100% very rapidly. Over time, the initially small dislocation loops grow by migrating and coalescing with other nearby dislocation loops, eventually resulting in only a few large loops as seen in Fig. 8(b). The coalescence of loops is in Fig. 8(a) reflected by the rapid decrease in the dislocation density. Compared to interstitials, vacancies in W migrate slowly (migration energy 1.7 eV [50]). Additionally, the binding energies of small vacancy clusters in W is close to zero or even repulsive [20], which further limits the formation of vacancy clusters. Only around 20% of the vacancies in W are in clusters of two or more, as shown in Fig. 8(a).

The defect evolution in the HEA is in many ways different. Fig. 8(a) shows that the recombination of defects is more efficient, especially during the early stage of the annealing simulations. Fig. 8(b) shows that the final defect structure in the HEA only contains small dislocation loops. These dislocation loops form already in the early part of the simulations but, unlike in W, then remain stationary and similar in size throughout the rest of the simulation. This is evidenced by the constant dislocation density over time in Fig. 8(a).

Contrary to W, vacancies in the HEA are more mobile (Fig. 4) and about 40% of them are in clusters of two or more. As observed in Sec. III.3, single interstitials are most stable as V–V dumbbells and migrate three-dimensionally, preferentially through connections of V–V bonds. This is fundamentally different from W, where interstitials mainly migrate one-dimensionally along close-packed \hkl¡111¿ directions with a very low migration energy. The consequence of this difference is the much more efficient defect recombination in the HEA seen in Fig. 8(a), as both interstitials and vacancies are mobile and can explore their surrounding through 3D migration. Furthermore, because the dislocation loops are unable to move and interstitials are more likely to recombine with nearby vacancies, the overall defect structure does not significantly change over time. In particular, the dislocation loops remain very small. All loops are around 1–1.5 nm in diameter with the Burgers vector 1/2\hkl¡111¿, with most of the loops around 1 nm and consisting of 20–30 interstitials. Smaller interstitial clusters are also mostly parallel \hkl¡111¿ interstitials, but cannot be classified as dislocation loops and are thus not identified by the dislocation extraction algorithm. Additionally and in line with our observations in Secs. III.3 and III.4, all interstitial clusters are enriched by V.