Effects of electron–phonon coupling on damage accumulation in molecular dynamics simulations of irradiated nickel

ABSTRACT

The role of the electronic system in high energy displacement cascades is explored. The energy exchange between the electronic and the atomic subsystem is described by the electron–phonon coupling. The electronic effects on the damage accumulation due to 100 keV Ni ion cascades in nickel, a prototype system to a large group of nickel-based high entropy alloys, are investigated for overlapping cascades. It is shown that the energy exchange between the two subsystems affects microstructure evolution, resulting in the formation of smaller clusters and more isolated defects. This effect is more significant for the vacancy cluster formation and size distribution.

IMPACT STATEMENT

The effect of electron–phonon coupling on the damage accumulation due to ion irradiation is investigated for the first time, with results showing significant effects on the vacancy cluster formation.

GRAPHICAL ABSTRACT

KEYWORDS:

• Molecular dynamics
• two-temperature model
• electronic effects
• electron–phonon coupling
• Nickel

Molecular dynamics (MD) methods are traditionally used to describe the dynamic evolution of large systems. Within the MD approach, the motion of the atoms of the system is described by the Newtonian equations of motion, and the presence of electronic subsystem is incorporated in the interatomic potentials. However, this adiabatic approach is no longer applicable when atomic energies start to be comparable with typical energies of electronic excitations. The latter takes place in displacement cascades produced by high energy primary knock-on atoms (PKA). When a fast moving particle interacts with a target material, its energy is partitioned between the atoms and the electrons of the target, resulting in atomic displacements and electron excitations. The energy is dissipated in each subsystem, and depending on the temperature difference between them, energy is exchanged further between the electrons and the atoms, with the electrons acting as a heat source or sink. The energy transfer to the atomic subsystem, known as nuclear or elastic energy loss, creates ballistic displacements and collision cascades, which results in single defects and defect clusters. The effect of the electronic energy loss (energy transfer to the electrons) is twofold: (i) energy is removed from the ionic system due to inelastic electronic collisions (electronic stopping) and electron–phonon (e–ph) coupling and (ii) this energy is transferred to the electronic subsystem and diffuses further, and part of it is redeposited to the ions via e–ph coupling.

In atomistic modeling of ion irradiation, these energy loss mechanisms can be described with the two temperature (2T-MD) model [1,2], where the atomic and electronic systems are two coupled subsystems and exchange energy via the electron–phonon coupling. MD simulations of single ion events that include the 2T-MD model have shown that the e–ph interactions affect the damage morphology [1–9]. The e–ph coupling in Ni ion irradiation of nickel and fcc nickel-based concentrated nickel alloys results in a smaller number of defects, larger number of isolated defects and smaller clusters. Thus accurate accounting for the energy dissipation due to atomic interaction with the electronic system is necessary, as the damage morphology is an important feature of the response of materials to radiation. The cluster size can affect the damage evolution at longer times. Local damage can induce local strains that affect the material's behavior at longer times.

In this work, we investigate the effect of the e–ph coupling on the damage accumulation due to 100 keV Ni ion cascades in Ni. Nickel is a model system, as it is the main component of fcc nickel-based concentrated solid solution alloys. These alloys, consisting of two or more components (high entropy alloys), are promising materials for nuclear applications because of their good mechanical, thermal and electric properties [10–16]. Due to their chemical complexity, they exhibit different e–ph coupling parameter and electronic thermal conductivity, parameters which affect the damage production and morphology [17]. We study the effect of increasing the number of cascades with and without the e–ph coupling. The cascade events, a total of 20 for each case, are initiated close to each other so the cascades can interact and resulting damage can propagate through the system. The time between the events, about 60 ps, is sufficient for the system to relax.

