Molecular dynamics simulations of high-dose damage production and defect evolution in tungsten

Tungsten has been chosen as the plasma-facing wall material in fusion reactors, due to its high density and melting point. The wall material will not only be sputtered at the surface, but also damaged deep inside the material by energetic particles. We investigate the high-dose damage production and accumulation by computational means using molecular dynamics. We observe that the choice of interatomic potential drastically affects the evolution. The structure and stability of the obtained defect configurations are validated using a quantum-accurate Gaussian approximation potential. © 2021 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
Tungsten has due to its high density and melting temperature been chosen as one of the materials to be used in parts of fusion reactors that are exposed to the highest heat and particle fluxes. These parts are the first wall material and the divertor of the reactor. The wall materials will in addition to being sputtered by the energetic ions and neutrals impacting the surface also be damaged deeper inside the material, far away from the surface. As the bulk damage and defect production ultimately will affect the mechanical properties of the material, both experimental and computational studies have been carried out to investigate the defect production in tungsten [1,2] .
A recent study by Reza et al. [3] investigated the defect concentration as a function of dose in W at room temperature at both low and high doses. They found, using Transient Grating Spectroscopy (TGS), that Frenkel pair saturation occurs at doses between 0.06 and 0.1 dpa, with defect fractions between 0.004 and 0.005, similar to what has been found previously in other metals [4][5][6] . Another study investigated the vacancy clustering by Positron Annihilation Spectroscopy (PAS), and showed that vacancies at low doses and room temperature will be found mainly in small clusters or isolated [7] . Transmission Electron Microscopy (TEM) studies have also been carried out on irradiated tungsten, and investigated the dislocation density evolution as a function of dose and tempera-ture [8] . These studies showed a larger dislocation density at lower temperature, which increases as the dose increases in the sub-dpa regime. Several studies have investigated the swelling of irradiated tungsten and shown that swelling is a complex phenomena depending on many factors. It has for instance been shown that an increase in temperature will first increase swelling, but then if a high enough temperature is reached, the swelling decreases [9] . Another study showed that the swelling can decrease by increasing the dose [10] .
On the simulation side, different Molecular Dynamics (MD) studies have investigated the effect of the primary knock-on (PKA) energy [11][12][13][14] , temperature [11,13] , electronic stopping [15,16] , and interatomic potential on the defect production and morphology [14,17] . Studies have found that different energies will yield different amount of defects and differences in cluster size distributions [13,14] . It has also been observed that the potential choice will affect both of these [14] . Early studies of damage overlap in metals showed that heat spikes in dense metals can recombine pre-existing damage [18,19] , in qualitative agreement with the experimental observations [4,5] . Recent studies have focused the effects of pre-existing defects on cascade damage in tungsten, and analyzed the defect morphology change in much greater detail [20][21][22] . It has been shown that a defect close to the cascade region reduces the number of surviving new defects compared to a pristine sample. It was also seen that the defect morphology can change. For instance, the Burgers vector of an existing dislocation loops can change due to a nearby cascade.
As the experiments usually are carried out at high doses, where the same region is affected by several cascades, and with the knowledge that the effect of a cascade overlapping with preexisting damage can change the morphology of the defect, we here carry out massively overlapping cascades to bridge the gap between simulations and experiments. These simulations will both accumulate a higher dose as well as account for the overlapping effects. Similar simulations have been carried out in Fe and Fe-Cr-alloys [23][24][25] and in equiatomic multicomponent alloys [26][27][28][29][30] , where in the latter a good agreement with experiments was obtained by comparing experimental and simulated Rutherford Backscattering Spectroscopy in a channeling direction (RBS/c) spectra [31,32] .
It is widely known that in molecular dynamics simulations, the results depend on the choice of interatomic potential. For many materials there exist several different potentials, fitted to reproduce different properties accurately, however, usually at the expense of describing other properties well. In recent years, new interatomic potentials using machine-learning methods have become popular. These potentials are often flexible enough and systematically improvable towards the accuracy of density functional theory (DFT), enabling MD simulations with almost quantum accuracy. One of the most successful types of machine-learning potentials is the Gaussian approximation potential (GAP) [33] , which has recently been applied to model radiation damage in tungsten [34] .
In this article, we carry out 20 0 0 cumulative cascades reaching a dose of 0.18 dpa to study the defect production and evolution in tungsten. We use three different commonly used interatomic potentials and study the point-defect evolution and dislocation formation and evolution and discuss the differences between the potentials. To validate our results and to observe possible changes in defect morphology, we relax simulation cells at various doses using the GAP W potential.

