Atomistic material behavior at extreme pressures

Computer simulations are routinely performed to model the response of materials to extreme environments, such as neutron (or ion) irradiation. The latter involves high-energy collisions from which a recoiling atom creates a so-called atomic displacement cascade. These cascades involve coordinated motion of atoms in the form of supersonic shockwaves. These shockwaves are characterized by local atomic pressures >15 GPa and interatomic distances <2 A. Similar pressures and interatomic distances are observed in other extreme environment, including short-pulse laser ablation, high-impact ballistic collisions and diamond anvil cells. Displacement cascade simulations using four different force fields, with initial kinetic energies ranging from 1 to 40 keV, show that there is a direct relationship between these high-pressure states and stable defect production. An important shortcoming in the modeling of interatomic interactions at these short distances, which in turn determines final defect production, is brought to light.


INTRODUCTION
As scientific computing becomes evermore ubiquitous, it is common practice to simulate the effect of extreme environments on materials using molecular dynamics (MD). The ability to capture a material's response atom by atom has helped understand and complement experiments. For example, thanks to MD, the thermodynamic processes that control the production of nanoparticules after pulsed laser irradiation are relatively well understood. 1 Notably, such computational results explain the Newton rings observed in experiments. 2,3 MD also led to a better understanding of the nature of the destructive shockwaves that follow high-energy impact [4][5][6] and the behavior of materials at high pressure in diamond anvil cells. [7][8][9] It can predict many properties, such as the Hugoniot, dislocation densities after impact and fracture behavior. In the case of neutron and ion irradiation, MD can be used to predict the evolution of ballistic collisions that occur, and the nature of primary damage produced by the high-energy recoils. 10 It can also shed light on a plethora of mechanisms that take place as the kinetic energy of the recoil is dissipated: fractal replacement sequences, 11 supersonic shockwaves, 12 sonic waves, 12 the formation of liquid-like regions and their recrystallization or amorphization. 13 The aim of many MD studies is to describe the mechanisms and general trends, not to make precise predictions. However, as the scientific community moves towards materials by design and systematic computational characterization, increasingly precise and accurate models are required. In the case of atomistic modeling, methods that explicitly account for electronic structure, such as ab initio MD, are considered the gold standard. Unfortunately, their use has very serious limitations, given their high computational cost. In particular, the simulation of high-energy perturbations such as those involved in neutron or ion irradiation requires large simulation cells: a few hundred thousands of atoms at a minimum, but sometimes reaching billions of atoms. Most electronic structure codes do not scale well past a few hundred atoms, especially if modeling metals.
In addition, the high initial atomic velocities involved in these simulations limit timesteps to a fraction of those used in equilibrium MD (~10 − 18 s rather than~10 − 15 s). For these reasons, MD based on classical (or semi-empirical) force fields is currently the only technique used to simulate these energetic events.
These classical force fields are generally parameterized to obtain agreement with equilibrium and near-equilibrium properties, such as elastic constants, phonon spectra, phase transitions, and point-defect formation and migration energies. Given the high velocities involved, atoms move much closer to each other than under equilibrium conditions. Figure 1a illustrates this point for iron. The shortest atom-to-atom distance in a displacement cascade is plotted as a function of the energy of the primary knock-on atom (PKA). Distances approaching 0.7 Å, i.e., 0.24 of the equilibrium lattice parameter, are observed. To adequately model the interaction energies at these short distances, one cannot extrapolate the force field based on equilibrium properties. For very-short-distance interactions (less than~1.2 Å), it is common practice to use the 'universal' screened Coulomb force fields developed by Ziegler, Biersack and Littmark, which were calibrated using Hartree-Fock calculations of two atoms in vacuum. 14 The potential energy and forces for distances between 1.2 Å and near-equilibrium distances of~2 Å are obtained from an arbitrary interpolation of the equilibrium and Ziegler, Biersack and Littmark force fields. The potential energies corresponding to these short distances are given for the case of Fe in the secondary y axis of Figure 1a, where values as high as 210 eV per atom are seen for the case of a 0.9-Å separation.
The transition from the supersonic shockwave regime to the sonic regime is of particular interest. Indeed, the supersonic shockwave is considered destructive, whereas the sonic is considered non-destructive. 12 The former involves high kinetic energies at the shockwave front, larger than the displacement energy threshold. In other words, these atoms have enough momentum to potentially permanently dislodge their neighbors 1 from their original lattice sites. The later involves primarily small elastic displacements, which eventually lead the atoms back to their original lattice position. Characterizing the system at the point in time of this transition should provide insights into the nature of primary damage formation.
In this report, the results of simulations involving energetic knock-ons in Ni are presented and analyzed. The appearance of coordinated atom motion, leading to a supersonic shockwave, is presented and characterized. It is shown that the displacements created by this shockwave are highly correlated with the number of final stable defects created by the PKA in the Ni crystal. At this stage of cascade development, the local material state is described by the least accurate part of classical force fields. Furthermore, the crucial importance of correctly modeling short-range interatomic interactions is highlighted.