In MD, the 2T-MD model is implemented in two steps: (i) by including the electronic stopping and (ii) by including the e–ph interactions. The electronic stopping is included as a friction term ${\gamma }_s$ in the equation of motion (Equation 1(1) $m_i\frac{\mathrm{\partial }{\mathbf{v}}_i}{\mathrm{\partial }t}={\mathbf{F}}_i\left(t\right)-({\gamma }_s+{\gamma }_p){\mathbf{v}}_i+\tilde{\mathbf{F}}\left(t\right).$(1) ), which removes energy by slowing down highly energetic atoms, with velocity higher than a specified cut-off. The e–ph interactions are included with an extra friction term ${\gamma }_p$, and a stochastic $\tilde{\mathbf{F}}\left(t\right)$ force that returns energy from the electronic to the atomic system. The stochastic force depends on the electronic temperature via the fluctuation–dissipation theorem $⟨\tilde{\mathbf{F}}\left(t^{\prime }\right)\cdot \tilde{\mathbf{F}}\left(t\right)⟩=2k_B{T}_e{\gamma }_s\delta (t^{\prime }-t)$, whose evolution is described via the heat diffusion equation (Equation 2(2) $Ce\frac{\mathrm{\partial }{T}_e}{\mathrm{\partial }t}=\mathrm{\nabla }\left({\kappa }_e\mathrm{\nabla }{T}_e\right)-g_p(T_e-T_{\alpha })+g_s{T}_{\alpha }^{\prime },$(2) ). As seen in Equation 1(1) $m_i\frac{\mathrm{\partial }{\mathbf{v}}_i}{\mathrm{\partial }t}={\mathbf{F}}_i\left(t\right)-({\gamma }_s+{\gamma }_p){\mathbf{v}}_i+\tilde{\mathbf{F}}\left(t\right).$(1) , the term $g_p\cdot (T_e-T_{\alpha })$ describes energy exchange between the two subsystems, which depends on the local temperature difference. (1) $m_i\frac{\mathrm{\partial }{\mathbf{v}}_i}{\mathrm{\partial }t}={\mathbf{F}}_i\left(t\right)-({\gamma }_s+{\gamma }_p){\mathbf{v}}_i+\tilde{\mathbf{F}}\left(t\right).$(1) In order to use the 2T-MD model and include the electronic subsystem in MD, we need to include properties of the electrons, such as the electronic specific heat $C_e$ and the electronic thermal conductivity $k_e$, as well as the e–ph coupling constant $g_p$. In Equation  (2(2) $Ce\frac{\mathrm{\partial }{T}_e}{\mathrm{\partial }t}=\mathrm{\nabla }\left({\kappa }_e\mathrm{\nabla }{T}_e\right)-g_p(T_e-T_{\alpha })+g_s{T}_{\alpha }^{\prime },$(2) ), $g_s$ accounts for the electronic stopping and $T_{\alpha }$ is the atomic temperature. $T_{\alpha }^{\prime }$ has units of temperature and is calculated from the average kinetic energy of the subset of atoms that slow down due to the electronic stopping [1]. (2) $Ce\frac{\mathrm{\partial }{T}_e}{\mathrm{\partial }t}=\mathrm{\nabla }\left({\kappa }_e\mathrm{\nabla }{T}_e\right)-g_p(T_e-T_{\alpha })+g_s{T}_{\alpha }^{\prime },$(2) ${\gamma }_s$ and ${\gamma }_p$ in Equation  (1(1) $m_i\frac{\mathrm{\partial }{\mathbf{v}}_i}{\mathrm{\partial }t}={\mathbf{F}}_i\left(t\right)-({\gamma }_s+{\gamma }_p){\mathbf{v}}_i+\tilde{\mathbf{F}}\left(t\right).$(1) ) are related to the electronic stopping constant $g_s$ and the e–ph coupling constant $g_p$ in Equation  (2(2) $Ce\frac{\mathrm{\partial }{T}_e}{\mathrm{\partial }t}=\mathrm{\nabla }\left({\kappa }_e\mathrm{\nabla }{T}_e\right)-g_p(T_e-T_{\alpha })+g_s{T}_{\alpha }^{\prime },$(2) ) via the relations: (3) $g_p=\frac{3N{k}_B{\gamma }_p}{\mathrm{\Delta }V{m}_i}\mathrm{a}\mathrm{n}\mathrm{d}{g}_s=\frac{3N^{\prime }{k}_B{\gamma }_s}{\mathrm{\Delta }V{m}_i},$(3) where N is the number of atoms in volume $\mathrm{\Delta }V$, $k_{\mathrm{B}}$ the Boltzmann constant, m the atom mass and $N^{\prime }$ is the number of atoms with velocities larger than $v_c$ within volume $\mathrm{\Delta }V$.