Cascade simulations
All cascade simulations were carried out with the MD code PARCAS, implemented with an adaptive timestep to accurately follow the trajectories of atoms [35][36][37] . The simulation cell was 60 × 63 × 66 unit cells in the x -, y -and z-directions, resulting in a cell with the side length of about 20 nm and containing approximately 50 0 0 0 0 atoms. Before any cascades were initiated, all cells where relaxed and thermalized at 300 K. Electronic stopping was applied as a friction force [37] on all atoms with a kinetic energy over 10 eV [15] .
The cascade simulations were carried out in a two step manner, with their own parameters; PKA simulation and relaxation simulation. These two steps were carried out 20 0 0 times in the same simulation cell to accumulate damage and reach a dose of about 0.18 dpa, according to the NRT-equation [38][39][40] with a threshold displacement energy of 70 eV, without the arc-dpa correction [6,41] . 70 eV was chosen instead of the commonly used 90 eV for tungsten, as the experimental study used a similar value in their study [3] . The use of 90 eV would not affect any of the results nor the evolution, only the maximum dose in our simulations would be 0.14 dpa instead of 0.18 dpa. During the PKA simulation, an atom in the centre of the cell was given 10 keV energy in a random direction and followed for 20 ps. In experiments, a wide spectrum of recoil energies would be seen, ranging from eVs to hundreds of keVs. We chose 10 keV as it is below the subcascade splitting energy and still high enough to accumulate a high enough dose. Previous studies by Gilbert et al. have shown that under fusion-relevant conditions a large fraction of PKAs are still seen at tens of keV in W [42] and in Fe it was found that most of the dose was accumulated by PKAs in the keV regime [43] . Previous simulations with monoenergetic PKAs have also shown very good agreement with experiments [32] . A thermostat was applied to atoms close to the periodic boundaries to cool down the system, but far enough from the cascade region to not affect its evolution. A barostat was not applied during this simulation step. During the relaxation part, the thermostat was applied to all atoms to cool the cell down to 300 K, which avoids temperature increases during the series. This step lasted 10 ps, enough to equilibrate the system back to 300 K. A pressure control was also applied to remove any excess pressure that the damage build-up could introduce. After these two steps, the cell was randomly shifted in all directions to obtain a homogeneous irradiation of the cell. The thermostat and barostat were those by Berendsen et al. [44] .
Three separate simulations sets were carried out for each of the potentials to obtain statistics. All given numerical results are the average value of the three runs.

GAP relaxation
To validate the results from the EAM potentials, we run short high-temperature annealing simulations on selected systems using the W GAP from Ref. [34] . The GAP was trained to density functional theory data for a wide range of structures and was shown to accurately reproduce properties relevant for radiation damage, such as the temperature dependence of elastic properties and the energetics of vacancy and self-interstitial clusters and dislocation loops [34] . A set of damaged systems from the cascade simulations are selected at various doses and from all three EAM potentials. The systems are initially relaxed by rescaling the lattice to the GAP-predicted lattice constant and optimising the positions, followed by an annealing MD run at 20 0 0 K and zero pressure in the NP T ensemble using the Nosé-Hoover thermostat and barostat [51,52] . Finally, we performed an energy minimization at zero pressure to allow for reliable analysis. The computational cost of the GAP ( > 10 0 0 times heavier than EAM) limits the simulations to short annealing runs, but 40 ps proved to be enough for possible unstable small and medium-sized defect clusters to transform to reach a more realistic defect structure.

