Jiting Tian(田繼挺)

Institute of Nuclear Physics and Chemistry,China Academy of Engineering Physics,Mianyang 621900,China

Keywords: displacement cascades,molecular dynamics,potential energy landscape,metals

1. Introduction

Understanding and predicting materials’ response to radiation are crucial for many important applications like nuclear power safety, ion implantations in semiconductors, device performance in space, as well as radiotherapy. From a physical viewpoint, irradiation effects can be regarded as an energy deposition into the materials through strong interactions between the energetic particles and the substance, followed by a complex relaxation of the materials under a certain environment, finally leading to some obvious changes in the materials’ structures and/or properties. Despite decades of intense investigations,[1–12]a full understanding of radiation damage in solids is still inaccessible. This is mainly because the material’s evolution during and after irradiation is a typical multi-scale phenomenon expanding from nm and ps to m and year,[11,13,14]making the dynamics of the process quite difficult to capture and understand. Several theoretical models have been proposed to describe irradiationinduced damage in solids,such as the KP model,[15]the thermal spike model,[16,17]and the shock waves model.[18–20]However, most of them focus on only one particular aspect of the multi-scale process and the remaining contents are addressed in a quite empirical way or even completely ignored.Overall,a comprehensive and uniform framework to describe and understand radiation damage is still absent.

Taking another perspective, properties of a condensed matter system (such as a metal crystal) mainly depend on its potential function Φ(r1,r2,...,rN),which can be viewed as a topographic map (usually called potential energy landscape,PEL)[21–23]depicting the ‘a(chǎn)ltitude’ at any ‘position’ R =(r1,r2,...,rN) in the 3N-dimensional configuration space.The PEL of a many-body system is mainly determined by the atomic interactions,and usually consists of energy basins(corresponding to stable or metastable states)separated by activation energy barriers. At 0 K, the system must locate at the bottom of one basin (i.e., an energy minima) so that the total force on every atom is zero. At low temperatures (far below the melting point),the system usually spends a lot of time wandering around in one energy basin (corresponding to the small thermal vibrations of the atoms away from their equilibrium positions)so that the material’s properties do not change.However, when disturbed by an external stimulus (such as heating or strain), the system may obtain enough energy to cross the surrounding barriers to reach new regions in the PEL,inducing changes in the material’s structure and/or properties.This energy landscape framework has been widely applied in many important topics,such as protein folding,[24–26]the mechanical properties of solids,[27–29]and the dynamics of glassy transitions.[21–23]

As irradiation can also be regarded as a stimulus,we can expect that analyzing the evolution of the system in the PEL under irradiation may provide some fundamental insights on the dynamic process of radiation damage. Actually,a very recent study[12]has applied this idea to the irradiation-induced disordering of quartz,finding that the damage buildup is governed by the topography of the enthalpy landscape (enthalpy replaces potential energy when the system’s volume is changing). However, while defect accumulation and saturation in a semiconductor has been clearly elucidated by this inspiring work, the dynamic process of a single irradiation event(i.e., a displacement cascade induced by an incident particle) is still ambiguous. Besides, metals response quite differently compared with semiconductors upon irradiation,[3,4]and the understanding of radiation damage in metals from the PEL perspective is still lacking. As metals and metallic alloys are widely used as structural materials in nuclear power systems,[9]a fundamental understanding of irradiation effects in them is of great importance.

In this paper,using hcp Zr as a model material,we utilize molecular dynamics (MD) method to simulate single irradiation events as well as continuous irradiation in a metal crystal,and interpret the results from the PEL perspective.We find that prolonged irradiation can be understood as a discrete Brownian motion in the PEL, with one single irradiation event as a step whose trajectory contains an uphill journey and a downhill one. The dynamics of the trajectories and the topographic features of the PEL can explain some important results in the research of radiation damage in solids. Based on these results,we discuss the advantages of the PEL method as well as its importance in future work.

2. Methods

2.1. Simulating single cascades

All MD simulations in the present study are performed using the LAMMPS code.[30]We adopt a well-established protocol[4,12,31–33]to simulate 5-keV single cascades in hcp Zr at three different temperatures (1 K, 300 K, and 600 K).Its details have been provided in our previous works,[31–33]so here we only give a rough description. A size of 230 ?A×230 ?A×230 ?A hcp Zr crystal (containing 532224 atoms) is created and relaxed at a certain temperature for 5 ps using an NVT thermostat.[34]Then a primary knock-on atom(PKA,also called recoil)is randomly chosen from a 2-?A-thick spherical shell located at the center of the sample and its velocity is initialized as 1026.2 ?A/ps(corresponding to a kinetic energy of 5 keV)towards the sphere center.According to this kinetic energy,we set the radius of the spherical shell to be 50 ?A,which is able to ensure that the triggered cascade mainly expands at the center,as certified by our test simulations.