We use the DL_POLY MD code [18] for the simulations. The systems consist of about 7 million atoms and the MD box has 426 Å length. The MD box has periodic boundary conditions. The primary knock-on atoms are initiated near the $(-x,-y,-z)$ corner of the MD box, directed towards the center of the box (0,0,0). We simulate twenty different recoil directions, with the same starting point for the electronic stopping and 2T-MD cascades, within $\pm 2$ Å. We use an Embedded Atom Method (EAM) type potential by Bonny et al [19]. The systems were initially equilibrated in the constant pressure-constant temperature ensemble (NPT) at 300 K, with a timestep of 1 fs. For the cascades simulation, the atomic velocities within 8 Å of the MD box boundaries are scaled according to a Gaussian distribution to emulate the bulk. We use a variable timestep in cascade simulations to describe the cascade evolution. The electronic stopping mechanism is applied to fast moving atoms with velocities larger than the cut-off v${}_c$. This is equal to 54 Å/ps and corresponds to energy twice the cohesive energy of the system [1,20], which is about 8.8 eV. The corresponding relaxation time for electronic stopping is ${\tau }_{\mathrm{s}}=m_i/{\gamma }_{\mathrm{s}}$ [1] where ${\gamma }_s$ is 0.6 ps${}^{-1}$ and it is calculated using SRIM tables [21], as described in [9]. The e–ph coupling is activated at 0.2 ps simulation time, determined by the thermalization time in cascades with electronic stopping only. The e–ph coupling parameter $g_{\mathrm{p}}$ is $10.85\times {10}^{17}$ W m ${}^{-3}$ K${}^{-1}$, obtained using electronic structure results [23]. The timescale for the energy loss due to the e–ph friction term is ${\tau }_{\mathrm{p}}=m_i/{\gamma }_{\mathrm{p}}$ [1], and the e–ph coupling friction term ${\gamma }_p$ = 0.287 $p{s}^{-1}$. The specific heat $C_e$ was obtained from electronic structure calculations [23] as a function of the electronic temperature. The experimental thermal conductivity ${\kappa }_{\mathrm{e}}$ is 88 W m${}^{-1}$ K${}^{-1}$ [22]. The electronic grid extends beyond the MD box and it has Robin boundary conditions, to represent thermal electronic conduction into the bulk, allowing the electronic temperature to converge towards the initial system temperature.

Figure 1 shows the surviving number of Frenkel pairs (FPs) after each successive cascade event, with and without e–ph coupling. The defects were identified with the Wigner–Seitz method. As expected, in both cases the surviving damage increases after each cascade event. However, for each event, there is a smaller number of defects when the e–ph coupling is taken onto consideration, compared to the cascades where it is ignored. The difference in the contribution to the damage becomes more obvious and significant as the damage level increases. The damage accumulation clearly shows the evidence of approaching saturation due to the creation of defects within the recombination volume of existing defects.

Figure 1. Number of Frenkel pairs for 20 100 keV Ni cascade events in Ni, with and without the e–ph coupling.