Analysis
To obtain defect statistics and evolution throughout the series, we analyzed the cells for point defects and extended defects. Point defects, i.e. interstitials and vacancies, were identified with the Wigner-Seitz method [36] . Extended defect were identified by cluster analysis of the point defects and dislocations were detected with the DXA method implemented in OVITO [53,54] . The cutoffs for the cluster analysis were the mid-point between the 2nd nearest neighbour (2NN) and 3NN distances for vacancies and between the 3NN and 4NN distances for interstitials [55,56] . For all clusters, the net balance of vacancies and interstitials were calculated in order to obtain the correct defect size. The volumetric swelling was obtained by comparing the equilibrated cells with the initial undamaged relaxed cell.

Results
The number of defects (interstitials + vacancies) as a function of dose can be seen in Fig. 1 for all three interatomic potentials. The results are the average of three separate runs for each interatomic potential. From the figure we can see that the choice of interatomic potential drastically affects the damage build-up in the material. Even though the single cascades and very low-dose results are quite similar, the difference in defect stabilities in different interatomic potentials will dictate the evolution. The AT-ZN potential shows a saturation of defect concentration, whereas the two other potentials show an almost linear increase still at the higher investigated doses. Looking at the defect concentration in the GAP relaxed simulation cells (coloured markers), we observe almost no change apart from a slight decrease in the defect concentration. The black stars (Exp) are the data from Ref. [3] , converted into units of defect concentration, which was obtained by 20 MeV self-ion irradiation at different doses. The data are obtained from Fig. 4 therein, and is based on TGS measurements, which are converted into Frenkel pair concentrations with some assumptions. The main assumptions are that all defects are point defects and that all defects are isolated, which is not true especially for interstitials, which form dislocation loops. The TGS measurement is sensitive to the dislocation line, not the whole area, which may lead to errors in the dislocation densities utilizing these assumptions.
The more detailed evolution and clustering behaviour of the point defects, the GAP relaxation results and swelling are discussed in more detail separately in the next subsections.

Vacancy evolution
Looking at the overall vacancy number evolution (half the value given in Fig. 1 ) and the detailed vacancy cluster evolution in Fig. 2 , for all three potentials, several trends can be seen. As noted above, the evolution for the different potentials are different, apart from at very low doses. At low doses, mainly single vacancies and small vacancy clusters are formed in the simulations cells, for all three potentials.
In the AT-ZN potential, we observe that after a dose of 0.05 dpa, not too much is happening to neither the overall number of vacancies nor to the vacancy clusters. Additionally, the clustered fraction of vacancies, about 0.4, is stable at doses higher than this, as seen in Fig. 3 .
For the two other potentials, DND-BN and MS, a different trend is seen at higher doses compared to the AT-ZN potential. We observe an increase in the number of vacancies in all cluster sizes, explaining the overall trend of increasing defect concentration as the dose increases. The clustered fraction of vacancies is seen to increase as the dose increases, whereas the DND-BN shows a much lower fraction than the MS (and also the AT-ZN) over the interval investigated.
Overall, a clustered fraction of 0.2-0.4 is observed at higher dose. This shows that most vacancies are isolated, and looking at Fig. 2 , second most are in di-vacancies or small clusters. Only a few vacancy clusters (under 10) with sizes larger than ten vacancies are seen in all potentials.
The number of vacancy clusters as a function of dose can be found in the Supplementary material online for all three potentials, Fig. S1.