RESULTS
The evolution of high-energy displacement cascades can conveniently be subdivided into three distinct phases: a supersonic phase, a sonic phase and a thermally enhanced recovery phase. 12 A typical time profile of the number of atomic displacements is plotted in Figure 1b. The three phases are indicated by different colors. Note that the transition from the supersonic regime to the sonic regime is associated with an inflection point in Figure 1b. At this moment in time, the fastest atom in the system typically has a kinetic energy lower than that of the 'displacement threshold' of~40 eV. As a result, beyond this time, atoms no longer have sufficient energy to permanently displace a neighbor from its lattice site. Therefore, the increasing number of atoms displaced during the sonic regime are short-lived transient defects, which will recover their original lattice position shortly thereafter. This observation is consistent with those of Calder et al. 12 A noticeable feature of the supersonic regime is that the PKA creates many secondary and higher generation knock-ons. The coordinated motion of these energetic atoms creates a supersonic pressure wave. This is illustrated in Figure 2. The pressure waves are easily identifiable using the local atomic stress or atomic volume. Three moments in time are illustrated: at 0.10, 0.15 and 0.25 ps after cascade initiation. These correspond to times before, during and after the supersonic to sonic transition. The pressure wavefront is characterized by small atomic volumes and high local pressure (i.e., black spots in Figure 2). The supersonic wave leaves behind a region of high atomic volume and small or negative local pressure. In this cross-section, one can identify two interacting shockwaves. The low-density pockets are well known to be precursors of vacancy-type defects formed at the end of cascade evolution, and the interaction of nearby shockwaves is known to produce large interstitial-type clusters. 12 Likewise, the interactions of these shockwaves with a free surface can lead to the surface roughening and creation of huge vacancy dislocation loops. 15,16 In the 0.25-ps snapshot, one can see that the sonic wave propagates without creating low-density regions. For this reason, the supersonic regime is often considered destructive, whereas the sonic regime is considered non-destructive. Finally, during the thermally enhanced recovery regime, the highly disordered region created by the supersonic wave cools down and recrystallizes. The atoms displaced by the sonic wave recover their original positions. During this phase, temperatures can be very high in the core of the cascade (reaching over 4,500 K in the low-density pockets of our 10 keV simulations); atoms that were displaced during the supersonic phase remain very mobile and many recombine with vacant lattice sites. By~10-20 ps, thermal equilibrium is restored and a number of stable atomic displacements (vacant lattice sites and atoms in interstitial positions) remain. A complete sequence of snapshots illustrated in Figure 2 can be found in Supplementary Figure S1 in the Supplementary Materials, along with Supplementary Movie S2-S4.
In Figure 3, the correlation between the number of atoms displaced at the end of the supersonic stage and the number of stable defects at the end of the cascade is demonstrated. The dashed line is the best linear fit to the data. Each data point is averaged over 16 runs, using a given initial PKA energy (1, 5, 10, 20 or 40 keV) and force field (data obtained using two Ni force fields and a Ni-Co and Ni-Fe force field). The trend is very clear. It is possible to predict, on average, how many stable defects will be present in the system by measuring the number of atoms displaced by the supersonic wave. This shows the crucial importance of correctly modeling these supersonic waves. One should note that this is not only a linear log-log relationship, but a true proportionality. In other words, the line in Figure 3 is a power law with an exponent of 1.
The magnitude of the stress involved in these waves is remarkable. As shown in Figure 2, some atoms exhibit local pressures as high as 50 GPa in 10 keV cascades. Such pressures are comparable to those observed during short-pulse laser ablation, high-energy impacts and diamond anvil cells. The interatomic distances involved are equally remarkable features of this process. As shown in the Supplementary Materials, they are typically

