Dislocation loop formation by swift heavy ion irradiation of metals

We employ the two temperature molecular dynamics (2T-MD) model described in [30], implemented in DL_POLY_4 [31], and used previously for modelling for SHI irradiation in metals [32, 33] and semiconductors [34, 35]. In this model the electronic temperature (T_e) evolution is described by a thermal diffusion equation,

$$\begin{eqnarray}&&\frac{\partial {{T}_{\text{e}}}}{\partial t}=\underset{(a)}{{\underbrace{{\frac{{{\kappa}_{\text{e}}}\left({{T}_{\text{l}}}\right)}{{{C}_{\text{e}}}\left({{T}_{\text{e}}}\right)}}}\,}}\,{{\nabla}^{2}}{{T}_{\text{e}}}-\underset{(b)}{{\underbrace{{\frac{{{G}_{\text{e}}}\left({{T}_{\text{e}}}\right)}{{{C}_{\text{e}}}\left({{T}_{\text{e}}}\right)}}}\,}}\,\left({{T}_{\text{e}}}-{{T}_{\text{l}}}\right)+\underset{(c)}{{\underbrace{{\frac{A(r,t)}{{{C}_{\text{e}}}\left({{T}_{\text{e}}}\right)}}}\,}}\,,\end{eqnarray} \tag{ 1 }$$

with T_e dependent electronic specific heat capacity, ${{C}_{\text{e}}}\left({{T}_{\text{e}}}\right)$, and electron–phonon coupling parameter, ${{G}_{\text{e}}}\left({{T}_{\text{e}}}\right)$. The electronic thermal conductivity, ${{\kappa}_{\text{e}}}\left({{T}_{\text{l}}}\right)$, is assumed to depend only on the lattice temperature (T_l) and to be spatially independent. The first term on the right hand side of equation (1), (a), represents electronic heat diffusion and the second term, (b), represents energy exchange with the lattice via electron–phonon coupling. The third term, (c), is a source term representing the energy deposited in the electronic system by the SHI, at a distance r from the ion path at time t. We assume that the SHI travels normally (z) through the material and that the energy deposition is independent of z. A(r,t), has a Gaussian spatial dependence and an exponentially decaying time variation;

$$\begin{eqnarray}&&A(r,t)=AD(r)\alpha {{\text{e}}^{-\alpha t}}.\end{eqnarray} \tag{ 2 }$$

Here α equals $\frac{1}{\tau}$ where τ is the characteristic temporal deposition time (assumed to be 1 fs) and D(r) is a Gaussian function describing the spatial variation of the energy deposition.

$$\begin{eqnarray}&&D(r)=\frac{{{S}_{\text{e}}}}{\sqrt{2\pi {{\sigma}^{2}}}}\exp \left[-\frac{{{r}^{2}}}{2{{\sigma}^{2}}}\right].\end{eqnarray} \tag{ 3 }$$

Here r is the lateral distance from the centre of the projectile's path and σ is the standard deviation of the Gaussian spatial deposition, equal to 1 nm in our simulations. A is a normalization constant that ensures that the spatial and temporal integration equates to the energy deposited into the electronic system (the electronic stopping power, S_e, of the swift heavy ion [36]).

$$\begin{eqnarray}&&{{S}_{\text{e}}}={\int}_{t=0}^{5\tau}{\int}_{r=0}^{{{r}_{\text{max}}}}2\pi A(r,t)r\,\text{d}r\,\text{d}t.\end{eqnarray} \tag{ 4 }$$

In this study we have considered stopping powers ranging from 10 keV nm⁻¹ to 100 keV nm⁻¹, with increments of 10 keV nm⁻¹.

The structural evolution of the lattice is calculated using molecular dynamics (MD) and the energy exchange with the electronic system in equation (1) is included via a Langevin thermostat, with the local electronic temperature (T_e) used as the thermostatting temperature [30].The MD equation of motion for an atom with mass m is given by;

$$\begin{eqnarray}&&m\frac{\partial \mathbf{v}}{\partial t}=\mathbf{F}(t)-\gamma \mathbf{v}+\mathbf{\tilde{F}}(t),\end{eqnarray} \tag{ 5 }$$