Interstitial evolution
Looking at the interstitial evolution in our simulations, again we observe different trends for the different potentials, Fig. 4 . In contrast to vacancies, we observe mainly large interstitial clusters, explained by their much higher mobility and therefore more efficient clustering. At doses above about 0.05 dpa, almost all interstitials are in clusters of size larger than 50. Only in the very beginning do we observe single interstitials and small interstitial clusters. The difference between potentials is seen in the intermediate range, clusters between 10 and 50, where DND-BN and MS show a similar trend, which differs from the AT-ZN one. In AT-ZN, almost no clusters in this regime exists, only for a brief moment before the large (50+) clusters overtake the evolution. For the two other potentials, we observe a large amount of clusters in the size range of 10 to 50 interstitials (up to a dose of 0.05 dpa). Looking more closely at the number of clusters in these larger size ranges, shown in Supplementary material Fig. S2, we can see that there are very few in the AT-ZN potential (2-3 size 50+ clusters) and a bit more in the other potentials (about 10 size 50+ clusters). The clustered fraction for all potentials is above 0.9 at all doses except during the really first cascades.

Dislocation evolution
To understand the evolution of vacancies and especially interstitials in the samples, we analyzed the length of the two commonly seen dislocations in BCC materials, the 100 and the 1/2 111 loop. In tungsten, the 1/2 111 loop is the most stable one according to DFT [57] and also the one seen in largest fractions in experiments [8] . In Fig. 5 , the total length of these two dislocation types as a function of dose is plotted. All dislocations identified by the DXA, except some very rarely seen, were of interstitial type. We observe in the results obtained with the AT-ZN potential that almost all dislocations or dislocation segments are of 1/2 111 type. In the MS potential the opposite is seen, almost entirely 100 dislocations or dislocation segments are seen. In the DND-BN potential, we initially observe 1/2 111 dislocations or segments, but as dose increases, more 100 dislocations or segments are seen. It is known that the 1/2 111 loop is very mobile, and the 100 loop is practically immobile. This can be related to the different sizes of interstitial clusters observed for the different potentials as well as the differences in evolution seen. The average dislocation density of both types at the doses higher than 0.05 dpa (above which the mean dislocation density is constant) is about 0.009 nm −2 for the AT-ZN potential, 0.016 nm −2 for the DND-BN potential and 0.015 nm −2 for the MS potential.  Looking at the same graphs and the results from the GAP relaxations, a few trends can be seen. In the cells originated from the AT-ZN potential, nothing happens, whereas for the two other potentials the length of 100 dislocation segments decrease and the length of 1/2 111 stays stable or increases.
Considering the overall dislocation evolution, we can again observe some differences between the potentials. The common thing between potentials is that almost exclusively interstitial-type dislocation loops were observed, except for some occasional vacancy loop. Snapshot of the dislocation structure for one of the runs for all potentials are found in the Supplementary material online, Sec. S3 and Fig. S3. For the AT-ZN potential, as seen in Fig. 5 , small 1/2 111 loops are formed, which grow and form a large 1/2 111 loop in the end. Some segments of 100 type are seen, mainly when the mobile loops are interacting, and occasionally small 100 loops are seen. The same evolution was seen for all three independent runs. For the DND-BN potential we observe mainly small 1/2 111 loops, with a few 100 loops, that start growing with increasing dose. However as the dose increases and the sizes of the loops grow, the 1/2 111 loops transform into 100 loops either via loop interaction or cascade overlap effects. We observed that this transformation could require higher doses in some of the runs compared to the other runs. In one of the cases, all loops were transformed into 100 loops, whereas in one it was still ongoing (see Supplementary material online for more details). At low doses for the MS potential we can see a mixture of 100 and 1/2 111 loops, and as the dose increases only the 100 loops survive and continue to grow. The same was seen in all three runs. Fig. 6 shows one representative system at ∼ 0 . 046 dpa (500 cascades) for each potential, before and after GAP relaxation. Snapshots at other doses are found in the Supplementary material online, Sec. S4 and Fig. S4. In all figures, as also seen for the dislocation segment lengths and types, when the cells are relaxed with GAP, small 100 loops and segments transform into energetically more favoured 1/2 111 dislocations, collapse and disappear, or shrink. Some of the large 100 loops are stable, as the relaxation with the GAP is short.