In 100 keV Ni ion irradiation of Ni, about 30 % of the projectile's energy is transferred to the electrons, according to calculations with the SRIM code. Figure 2 shows the atomic and electronic temperatures along the y-direction in the center of the system's xz-plane, for a representative 2T-MD cascade event. Each curve corresponds to the temperature of a voxel of 18.5 Å size, along the y direction. The temperature is higher at the cascade evolution area, while away from the cascade, the temperature is lower. As seen here, the local electronic and atomic temperatures can be very high in such high energy events. Depending on the local temperature difference between the two subsystems, energy can flow both ways. At the start of the simulation, the energy from the primary knock-on atoms dissipates, causing an increase of the atomic temperature. The electronic stopping removes energy from the cascade by slowing down fast moving atoms. This energy is transferred to the electronic subsystem, increasing the temperature of the electrons. In addition to the electronic stopping, the friction term due to the e-ph interactions removes energy from the atomic to the electronic system. The energy dissipates further in the electronic subsystem, and depending on whether the temperature is locally higher than that in the atomic system, energy returns to the atoms. This energy feedback contributes to the recombination of displaced atoms produced both due to the current event and to cascades that occurred earlier. Consequently, the number of surviving defects when the e–ph coupling is included in the simulations is smaller.

Figure 2. (a) Atomic and (b) electronic temperature temporal evolution for a representative 100 keV Ni 2T-MD cascade in nickel. The curves correspond to the local temperature of voxels in the atomic and electronic grid, respectively, along the y direction of the MD box. Each voxel has length $y_2-y_1$ ≈ 18.5 Å. In the legend, we show the $y_2$ value of each voxel.

Figure 3 shows the distribution of interstitial and vacancy clusters found after the each cascade event, with and without the e–ph coupling. The clusters were identified using the next nearest neighbor criterion, and we have sampled them as single interstitials or vacancies, small (size 2–20), medium (size 21–50) and large (size larger than 50) clusters. For single defects, the e–ph coupling results in a larger number of isolated vacancies, with this effect being more profound as the number of cascade events increases (Figure 3 f), while the e–ph coupling effect on the single interstitials is smaller (Figure 3 a). For clusters consisting of 2–20 interstitials or vacancies, the effects of the energy feed-back to the lattice are more significant: the e–ph interactions affect the cluster formation, resulting in the formation of a larger number of small interstitial (Figure 3b) and vacancy clusters (Figure 3 g). For medium size clusters (size of 21–50 interstitials/vacancies), the effects on the interstitial clusters are less significant, while there is a strong effect in the number of vacancy clusters (Figure 3 h), which is more profound as the number of events increases. The same important effect on the vacancy clusters is also observed for large size (size larger than 50) clusters (Figure 3i and j). The e–ph interactions also affect the formation of large interstitial clusters (size 51–100) (Figure 3 d), while for interstitial clusters with size larger than 101 the effect is not significant. Overall, in this figure it can be seen that, when energy is allowed to return to the atomic subsystem due to the ionization of the electrons, a larger number of small vacancy clusters is formed. Therefore, less vacancies are available for the formation of large clusters, meaning dramatic suppression of the formation of large vacancy clusters. In Table 1, the largest interstitial and vacancy cluster found in each simulation are listed. Both the interstitial and vacancy clusters found in the 2T-MD model cascades are significantly smaller compared to the electronic stopping only cascades. As seen in this table, the largest interstitial cluster found in each 2T-MD cascade has size about half of the largest cluster found in the corresponding cascade event when the e–ph coupling is ignored. The difference is more dramatic for the vacancy clusters, where the size of the largest clusters when the e–ph coupling is considered is some tens of vacancies, while when the e–ph coupling is ignored the largest vacancy clusters consist of hundreds of vacancies. These statistics indicate that ignoring the e–ph coupling in simulations of irradiation in metals results in overestimating the size of the defect clusters. This is an important factor to consider, as larger clusters evolve differently than smaller defect clusters at longer times.