DISCUSSION
The simulations herein exhibit extreme conditions. They raise serious questions about the accuracy of results obtained with the force fields used by the scientific community to model such highenergy phenomena. Indeed, as explained in the introduction, the parameterization of these potentials in the critical between 1.2 and 2.0 Å is very loosely constrained in practice. However, these are precisely the distances involved in these high-pressure processes. In the case of displacement cascades, Figure 3 shows that configurations involving short distances are directly related with final defect production. Although current modeling techniques are probably adequate to identify general features of these high-energy phenomena, our analysis indicates that the force fields used to compute short-range interactions require more physically based constraints to generate quantitatively accurate simulation results. It is also important to note that the two Ni potentials showcased in this study 17,18 have very similar equilibrium properties, but nonetheless exhibit different number of stable defects. Likewise, both Ni potentials and their respective alloys (Mishin's Ni-NiCo 18,19 and Bonny's Ni-NiFe 17 ) have largely similar interstitial formation energies and displacement energy thresholds. Using the Mishin potentials, the o1004 dumbbell interstitial formation energy is 3.97 eV in Ni and has a range of 3.8-4.2 eV in NiCo. For the Bonny potentials, the o1004 dumbbell interstitial formation energy is 5.83 eV in Ni and has a range of 4.5-7.5 eV in NiFe. For displacement thresholds, the respective values are 37, 40, 55 and 51 eV. The divergence in these values is not enough to explain the divergence in the number of stable defects produced by each potential. This suggests that the variations in stable defect production from one potential to another are largely due to short-range interactions between 1.2 and 2.0 Å.
Formally, the effect of these short-term interactions on the final outcome was only shown for the case of neutron and ion irradiation. However, the pressures involved are comparable to those of other high-energy processes, such as laser ablation, high-energy impact and diamond anvil cells. [1][2][3][4][5][6][7][8][9] This suggests that further work is needed to improve the description of atomic force fields at short interatomic distances in order for future computational work using these force fields to be as quantitatively accurate as possible. Some efforts have previously been made in this direction. For example, Tiwary et al. 20 took special precautions to ensure that their uranium oxide potential was based on the ab initio data at these distances. Zhakhovskii et al. 21 showed that high-tension states should also be included in the fitting database of an embedded atom method (EAM) potential to properly describe short-pulse laser ablation. Likewise, Zhang el al. 22 demonstrated that properly reproducing generalized stacking fault energies is of crucial importance to model compressive shocks. Our study strongly indicates that such efforts are justified and should be pursued in a systematic fashion.

MATERIALS AND METHODS
Displacement cascades with 10-40 keV of initial kinetic energy were simulated; this energy range is relevant to the knock-ons produced during neutron or MeV ion-beam irradiation. MD simulations were performed Figure 3. Near-linear relationship between the number of atoms displaced at the end of the supersonic stage and the final number of stable defects (R 2 = 0.93). It is important to note that the power law in this log-log plot has an exponent of 1. Data obtained from cascades with 1, 5, 10, 20 and 40 keV of initial kinetic energy; different colors indicate different classical potentials. Figure 2. Snapshots of atomic volume and local atomic pressure during a 10-keV displacement cascades in Ni, at times 0.1, 0.15 and 0.25 ps after the creation of a PKA. A 5-atom-thick cross-section of a 256,000-atom system is shown. These times correspond to the transition between the supersonic and sonic phase.
using the LAMMPS package 23 in an atomic system thermally equilibrated at 300 K. The PKA was given initial velocity in the [135] direction. A variable timestep was used to maintain accurate integration of the equations of atomic motion. Interatomic interactions were modeled using four force fields: two for Ni, 17,18 NiFe 17 and NiCo. 18,19 The number of displaced atoms was determined using a Wigner-Seitz cell analysis. Local pressure was computed using the trace of the atomic virial stress tensor and dividing by the atoms's Voronoi volume. The Voronoi volume calculation and vizualisation are performed with the Ovito package. 24