where $\mathbf{v}$ is velocity of the atom at time t, $\mathbf{F}(t)$ is the deterministic force due to the interatomic potential, $\gamma \mathbf{v}$ is the frictional force (with friction coefficient γ) that represents energy loss due to electronic drag and $\mathbf{\tilde{F}}(t)$ is a stochastic force which is thermostatted at the electronic temperature. Thus energy lost/gained from the (b) term in equation (1) is gained/lost by the MD cell. Equation (1) is solved numerically by dividing the simulation cell into voxels and using a variable time step finite difference solver. The centre of the cell overlaps with the MD simulation cell to allow energy exchange between the two sub-systems (figure 1). The electronic simulation cell extends well beyond the atomistic simulation cell to enable energy to be transported away from the atomistic cell by electronic thermal conduction. The MD simulation cell has periodic boundary conditions in all 3 directions. For the electronic simulation cell we use von Neumann (zero flux) boundary conditions with zero normal derivative in the z direction to impose zero electronic heat flux in this direction. In the perpendicular (x,y) directions we aim to simulate heat transport through the bulk material beyond the electronic simulation cell, therefore we impose Robin (variable flux) boundary conditions to the electronic cells in the x,y directions. The Robin boundary conditions are a weighted combination of Dirichlet (infinite flux) and von Neumann boundary conditions and they provide a method for gradually dissipating electronic energy from the boundaries of the cell to simulate electronic thermal conductivity through the bulk. In our case we take $\frac{\partial {{T}_{\text{e}}}}{\partial t}$ to be $0.96\left({{T}_{\text{e}}}-300\right)+300$ on the cell boundaries in the x and y directions to simulate electronic energy dissipation by conduction through the bulk.

Zoom In Zoom Out Reset image size

The MD simulations cells for Fe and W have $140\times 140\times 60$ cubic unit cells (2350000 atoms). The simulation cell was subdivided into $40\times 40\times 17$ coarse-grained cells for Fe ($44\times 44\times 19$ for W). These coarse grained cells overlap with the voxels of the electronic system to enable energy exchange. The electronic system had $120\times 120\times 17$ voxels for Fe and $132\times 132\times 19$ for W, to extend beyond the boundaries of the MD cell in the x and y directions (figure 1). The fcc MD simulation cells have 2420000 atoms ($110\times 110\times 50$ cubic unit cells). The extended Finnis Sinclair interatomic potentials [37] were used for all metals and the simulation cell was equilibrated for 200 ps at 300 K and 1 atmosphere, using the Nosé Hoover thermostat and barostat, prior to initiating the SHI simulations by energy deposition to the electronic system. A time step of 1 fs was used for all MD simulations.

The temperature dependent electronic specific heat capacity, ${{C}_{\text{e}}}\left({{T}_{\text{e}}}\right)$, and electron–phonon coupling ${{G}_{\text{e}}}\left({{T}_{\text{e}}}\right)$ for all the metals (Fe, W, Cu and Ni) were calculated by Lin et al [38] using density functional theory. The ${{G}_{\text{e}}}\left({{T}_{\text{e}}}\right)$ used for W was modified slightly, based on the recent measurement of the ground state value [39]. The full T_e variation for W, and a comparison with [38], is shown in [40]. The electronic thermal conductivities were assumed to depend only on the lattice temperature and the temperature dependence was taken from Kaye and Laby [41].

The time evolution of the local lattice and electronic temperatures at the centre of the simulation cells during a 60 keV nm⁻¹ SHI irradiation is shown in figure 2. Fe has a higher electronic specific heat capacity at high T_e than W, therefore the peak electronic temperature is lower (due to (c) in equation (1)). However, the higher electron–phonon coupling strength of W results in a higher rate of energy exchange with the lattice, causing the lattice temperature of W to increase more rapidly than Fe. The electronic and lattice temperatures equalise at around 2 ps in Fe and relax to ambient temperatures at the same rate. In W the electronic temperature falls below the lattice temperature, with energy subsequently flowing from the lattice to the electrons. Both systems reach ambient temperatures after approximately 80 ps.