During a collision cascade, like previous studies,[31–33]we apply an NVT thermostat[34]on six 10-?A-thick atomic layers at the boundaries to absorb the heat. The remaining region is treated in the NVE ensemble to model realistic atomic collisions and energy transfer. The dynamics of the cascade is simulated for 42.2 ps(82.2 ps for a typical cascade event,see Subsection 3.1) with a timestep of 0.1 fs for the first 0.2 ps and another of 1 fs for the left. Our results will show that this time(42.2 ps)is long enough for the sample to cool down and reach a stable configuration. For each simulation run,(instantaneous)configurations at different moments in the whole process are stored for a later analysis(see Subsection 2.3). As the evolution of a cascade is stochastic in nature, 20 simulations (i.e., 20 different PKAs) have been performed for each temperature to obtain statistical results.

2.2. Simulating cascade overlap

In experiments, materials usually suffer from continuous irradiation, which leads to defect accumulation and reaction in a local region in the sample. To model this phenomenon, we should simulate overlapping cascades rather than a single cascade in the crystal.This kind of simulation has been performed for several materials,like SiC,[35]bcc-Fe,[36]fcc-Ni,[37]concentrated solid solution alloys (CSAs),[38]and quartz,[12]while a similar study for hcp Zr is still absent.

The basic MD setup for cascade overlap is similar with that for single cascades,so here we only describe the different aspects. More specially,an atom is randomly chosen from the whole sample (no longer from the spherical shell mentioned in Subsection 2.1) as the PKA and then is given a velocity(still 1026.2 ?A/ps) towards a random direction (no longer towards the sample center),as shown in Fig.1. Then we assume that the cascade center is 50 ?A away from the PKA’s original position along the velocity direction. As we mentioned in Subsection 2.1, test simulations show that this assumption is reasonable for 5-keV cascades. Based on this assumption,three 20-?A-thick atomic layers(the blue part in Fig.1)located 105 ?A(～half of the sample size)away from the cascade center are set as the constant-temperature walls(CTWs)using an NVT thermostat[34]to dissipate the heat.

Each cascade event is simulated for 42.2 ps with the same timesteps used for single cascades. When one cascade event ends, the (instantaneous) configuration is stored for the later analysis (see Subsection 2.3). Then a new PKA is randomly chosen in the same sample and the process mentioned above is repeated until a total of 1000 recoils for each temperature are simulated. Using the NRT equation[39]we find that one 5-keV cascade in an hcp Zr crystal containing 532224 atoms corresponds to a dose of ～0.0001 dpa (we use 40 eV as the displacement threshold energy for hcp Zr), thus 1000 events correspond to ～0.1 dpa. Our results below will show that at this dose the damage in the sample is already saturated.

Note that unlike simulations of single cascade events,here the CTWs are not fixed at the boundaries but are adaptively established according to the cascade center in each event.Benefitted from the periodic boundaries,the cascade always occurs at the center of the region surrounded by the three CTWs,no matter where the PKA is initiated.This kind of consecutive totally random PKA initializations(for both position and velocity direction)enable us to model homogeneous irradiation over the whole sample. However, we must stress that this method strongly depends on the assumption that the cascade center is approximately located along the PKA velocity direction with a certain distance(50 ?A for 5-keV PKAs in this work) away from the PKA’s original position, which is only reasonable for low-energy cascades. For a high-energy recoil(such as 20 keV or beyond[33]), the cascade region is often quite irregular and the cascade center is quite difficult to predict,thus our method would not work. It is worthy to mention that our procedure to simulate continuous irradiation is similar with those adopted in the literature,[12,28]although the specific techniques are different.

At last,we should point out that our simulations have several aspects quite different from experiments. First, in realistic situations the PKA energy is in a wide spectrum (usually from eV to MeV) instead of a single value, so some phenomena induced by high-energy PKAs (such as large defect cluster formation[33]) are missing in our simulations. However, here we focus on a fundamental understanding of the irradiation process rather than exact comparisons with experiments, thus 5 keV should be a good choice for both inducing a typical displacement cascade and saving a lot of computational time. Second, the dose rate in our simulations(0.0001/42.2 dpa/ps=2.37×106dpa/s)is extremely high comparing with that in experiments (usually on the order of 10?3dpa/s[40,41]), as the thermally activated relaxation between cascades in experiments has been ignored in the present simulations for the well-known fact that currently MD simulations can only reach a time scale of ps-ns.However,it is worth noting that the previous studies using similar methods to simulate the damage buildup in metals have shown a fine agreement with experimental results(measured by RBS/C),[37]implying that the difference in the dose rates between simulations and experiments has no significant influence, at least on the damage characterized by RBS/C.