Swelling
Looking at the swelling in the different potentials, shown in Fig. 7 , we again see different trends for the different potentials. The AT-ZN potential results show a plateau of swelling, in correspondence with the saturation seen in defect numbers and dislocation evolution, at a bit under 0.2%. Even though the other two potentials show a very similar defect concentration evolution, their swelling are completely different, except that it is linearly increasing throughout the whole investigated dose range. The DND-BN potential show only one quarter of the swelling compared to that of the MS potential.
Looking at the swelling of the cells after the GAP relaxation, Fig. 7 markers, we can see several trends. The AT-ZN cell, with the correct dislocation type according to DFT shows almost no difference at all. The other two, which show a bit similar dislocation structures, but completely different swellings, show after the GAP relaxation very similar swelling. This swelling is three times higher than the one predicted by DND-BN and a bit lower than the one predicted by MS.

Discussion
Looking at all obtained results, we observe that the point defect evolution is tied to the dislocation evolution in the system. In pure tungsten, without any impurities that will act as traps, the 1/2 111 dislocation loop is highly mobile and the 100 dislocation loop immobile. If we do have mobile dislocations, we can at quite low doses observe that most interstitials will be in the large clusters. These large clusters are the dislocations. The same is true if we have immobile dislocations, however, the clusters will remain smaller due to the non-existent dislocation mobility at lower doses. We can also observe that the mobile dislocation loops cause a saturation of the overall number of defects, whereas if they are not mobile, the number of defects keeps growing. In the former situation, the loops will on average collect similar amounts of vacancies and interstitials, leading to a zero net increase in defects.
In the case where most dislocations are immobile, the mobility of single interstitials and vacancies created in subsequent cascades becomes more important. Single interstitials are much more mobile than vacancies and can more efficiently find and merge with the immobile loops, which will then grow and a steady increase in number of defects is seen. If we observe mobile dislocations, we also observe a defect saturation of about 0.4% defect concen-tration (at doses around 0.05 dpa), close to experimental values in tungsten and other metals [3][4][5][6] . If only or mainly immobile dislocations are formed, we see an increase in defect concentration, far exceeding these values due to the high interstitial mobility and clustering rate.
For the vacancies, which are immobile at low temperatures, we observe that mainly single and small clusters are formed, with a clustered fraction between 0.2 and 0.4. This fraction is seen to increase as the dose increases, as cascades in close vicinity to preexisting vacancies are happening, which can cause some clustering. There are about a factor 10 more single vacancies than divacancy clusters, which are seen the second most. This shows that single vacancies will account for the vast majority of all vacancytype clusters. Experimentally, it has been observed by PAS that at similar temperatures and doses as we investigate, most vacancies are indeed single vacancies or in small clusters before thermal annealing [7] .
In all simulations, we observe that the dislocation density has saturated at values between 0.01 and 0.015 nm −2 , even though different potentials show different dislocation types. This is about a factor two higher than experiments done at room temperature in the interval of 0.1 to 1 dpa, showing dislocation densities between 0.0 03 and 0.0 05 nm −2 [8] . The difference can, at least partially, be attributed to (1) we count all dislocations, even small ones that would not experimentally be visible, and (2) we account for the three-dimensional length of all dislocations and also for the dislocation junction length when several dislocations of different orientations are joined. Both of these factors will overestimate our results compared to experiments.
The short GAP relaxation did not affect the number of defects, as the annealing was not long enough for any significant defect recombination (neither was it the aim). The relaxation did, however, change the Burgers vector of most of the 100 loops. The total dislocation length was not affected dramatically, but the trend is clear. The length would most likely go down for the cases where dislocation changed Burgers vector, as the newly formed mobile 1/2 111 loops would combine with each other and absorb vacancies, leading to lower overall dislocation density, close to 0.01 nm −2 .
Looking at the swelling, again, a large difference is seen depending on defect evolution, however, there are also very large discrepancies between the EAM potentials. If the defect concentration and dislocation structure has stabilized, the volumetric swelling is also constant. If this is not the case, the swelling is still increasing. However, the two potentials showing a similar increase in defect concentration show drastically different swelling, but after a short GAP relaxation, they show a very similar result. This shows that how the interatomic potential is describing the defect volumes is crucial for obtaining a reasonable description of swelling. We also observe that the potential with the correct stable dislocation show no difference after the GAP relaxation, with a swelling of about 0.2%. This is, however, not a proof of that all defect clusters relaxation volumes found for that potential are correct, as it can be a result of error cancellation. For example, the AT-ZN potential predicts significantly lower and higher relaxation volumes than DFT for a single interstitial and vacancy, respectively, but the combined Frenkel pair relaxation volume is similar to DFT and the GAP [58] . Experimentally, the swelling has been seen to depend on many different factors, which makes a direct comparison difficult. Nevertheless, the numbers found in the literature for different doses and/or different tem peratures vary in the range of 0.2-0.4% [9,10] , close to our obtained results.
Comparing the three investigated widely used interatomic potentials, we can conclude that all aspects of the high-dose damage evolution depend strongly on the choice of potential. The DND-BN potential has previously been shown to produce similar cluster size distributions as experiments in single high-energy cascades [59] . Also previously seen for single cascades is that the AT-ZN and DND-BN show similar number or defects produced and similar fraction of Frenkel pairs surviving fraction compared to displaced atoms [14] . The MS potential shows a bit higher number in both regards [14] . This is in line with what we observe at low doses, where overlap effects are small, where the MS shows the highest rate of defect build-up in the early stages. We observed that even though two potentials showed a similar defect concentration evolution, their vacancy fraction evolution, dislocation structure evolution and swelling could differ by a factor two. Additionally, it was observed that depending on how the interatomic potential reproduced the relative stability of dislocations, the end result would either be defect saturation or a linear increase in defect concentration as a function of dose.
As our results showed a wide variety of results depending on choice of interatomic potential, the GAP simulations were carried out. We observed that the potentials with incorrect dislocation stability will yield the incorrect results, as the dislocation did change Burgers vector even in the short relaxation. Also the swelling was observed to dramatically increase and decrease respectively for these two. Our results shows that the DND-BN and MS potentials should not be used for high-dose simulations, due to their wrong defect stability and incorrect swelling. The AT-ZN potential did give very good agreement with the GAP, which indicates that among these three potentials this is the best choice for studying radiation damage buildup.