Zoom In Zoom Out Reset image size

The temperatures at the central region of the cells highlight how the peak temperatures evolve, however, greater understanding of the structural evolution is obtained from the temperature evolution of a cross section of the simulation cell, as shown in figure 3. Significant differences are observed between the time evolution of T_e for W and Fe. In both cases the initial electronic temperature rises in the centre of the simulation cell and spreads via diffusion and, simultaneously, energy transfers to the lattice via electron–phonon coupling. The most striking difference observed is due to the significantly higher electronic diffusivity of W. As a result of the high diffusivity, energy rapidly diffuses to the edges of the cell, resulting in a lower T_e at the centre of the cell for W than for Fe. Interestingly, despite these extreme differences in the electronic temperature profiles, the corresponding lattice temperature cross-sections, shown in figure 4, are remarkably similar. This is due to a combination of the differences in the thermal properties of the lattice and the electron–phonon coupling strength.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The melting temperature is a key property for resistance to damage by SHI irradiation. The lower melting temperature of Fe means that, despite similar lattice temperature profiles, a much wider region of the cell exceeds the melting temperature for Fe than for W, and the melting temperature is exceeded for a longer period of time. At 5 ps there is a cylindrical region of about 30 nm diameter in the Fe cell that is above the melting temperature and the corresponding diameter in the W cell is only 5 nm. Thus, it is expected that W will be significantly more resistant to SHI radiation damage than Fe.

A view along the [100] direction of the simulation cell of Fe, at various times following a 60 keV nm⁻¹ SHI irradiation event, is illustrated in figure 5. The figure shows both the time evolution of the Wigner–Seitz (WS) defects (left-hand column) and the evolution of the amorphised volume and dislocations (right-hand column). The Wigner–Seitz defects were calculated using Voronoi cell analysis [42] in OVITO [43]. The dislocations were also detected using OVITO [44]. 10 ps after irradiation a large amorphous or molten region is established at the centre of the simulation cell. As the simulation proceeds, and the lattice temperature cools, the metal starts to recrystallize and a cylindrical recrystallization front shrinks the molten region. The rapid recrystallization leaves a number of isolated vacancies in its wake, and a much lower number of interstitials. The excess of vacancies produced results in an excess of atoms at the centre of the cell, after complete recrystallization. The excess atoms create interstitials, which cluster and collapse to form dislocation loops that are aligned approximately parallel to the direction of the ion path. Two types of loops are formed, one type has Burgers vectors $\frac{1}{2}\langle 1\,1\,1\rangle$ and the second type have $\langle 1\,0\,0\rangle$ Burgers vectors, which generally have higher energies than $\frac{1}{2}\langle 1\,1\,1\rangle$ loops in bcc metals. The total length of dislocations with $b=\frac{1}{2}\langle 1\,1\,1\rangle$ is 336 Å, and for $b=\langle 1\,0\,0\rangle$ this length is 106 Å.

Zoom In Zoom Out Reset image size

We note that this defect configuration, with a halo of vacancies surrounding aligned dislocation loops, is quite distinct from that observed in cascade simulations that model low to moderate ion irradiation. In cascades, vacancies tend to reside near the ion path and dislocation loops migrate further from the path. In contrast to cascade simulations, the cross-section for nuclei interactions is extremely low at the stopping powers considered here so there are no knock-on atoms in our simulations. The residual defects are all created by the recrystallization of the molten region. Vacancies are created at the recrystallization front and these vacancies are immobile due to the high activation energy for migration. As the cylindrical molten region shrinks, there is an excess of atoms and these atoms, necessarily, form interstitials close to the ion path. For sufficiently high interstitial densities, these interstitials cluster to form pure edge (prismatic) dislocation loops. The conformation of the dislocations bear a striking resemblance to TEM characterization of polycrystalline Fe irradiated with C60 fullerenes [14], where 'two quasi-parallel dislocation lines joined near the surface' were observed. Experiments involving cascade damage in Fe [45] also create $\frac{1}{2}\langle 1\,1\,1\rangle$ and $\langle 1\,0\,0\rangle$ interstitial dislocation loops but in these cases the loops are not elongated and they are oriented in random directions.