2.3. Data analysis

Two traditional methods are utilized to analyze the defects produced in the sample during and after a cascade. The Wigner–Seitz method is used to identify and count the point defects(i.e.,self-interstitial atoms(SIAs)and vacancies). The common neighbor analysis(CNA)method is adopted to identify the lattice type(cut-off for hcp Zr is 3.887 ?A).If an atom does not belong to the hcp type, we regard it as a disordered atom. To quantitatively describe the irradiation-induced disorder in a crystal,we define the disorder fraction(fd)as

where Ndisrefers to the number of disordered atoms, equaling the number of all atoms N reduced by the number of hcp atoms Nhcpwhich can be counted by the CNA method. Both of the two methods mentioned above have already been implemented into LAMMPS as well as the visualization software OVITO.[42]

We mainly use potential energy (PE) per atom and displacement(D)in the 3N-dimension configuration space to describe the system’s evolution in the energy landscape. PE is easy to calculate in MD simulations and D can be obtained using the following equation

where ri,0and ri,tare the initial and the instantaneous position of the i-th atom,respectively. For a given(instantaneous)configuration of the sample stored during a simulation run (see Subsections 2.1 and 2.2), we perform an energy minimization(using the conjugate gradient algorithm[43])to obtain the ground-state configuration(i.e.,inherent structure),which corresponds to a local minima in the energy landscape. PE and D are also calculated for these ground-state configurations. To distinguish the results calculated before and after a minimization, we add I (for instantaneous) or G (for ground-state) to the variables. For example,IPE and GPE represent the instantaneous potential energy (IPE) and the ground-state potential energy (GPE), respectively. From a physical viewpoint, IPE represents the potential energy of the real system at a finite temperature and contains contributions from the thermal vibrations which are stochastic in nature and strongly depends on the environment temperature. On the other hand, GPE is the potential energy of the inherent structure(the equilibrium positions of the vibrating atoms)where the thermal noises are already removed.

3. Results

3.1. A typical cascade event

First we use a typical irradiation event induced by a 5-keV recoil at 1 K to show the evolution of the cascade process.This event is simulated for 82.2 ps to ensure that the sample completely cools down at the end of the simulation. Figure 2(a)exhibits how the number of defects changes along with time.We can see that either Ndef(number of point defects) or Ndis(number of disordered atoms)increases rapidly at first,reaches its peak value at 0.5 ps–1.0 ps,and then decreases with a comparatively slow rate, reaching a steady value at around 10 ps with no changes any more. The inset snapshots in Fig.2(a)demonstrate the disordered atoms in the sample at three moments (colors represent IPE, see the caption of Fig.2). We can see that a nanoscale disordered region with an approximatively sphere shape evolves at the center of the sample and then disappears at the end of the irradiation event, leaving a few dispersed disordered atoms(corresponding to several point defects).Previous studies[1–4]have shown that this evolution of a disordered region in an irradiation event actually corresponds to a local short-time melting(usually called a heat spike).This kind of ps-scale melting-and-cooling process is a typical feature of a collision cascade and has been extensively studied in the past.[1–4]

Fig.2. Evolution of a typical displacement cascade induced by a 5-keV recoil at 1 K in hcp Zr: (a) numbers of defects and disordered atoms, (b)potential energy, and (c) displacement in the configuration space. The inset snapshots in(a)show the disordered atoms with colors from blue to red presenting the potential energies from low values to high ones. The grey rectangle in panels (b) and (c) represents the peak moment of the damage during the cascade,and the grey dashed line refers to the moment at which the system reaches a stable state.