Conclusions
We have investigated the high-dose damage production and evolution in tungsten by computational means using three different EAM potentials. We found that the choice of interatomic potential will affect almost all measurable results, defect saturation, clustering and density, dislocation morphology and evolution, and swelling. We validated the results by performing short annealing and relaxation simulations with a recent GAP machine-learning potential, in order to obtain near-quantum-accurate relaxed structures. These simulations revealed that one of the three EAM potentials produces high-dose defects structures that are consistent with the GAP. This potential, the AT-ZN potential [45,46] , showed good agreement in both the dislocation structure as well as the volumetric swelling as a function of dose. Therefore, in further high-dose investigations, the Ackland-Thetford potential with the short-range modification by Zhong et al. should be preferred.
The results obtained with the AT-ZN potential also show good agreement with experiment in both saturated defect concentration (concentration of 0.4%) and dislocation density (0.01 nm −2 ). The reason for seeing this defect saturation, which was not observed for the other two potentials, is due to the incorrect stability of the two different dislocation loops in the DD-BN and MS potential. The cells obtained in the AT-ZN potential, and especially the ones from this potential that were relaxed with the quantumaccurate GAP potential, can be used in further investigations to study how these high-dose defect structures will affect other properties, for instance mechanical properties. The GAP-relaxed structures are freely available from the link found in the Supplementary material.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.