Figure 3. (a)–(e) Interstitial cluster count found in each cascade event, shown according to their size. (f)–(j) Vacancy clusters found in each cascade event, shown according to their size. The green bars represent simulations without the e–ph coupling, while the magenta bars represent the cascades with e-ph coupling taken into account.

Table 1. Largest interstitial and vacancy cluster found in each irradiation event in nickel, for both irradiation models.

Event	Largest int. cluster	Largest vac. cluster	Largest int. cluster	Largest vac. cluster
	Elec. Stopping		2T-MD
1	66	21	32	13
2	67	22	34	14
3	69	60	35	20
4	70	63	36	21
5	238	603	143	51
6	264	410	120	65
7	288	259	146	63
8	322	415	149	118
9	321	413	146	91
10	266	176	139	61
11	244	188	145	31
12	273	188	176	32
13	266	283	177	32
14	227	526	176	35
15	417	204	178	46
16	434	363	186	61
17	412	279	183	61
18	371	529	189	79
19	370	515	184	61
20	409	448	184	61

Figure 4 shows the final damage, after 20 events for each irradiation model. As seen here, the e–ph coupling affects the microstructure, resulting in more isolated defects and smaller clusters. The e–ph coupling significantly affects the dislocations that are formed during irradiation, resulting in a much smaller dislocation network, as shown in Figure 4 (c) and (d). The dislocations shown here are calculated with Ovito software [24,25], with a total of 503 and 320 dislocation segments for the electronic stopping and the 2T-MD model cascades, respectively. The green color corresponds to Shockley dislocations, the turquoise to Frank (interstitial dislocations), stair-rod are shown in magenta, Hirth in yellow and perfect dislocations in blue. The purple structures represent defect clusters.

Figure 4. Final damage after 20 cascades (a) with electronic stopping and (b) with the 2T-MD model. Interstitials are shown in green and vacancies in orange. (c) and (d) show the dislocation network formed in each irradiation model, calculated with Ovito software [24]. The green color corresponds to Shockley dislocations, the turquoise to Frank, the magenta to stair-rod type, the yellow to Hirth and the blue to perfect dislocations. In these figures, the defect clusters are shown in purple. Some SFTs and partials SFTs are shown in orange.

In conclusion, the energy exchange between the atomic and electronic subsystems in high energy irradiation has significant impacts on the microstructure, resulting in a larger number of isolated defects and smaller size clusters. These effects are more profound for increasing number of irradiation events, particularly for vacancy clusters, the size of which is dramatically smaller. These microstructure features play important roles in the material's response to radiation. Taking the e–ph coupling into account is necessary in order to fundamentally understand primary radiation damage, as a large part of the projectile's energy is transferred to the electrons. The 2T-MD model is a more realistic approach to large energy irradiation, as includes both the energy transfer to the nuclei and to the electrons. Our results indicate that properties of the electrons such as the e–ph coupling and the thermal conductivity should be taken into account in the design of new materials for radiation environments. Further investigation of the e–ph coupling effects in microstructure and damage accumulation is needed in order to fundamentally understand the primary stages of radiation effects. Investigation of the resulting structures at longer times using accelerating methods is needed in order to assess these effects at longer time, relevant to experimental timescales.

Acknowledgments

This manuscript has been authored by UT-Battelle, LLC under Contract No. DE-AC05-00OR22725 with the U.S. Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan).

Disclosure statement

No potential conflict of interest was reported by the authors.

ORCID

Eva Zarkadoula http://orcid.org/0000-0002-6886-9664

German Samolyuk http://orcid.org/0000-0001-6877-8255

William J. Weber http://orcid.org/0000-0002-9017-7365

Additional information

Funding

This work was supported by Energy Dissipation to Defect Evolution (EDDE), an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under Contract DE-AC05-00OR22725. The simulation used resources of the National Energy Research Scientific Computing Center, supported by the Office of Science, US Department of Energy, under Contract No. DEAC02-05CH11231.