Now we adopt the energy landscape perspective to understand the irradiation process. Figure 2(b) depicts the evolution of GPE and IPE during a cascade event. We can see that the change in GPE (on the order of 10?4eV) is significantly smaller than that in IPE (on the order of 10?3eV). Note that we use a double-Y plot to show the details of the two variables at the same time. It is also observed that like Ndefand Ndis,the two PEs also increase rapidly at first and then decrease, indicating that the system at first moves to a high-altitude region in the PEL and then slides down to a low-altitude one. While the IPE changes quite smoothly,the GPE shows some fluctuations during the cascade,implying that the system is jumping between basins in the PEL,as different GPEs indicate different basins. After ～10 ps,the GPE stays unchanged,meaning that the system is already trapped in a basin,although the IPE still decreases slowly due to the gradually-declining thermal vibrations during the cooling process. Note that the moment at which the GPE becomes stable(～10 ps)is in exact accordance with those for Ndefand Ndis(see the grey dashed line in Fig.3),indicating that GPE is good at describing the damage evolution in the sample. Figure 2(c)demonstrates the 3Ndimensional displacement D defined by Eq. (2) as a function of time. The ground-state displacement (GD) is close to the instantaneous one, as the two corresponding positions of the system in the multidimensional configuration space are in a same basin in the PEL(the ground-state configuration is at the bottom of the basin while the instantaneous one usually not).We can see that after an irradiation event,the system has gone～200 ?A far away from its initial position in the configuration space. The grey rectangle at 0.3 ps–0.5 ps(when the GPE peaks)roughly separates the system’s trajectory into two parts:the uphill one and the downhill one. It can be observed that the two journeys are similar in length(～100 ?A),but the time the system spends on them are quite different. While the system only spends ～0.4 ps to reach the high-elevation region,it spends ～10 ps to slide down to the low-elevation region.

To examine the details of the fluctuations in Fig.2(b),we sample more data points in two periods (0 ps–0.2 ps and 0.2 ps–2.2 ps)and plot them in Fig.3(each period containing 1000 equidistant points). From Fig.3(a), it is observed that the IPE features several narrow peaks in 0-0.1ps, indicating strong short-distance atomic collisions(usually termed ballistic phase). On the contrary, in the first 0.1 ps the GPE shows a ladder shape,implying occasional jumps from low basins to higher ones in the PEL.After 0.1 ps,the peaks in the IPE become wider and less-frequent,meaning that drastic atomic collisions gradually disappear. The GPE is still like a ladder,but the steps become shorter and shorter,indicating that the transitions between basins become more and more frequent. From Fig.3(b),we can see that the IPE evolves quite smoothly,but the GPE shows highly frequent fluctuations,meaning that the system is constantly exploring a vast number of basins in the PEL,which is a typical feature of a liquid.[22]For the cascade process after 2.2 ps, we can expect that as the sample cools down,the system will have less and less kinetic energy to cross the barriers,so it will finally be trapped in a low-energy basin,as implied by the steady value of GPE in Fig.2(b). As defects are already introduced into the sample (7 point defects in this event),the final state should be a little higher(actually～5.7×10?5eV/atom in this event)than the initial state(the pristine crystal)in GPE.

Fig.3. Detailed evolution of the potential energy of the system during a typical 5-keV cascade in hcp Zr at 1 K: (a) 0 ps–0.2 ps and (b) 0.2 ps–2.2 ps.1000 data points are sampled for each period. Note that we use double-Y to show the two different scales.

Overall, using a typical event as an example, we show that from the PEL perspective,an irradiation event can be understood as a trajectory of the system in the PEL triggered by an energy deposition on an atom (the PKA). The trajectory comprises a fast uphill journey and a slow downhill one with a similar length.When the system is close to the damage peak,it exhibits a liquid feature associated with a local melting in the sample. In the next section, we will investigate how ambient temperature affects the system’s trajectory in the PEL.

3.2. Effect of temperature

It is well-known that ambient temperature can have a significant influence on the irradiation-induced damage in metals.[8,44]Here we plan to understand this phenomenon from the PEL perspective. To establish this goal,we perform simulations of single cascades at three ambient temperatures(see Subsection 2.1 for details). As figure 2 has shown that the system has reached a steady state at 40 ps, we simulate each cascade for 42.2 ps instead of 82.2 ps used in the typical event.

We display the simulated results of single cascades at the three temperatures in Fig.4. Each data point is an averaged value over 20 runs,and the standard errors are very small(usually less than 1%) so we do not display them for simplicity(original data is available upon request). From Fig.4(a), we can see that averagely a cascade at a higher temperature creates more point defects at the peak moment but leaves less steady points finally(13,9,and 7 for 1 K,300 K,and 600 K,respectively), in line with the previous study.[44]Besides, Ndefdecreases more slowly at higher temperature,implying a slower cooling, which is understandable as the temperature gradient between the cascade region and the CTWs are smaller for a warmer environment. These two results can also be obtained from the evolution of Ndis(not shown here), although the strong thermal vibrations at 600 K introduces ～2000 background disordered atoms thus an exact comparison between 600 K and the other two temperatures is difficult.