W shows much less damage than Fe, due to the combination of higher diffusivity and higher melting temperature. The morphology of the defects created by a 60 keV nm⁻¹ ion is shown in figures 6(a) and (b). It is clear, by comparison with figure 5, that a lower number of defects was formed but that they show the same features as Fe, namely that there is a residual halo of isolated vacancies surrounding interstitial dislocation loops, which is further confirmed by the defect distribution plot shown in figure 7. The time evolution of the number of Wigner–Seitz defects is plotted for Fe and W in figure 8, which clearly shows the reduced damage in W.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

SHI irradiation experiments done on W [15] did not find any evidence of defect creation up to stopping powers of 70 keV nm⁻¹, however, these experiments were performed at low temperatures (25 K). At these temperatures the thermal conductivity of W is significantly enhanced over the room temperature value (2320 W m⁻¹ K⁻¹ [46] compared to 177 W m⁻¹ K⁻¹ at room temperature), which dissipates the deposited energy before energy transfer to the lattice. We repeated our simulations at 25 K, using the appropriate value for the thermal conductivity, and found no residual defects at an electronic stopping power of 70 keV nm⁻¹.

The number of Wigner–Seitz defects created in Fe and W as a function of stopping power is presented in figure 9. Such plots give an indication of the threshold stopping power, i.e. the stopping power below which no residual damage is created. For Fe the threshold is estimated to be 10 keV nm⁻¹ and for W the higher threshold of 30 keV nm⁻¹ is observed. The number of defects is observed to increase approximately linearly with stopping power. The variation in the length of residual dislocation lines is plotted in figures 9(b) and (c). It is clear that $\frac{1}{2}\langle 1\,1\,1\rangle$ dislocations dominate over $\langle 1\,0\,0\rangle$ dislocations and there is a general trend of increasing dislocation line length with stopping power. A decrease in length of $\frac{1}{2}\langle 1\,1\,1\rangle$ generally corresponds to an increase in the length $\langle 1\,0\,0\rangle$ dislocations. Further investigation revealed that $\langle 1\,0\,0\rangle$ dislocations were created by the overlap of two $\frac{1}{2}\langle 1\,1\,1\rangle$ loops, as previously observed in cascade simulations [47].

Zoom In Zoom Out Reset image size

A similar investigation was performed for the fcc metals, Cu and Ni. As expected from earlier experiments [14, 15], these metals proved to be very resistant to damage, with no residual defects being observed up to stopping powers of 60 keV nm⁻¹ in Ni and a very small number (32) of defects in Cu at this stopping power. However the reason for the resistance to damage was different for the two metals. In Cu it is due to the high electronic thermal conductivity, which enables energy to diffuse rapidly away from the ion path. The rapid diffusion of electronic energy ensures that energy is transported away from the track before the electron–phonon coupling has time to transfer a significant amount of energy to the lattice and melting is limited. In Ni, however, it is the decrease in the electron–phonon coupling strength with increasing electronic temperature that is primarily responsible for decreased energy transfer to the lattice. The electronic temperature dependence of the electron–phonon coupling constant for Cu and Ni [38] is shown in figure 10, which illustrates the very low coupling for Ni at high T_e. The lattice temperature evolutions for 60 keV nm⁻¹ simulations for Cu and Ni are shown in figure 11. There is no detectable rise in lattice temperature for Ni due to the weak electron–phonon coupling at high electronic temperatures. The lattice temperature of Cu rises above the melting temperature in a narrow region (radius 4 nm) and results in limited melting and the creation a small number of defects. As with Fe and W, the vacancies are isolated and the interstitials form very small clusters. These results highlight that the electronic temperature dependence of the material parameters (electronic diffusivity and electron–phonon coupling) can play a crucial role in the resistance of metals to damage by SHI irradiation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size