Fig.4. Evolution of single cascades at three ambient temperatures: (a)number of the disordered atoms,(b)ground-state potential energy,and(c)groundstate displacement in the configuration space. Each line is a result averaged over 20 runs. The meanings of the grey rectangle and the grey dashed line are similar with those in Fig.3.

We plot GPE and GD of the system as a function of time in Figs. 4(b) and 3(c), respectively. It is found that the peak GPEs for the three temperatures are close to each other and the three GDs at the corresponding peak moments are also very similar (see the grey rectangle), implying that the uphill journeys in the PEL among the different temperatures do not differ obviously from each other. However,as the sample cools down, the temperature begins to affect the trajectories.More specially,GD becomes greater for higher temperature at the same moment,indicating that a longer path and thus more minimas are explored by the system. As a result, at higher temperature the system has more chances to finally locate in a lower energy minima,as shown in Fig.3(b). A notable difference between 1 K and the other two is that after 20 ps(the grey dashed line in Fig.3), the GPE and the GD at 1 K stay unchanged while those at 300 K and 600 K still slightly decrease(GPE)or increase(GD),although the Ndefvalues do not change any more(or just decrease by 1 or 2). This is because at 300 K or 600 K the system does not completely cools down at 20 ps,thus can still cross some low barriers to explore more new basins. By visually observing the cascade evolution, we find that this behavior in the PEL corresponds to the defect migration,recombination and transformation in the real crystal at the end of the cascade,[8]which is similar with the long-term defect evolution induced by thermal fluctuations[45,46]but the frequency here is much higher.

3.3. Damage buildup by cascade overlap

In this section we further expand our research from single cascade events to massively overlapping cascades,as the latter is more close to the realistic experimental situations. We plot the disordered fraction fddefined by Eq. (1) as a function of irradiation dose in Fig.5(a).Note that Ndis(and thus fd)is calculated for the ground-state configurations to remove the influence of thermal vibrations. We can see that each fdincreases linearly in the very initial stage as almost every cascade at this dose is exploring a new pristine region in the crystal. As the dose increases, the subsequent cascades begin to overlap with the defect debris left by the previous recoils,thus the increase of fdbecomes sublinear. Finally, as the sample has been homogenously irradiated, the damage begins to saturate at ～0.02 dpa(see the dashed grey line)with a nearly constant fd. We notice that the largest saturated fdin the present work is only around 3%(at 1 K),meaning that most of the sample still remains crystalline,which is quite different from the case of semiconductors like quartz[12]where the sample can be irradiated into a completely amorphous state (i.e., fdis close to 100%). This difference between metals and semiconductors has been explained by their considerably different recrystallization abilities,[3,4]which will be revisited from the PEL perspective soon later in Subsection 3.4.Note that our result of the general trend of damage evolution along with the dose(including the low saturated damage level)in hcp Zr is consistent with the previous studies on other metals,[37,38,47]suggesting that elemental metals are difficult to be amorphized solely by irradiation,which is widely acknowledged among the nuclear materials community. Another finding from Fig.5(a) is that higher ambient temperature will result a lower saturated disordered fraction, which can also be easily explained by the longer downhill journey in the PEL we have proposed in Subsection 3.2.

We depict how the GPE changes with the dose in Fig.5(b). It is observed that the evolution of GPE is quite similar with that of fdin Fig.5(a),even the local fluctuations being almost identical. This result again supports our idea that GPE is a good indicator for damage level. The final increasing amount in GPE found here for hcp Zr (on the order of 10?3eV/atom)is much lower than that for semiconductors like quartz(on the order of 0.1 eV/atom),[12]indicating a much lower damage level in the fully-irradiated metal crystal,in line with the result of disorder fraction.

Fig.5. The damage buildup in an hcp Zr crystal under prolonged 5-keV recoil irradiation: (a)disorder fraction,(b)ground-state potential energy,and(c) ground-state displacement at 1 K in the configuration space. The grey dashed line in panels(a)and(b)represents the moment at which the damage begins to saturate. The inset in panel(c)depicts the squared displacements D2 during every n cascades(averaged over 1000/n values)at 1 K as a function of n.

The evolution of GD under continuous irradiation at 1 K is shown in Fig.5(c). The changing trends of the GDs at the other two temperatures are similar with this one except for the difference in the absolute values,so we ignore them in this figure. We can see that GD monotonically increases with a trend perfectly fitted by GD=29319×dose0.5,indicating that GD2increases exactly in a linear way with dose, which is easy to understand by examining the atomic displacement distribution in a cascade. In detail, as will be shown in Subsection 3.5,among the 532224 atoms in the sample only hundreds of them experience obvious displacements (no shorter than the firstnearest-neighbor distance so the atom has left its original position)after a cascade,and only several in these atoms experience large displacements(for example,greater than 10 ?A).As a result, according to Eq.(2)the increase in GD2after a cascade mainly comprises a linear addition contributed by a few new displaced atoms, thus GD2changes linearly with dose.As every consecutive cascade is randomly initialized, we can expect that the trajectory in the configuration space under continuous irradiation will evolve like a discrete Brownian motion(Markov chain). To verify this,we plot the mean squared displacements D2for every n cascades(averaged over 1000/n values)in the inset of Fig.5(c). It is clear that D2shows a perfect linear relationship with n, in agreement with the famous Einstein relationship for Brownian motion,thus our conjecture is verified.

3.4. Comparison between metals and semiconductors

To have a better understanding on the system’s evolution in the PEL under continuous irradiation,we estimate the roughness of the local energy landscape using the method described in detail in the work of Krishnan et al.[12]For an inherent structure obtained at a certain dose,we give the system an energy pump corresponding to an equilibrium temperature of 1200 K. Then the system is relaxed under an NVE condition for 100 ps,and finally the mean squared displacement(MSD)of the atoms after this relaxation is calculated, which can be used to indicate the roughness of the landscape: The higher the MSD, the lower the energy barriers in the PEL. Our results of MSD are depicted in Fig.6. We can see that the MSD for each temperature saturates at a very low dose(before or at 0.01 dpa,see the dashed grey line),even though the defects are still increasing at this moment(see Fig.5(a)). This is possibly because the MSD of the irradiated Zr crystal after an energy pump of 1200 K is mainly contributed by the migration of single SIAs during the 100 ps,and the number of single SIAs may be already saturated at a very small dose while the subsequent irradiation only increases the number of defect clusters.[38]Then from Fig.6 we notice that the highest saturated MSD here is only around 0.25 ?A2,much smaller than the value for a fully irradiated semiconductor(around 3.2 ?A2for quartz[12]),which means that the energy barriers surrounding the ground states of the irradiated Zr crystal are much higher than those for an irradiated semiconductor. It further means that the Zr system is always trapped in rough regions in the PEL under prolonged irradiation,which is in sharp contrast with the case of semiconductors where the system can reach a very flat region in the landscape when the damage saturates.[12]As elemental metals are generally good at recrystallization and hard to be amorphized by irradiation, we can expect that the topographic features of the local PEL under irradiation presented here for Zr can also be found in other metals.

Fig.6. Atomic MSD in the 100-ps NVE relaxation after an energy pump of 1200 K in an irradiated Zr crystal.

Now we provide a qualitative understanding of the different responses of elemental metals and semiconductors to continuous irradiation from the PEL perspective. As we mentioned in Subsection 3.1,in an irradiation event,the system is driven by the incident energy to a high-elevation region in the PEL and then slowly drops into a low-altitude region. Consequently, the resulted radiation damage should be strongly influenced by the topography of the PEL, which is intrinsically determined by the atomic interactions in the many-body system. For semiconductors like quartz, the loose structure and the abundant covalent bond angles provide a vast number of energy basins(i.e.,metastable configurations)near the end of the system’s trajectory, thereby resulting in a high probability for the system to find an entrance channel to a state more defective than its initial state before the cascade. As a result, under prolonged irradiation the system will experience constantly-increasing damage until a totally glassy state is achieved where the PEL becomes quite flat.[12]On the contrary, for metals, the energy basins available at the end of the downhill journeys are quite limited because of the limited closed-packed types(fcc,bcc,hcp)and bond angles,thus the system is easy to slide down into the deep valleys corresponding to the nearly-crystalline states(with only several point defects). As a result,the sample remains mostly crystalline after a cascade,only containing some distant defects that do not recombine in the spike phase. In this way repeated irradiation to a local region in a metal crystal merely causes repeated local melting and recrystallization,and the system is always trapped in deep valleys in the PEL,thus a completely disordered configuration is difficult to be formed. Putting in other words,in a semiconductor there exits a path from the crystalline state to an amorphous state with significantly dense defective states between them, so that the system can travel along it to reach a totally damaged configuration under continuous irradiation.However, such paths usually do not exist in metals, thus an elementary metal is very hard to be amorphized by irradiation.

It is interesting to compare irradiation with homogenous quenching, by which an elementary metal can be amorphized. Theoretically, a glass of any substance can be formed as long as the cooling rate is sufficiently high to avoid crystallization.[23]However, glasses of elemental metals are quite difficult to be created because of the strong recrystallization abilities of the metals. Until just recently,monoatomic metallic glasses of tantalum and vanadium are realized in experiments as well as in simulations using an ultrahigh cooling rate of 1014K/s.[48]The irradiation events simulated in the present work show a similar cooling rate(roughly 4500/42.2 K/ps=1.06×1014K/s), but the sample remains quite crystalline with only a small fraction of point defects,seeming contrary with the case of quenching mentioned above.This contradiction can be explained, again, by the locality of the melting (energy deposition). In detail, irradiation merely introduces a nanoscale melting region which can recrystallize by following the closed-packed pattern of the surrounding crystalline atoms. On the contrary, homogeneous quenching involves a melting of the entire sample thus the atoms have no crystalline neighbors to follow,leading to a totally disordered structure when the cooling rate is very high.

3.5. Comparison between irradiation and heating

Irradiation and heating can both be regarded as energy depositions to a sample, but the former is extremely local(to only one atom)and non-equilibrium. To figure out the difference between the two ways,we simulate an instantaneous local heating in hcp Zr at 1 K with a very similar MD setup used in single cascade events. The sole difference is that the incident energy(5 keV)is not assigned to a PKA but to a sphere region at the sample center using a Maxwell velocity distribution equivalent with a certain equilibrium temperature. As a 5-keV recoil usually introduces a disordered region containing ～8500 atoms(see Fig.2(a)),we set the sphere’s radius to be 37 ?A so that it contains a similar number of atoms (actually 8973). A temperature of 4500 K is assigned to this group of atoms and the exactly deposited energy is 5.2 keV, close to that (5 keV) used in our simulations of irradiation. Like single cascades, here we also perform 20 runs(with different initial velocity distributions) to get good statistics. By visually monitoring the evolution of the sample, we find that the assigned high temperature will induce a local melting in the sample,which will be slowly cooled down by the CTWs during the simulation(still 42.2 ps for each run),just like the heat spike phase in an irradiation event. By comparing the resulted changes in the sample caused by the two deposition ways,we can have a clear picture of the difference between irradiation and heating.

The main simulated results are listed in Table 1 together with those of single cascades at 1 K.We can see that while the peak values of Ndisare similar because of our deliberate setup,all the other results are quite different. A local heating introduces obvious less damage into the sample, as indicated by the steady values of Ndis, Ndef, and ?GPE. The averaged GD for a heating event is only 42 ?A, obviously smaller than that for a cascade,implying that the trajectory in the PEL explored by the system during heating is quite different from that during irradiation, although they both involve a similar-size and similar-lifetime local melting.We examine the distributions of atomic displacements in two typical events(see Fig.7), finding that heating mainly leads to short-distance exchanges of atomic positions (see the red part in Fig.7) whereas irradiation can induce a few large atomic displacements which contribute significantly to the large GD after a cascade. These farreaching atoms usually correspond to the first one or two generations of recoils which receive considerable kinetic energies from the incident particle, so they are a natural result of the extreme locality of the energy deposition in irradiation events.As the steady damage is totally determined by the atomic configuration of the sample(defects in metals usually do not stay in electronically excited states[49]), we can conclude that the large displacements are the key difference between the effects of irradiation and heating. The possible detailed mechanism for this difference is that the greatly displaced atoms in an irradiation event can form SIAs and vacancies which are too distant from each other to recovery during the cooling process,[50]while a similar phenomenon does not occur in a heating event.

Table 1. Comparison between irradiation and heating in an hcp Zr crystal at 1 K. Deposited energies are both 5 keV. Each data point is an average of 20 simulations.

Fig.7. Atomic displacement distribution in a typical heating event and a typical irradiation event in hcp Zr. Deposited energies are both 5 keV.

We notice that previous works have obtained similar conclusions with our study. Simeone and Luneville[51]theoretically investigated the fractal nature of collision cascades and found that the atomic displacements undergo Levy flights featuring heavy tails in the distribution. This result was further experimentally verified by Boulle and Debelle[52]in irradiated SiC crystals using means of x-ray diffraction. The authors[52]further compared the atomic configurations obtained by moving the atoms according to a Levy-stable or a Gaussian displacement distribution, and found that the former was much closer to the experimental results, indicating that the large atomic displacements (corresponding to the heavy tail in the Levy-stable distribution)in irradiation events are crucial to account for the radiation damage in solids,in agreement with our conclusion. Additionally,in their pioneering work,[12]Krishnan et al. found that under irradiation the system will explore the forbidden regions of the enthalpy landscape which are inaccessible by simply cooling and heating, but the underlying mechanism for this behavior is still unclear. Here by comparing irradiation and heating under an almost identical condition,we show that irradiation differs from heating mainly in the resulted atomic displacement distribution which stems from the extreme locality of the energy deposition in irradiation events,implying that the trajectory passing through the forbidden regions upon irradiation is caused by the large atomic displacements which create distant defects and thus introduce inhomogeneity into the atomic density distribution.

4. Discussion

In this work we utilize a PEL perspective to interpret the simulated results of single cascade events and cascade overlap in hcp-Zr, while in the previous works on defect production in metals the traditional viewpoint of crystallography is usually employed.[1–4]We argue that the PEL perspective owns some unique advantages and can be a good supplementary to the traditional crystallographic methods, as illustrated by the following examples:

(i) Nuclear materials scientists usually separate the cascade process into two stages:[8,50]the ballistic phase corresponding to the strong atomic collisions and the heat spike phase corresponding to the collective motion of atoms in the local melting region. In Fig.3(a) we use the IPE to clearly show that the ballistic phase ends at ～0.2 ps in a 5-keV cascade in hcp Zr, providing a way to locate the boundary between the two stages. Besides, in the previous studies, pair distribution function (PDF) is applied to prove the existence of a local melting during the cascade process.[1,44]In Fig.3(b)we provide a stronger evidence from the PEL perspective, as liquids and solids (such as glasses) can have similar PDFs whereas only liquids can have constantly-changing inherent structures.

(ii)Reduced irradiation damage at high ambient temperature is an important result in the research of radiation effects in materials. To explain it, Gao et al.[44]used the crystallographic analysis method to find that the melt region was larger and longer for high temperature (consistent with our result presented in Fig.3(a)), indicating a more pronounced heat spike which was able to facilitate the in-cascade defect recombination.[8,50]In Fig.3(c) we provide another elucidation using the longer downhill journeys in the PEL at high temperature,which may be inspiring for future work.

(iii)There are many defect types defined from the crystallographic viewpoint,such as point defects,dislocations,grain boundaries, voids and Stacking fault tetrahedras (SFTs), but some of them are sometimes difficult to clearly identify in experiments or simulations. Besides, some defects in covalent compounds may be even ill-defined by the traditional methods,such as stretched bonds or angles,disordered atomic networks,etc. However,all these kinds of defects are equivalent from the PEL viewpoint,as every defect just corresponds one point in the multidimensional configuration space and thus can be addressed in a universal way.

It should be emphasized that we do not regard the PEL perspective to be superior over the traditional crystallographic methods, as each of them has some unique advantages that the other one does not own. For example,the crystallographic viewpoint is directly connected with the experimental observations by TEM,while the PEL one is more suitable for providing theoretical understandings and reliable predictions from the fundamental principles of condensed matter physics. We suggest that a combination of the two methods can probably achieve some intriguing results in the research of irradiation effects in materials. Actually,some researchers have used this idea to investigate the thermally activated defect evolution in irradiated materials and obtained some novel and important results.[45,46,53]We believe that this idea can also be used to explore the complex mechanisms of defect production in materials upon irradiation,thus the entire multi-scale phenomenon(including both defect production and evolution) of radiation damage in solids can be investigated in a unifying and fundamental framework.

5. Conclusions

Overall,by simulating the defect production in an hcp Zr crystal and interpreting the results from the PEL perspective,we find that an irradiation event can be regarded as a trajectory in the PEL while continuous irradiation corresponds to a discrete Brownian motion composed by these independent trajectories. Two important phenomena in irradiation experiments of materials,i.e.,high-temperature-induced damage reduction and lower damage level in metals compared with semiconductors,can both be easily explained by the dynamics of the trajectories associated with the topographic features in the PEL.Our work highlights the fundamental importance of PEL in understanding the physical origin of radiation effects in solids,and suggests a combination of the PEL perspective and the traditional crystallographic ways in future research of irradiation tolerance in materials.

Acknowledgements

The author greatly acknowledges Yue Fan for the careful reading of the manuscript with instructive comments. The author declares no competing financial interests.