A metastable phase of shocked bulk single crystal copper: an atomistic simulation study

Structural phase transformation in bulk single crystal Cu in different orientation under shock loading of different intensities has been investigated in this article. Atomistic simulations, such as, classical molecular dynamics using embedded atom method (EAM) interatomic potential and ab-initio based molecular dynamics simulations, have been carried out to demonstrate FCC-to-BCT phase transformation under shock loading of 〈100〉 oriented bulk single crystal copper. Simulated x-ray diffraction patterns have been utilized to confirm the structural phase transformation before shock-induced melting in Cu(100).

The previous studies set the stage for an investigation of the possibility of phase transition in Cu at high temperature and pressure conditions. Hirth et al. 27 postulated analytically that shock loading of Cu along the 〈110〉 direction would lead to the formation of a body-centered orthorhombic (BCO) phase. Levitas and Ravelo 28 observed loss of coordination number of the crystals behind the shock front for a shock loaded Cu in the 〈110〉 and 〈111〉 directions, which they proposed as "virtual melting". The detection of phase transformation of bulk Cu under high pressure and temperature conditions (as observed in shock experiments) remained elusive. Within the context of a polycrystal, Bolesta and Fomin 29 demonstrated nucleation of solid-solid phase transition in Cu behind a shock wave front. The effects of grain size, grain orientation of the polycrystals of Cu under shock loads of different intensities with regards to structural phase transition has also been investigated 30 . However, it should be noted that a polycrystal material is a summation of different single crystal materials of different grain sizes in different orientations having distinct grain boundaries between them. Thereby identification of body-centered phase does not properly identify which direction of single crystal of bulk Cu is more prone to phase transition. Moreover, it is also not known from the polycrystal study whether the defect induced due to grains is responsible for the phase transformation. These questions are being specifically targeted in this study with regards to the structural phase transition of Cu under shock loading.

Results
In this work multi-million atom classical non-equilibrium molecular dynamics (NEMD) simulations (based on many-body embedded atom method (EAM) interatomic potential for Cu parameterized by Mishin et al. 31 ) has been carried out for simulating piston driven shock compression of single crystal FCC copper with different orientations, 〈100〉, 〈110〉 and 〈111〉 for the piston velocity of 0.8-3.0 km/s. For the shocks propagating along 〈100〉 and 〈110〉 a structural phase transformation, FCC-to-body centered phase of single crystal copper has been identified, as characterized through adaptive-common neighbor analysis (a-CNA) 32 , radial distribution function (RDF) and x-ray diffraction (XRD) pattern analysis which has been discussed in details in this paper. Typically, the NEMD simulations are carried out for several pico-seconds timescale; thus in order to test the stability of the newly identified phase for longer time periods multi-scale shock technique (MSST) 33 simulations have also been carried out using same EAM potential of Cu 31 for smaller target samples, up to 10 ns timescale. Questions could be raised with regards to whether the observed phenomenon of structural phase transition is primarily due to the choice of the force potential [34][35][36] ; thereby, ab − initio molecular dynamics (AIMD) simulations are also performed for the same range of piston velocities. Details of the different simulation methodologies utilized have been provided in the METHODS section.
The shock Hugoniot plane (pressure and density) have been plotted for the above three simulation methods and three different orientations in Fig. 1. The results demonstrate close comparison with experimental shock data of poly-crystalline copper 37 . Figure 2 represents volume fraction of different lattice configuration at different shock intensities in different directions. Similar results have been obtained from NEMD and MSST simulation studies. Body-centered phases has been observed for shock loading along 〈100〉 and 〈110〉 directions but not along 〈111〉 direction. It should be mentioned that the a-CNA is an improvement over conventional CNA in which a fixed cutoff distance used. Conventional CNA has been replaced by variable cutoff to account for positional fluctuation as a result of high temperature or inhomogeneous atomistic strain distribution. However, it should be mentioned that even though this method gives an approximate percentage of volume fraction for different phases in the mix, it is unable to distinguish the symmetry in the crystal lattice structures (such as triclinic, monoclinic, orthorhombic, tetrahedral, etc.) apart from the cubic phase.
In Fig. 3(a) and (b), shock induced deformed matrix (obtained using NEMD simulations) is presented for Cu(100) for the piston velocities of 1.0 km/s and 1.8 km/s. Based on utilizing neighbor counting based local lattice structure identification method, adaptive-common neighbor analysis (CNA) 32 , it can be observed that for low-intensity shock compression (e.g. U p = 1.0 km/s) dislocation based plasticity (resulting in the formation of Figure 1. Shock Hugoniot point, pressure and density of the target sample with desired orientation (i.e. 〈100〉, 〈110〉 and 〈111〉) has been measured for all sets of simulations [MSST (utilizing EAM potential), NEMD (utilizing EAM potential) and ab − initio molecular dynamics (AIMD)] and compared with the available experimental shock Hugoniot data points 37 for poly-crystalline Cu. The above figures represent shock data for different shock directions, (a) 〈100〉, (b) 〈110〉 and (c) 〈111〉.
HCP lattice due to the slip assisted change in stacking sequence of atomic layers) predominates. At piston velocity of 1.0 km/s, just behind the shock front a small portion of BCC lattice (blue colored atoms as depicted in Fig. 3(a)) is observed to be nucleated due to shock-induced uniaxial strain of around 5.4%; however as soon as the shock passes and moves further the observed BCC lattices reverts back to original FCC through post-shock relaxation of shear stress and associated strain. On the other hand, for higher piston velocities (e.g. U p = 1.8 km/s resulting in shock pressure of ~105 GPa) as soon as the shock interacts with the target the entire sample transformed into body-centered phase (plastic strain along the shock propagation direction, is ~23.8%) of copper (see Fig. 3(b)) as obtained through adaptive-CNA analysis and these body centered phases are observed to survive post-shock relaxation upto 18.5 ps. RDF plots (see Fig. 3(c)) and XRD plots (see Fig. 3(d)) also demonstrate a clear signature  of body-centered tetrahedral phase. More details regarding these two analysis methodologies have been presented in the next paragraphs.
It should be mentioned that RDF plots have been used to demonstrate structural phase transition in materials 38,39 . Thus, RDF plots are also considered to analyze in this work. Figure 4 shows RDF of the atoms which are calculated from the trajectory of the shocked specimen up to 8 Å and averaged out over 5 ps time scale to obtain a clear picture regarding the temporal dependency as well as probability to obtain the 1st, 2nd and 3rd neighbors inside the specimen under the shocking situation. Results shown in Fig. 4, correspond to EAM based MSST simulation results. Similar results have also been obtained for NEMD simulations and have been found to be identical. Typically a body-centered phase has certain typical signatures such as the presence of 'shoulder region' in the descending part of the first coordination shell which indicates the presence of two merged peaks. The peaks refer to the probabilities of atom neighbors near one another. In several situations, a combination of thermal broadening of the peaks and the 'shoulder region' results in the formation of 'kink' regions (or regions with two different slopes) in the descending part of the first coordination shell. These signatures of body-centered phase (typically the presence of the 'kinks') are observed from the RDF profiles corresponding to 〈100〉 shock loading direction for piston velocities of 1.8 to 2.0 km/s. The corresponding temperature and pressure range of these piston speeds are ~1788-2286 K and ~90-105 GPa, respectively. Another typical signature of the possibility of the presence of body-centered phase is merging of second and third coordination peaks resulting in a 'wavy' nature of the RDF curve. This feature could be observed for piston velocity of 1.5 km/s (in the 〈100〉 shock loading direction) which does not show any clear peak for second and third coordination shell, rather the curve at this region of 2-4 Å is 'wavy' in nature. Similar 'wavy' nature of the second and third coordination shell peaks could also be observed for 2.5 km/s piston velocity for the 〈100〉 shock loading direction and for piston velocities of 1.8 to 2.5 km/s for the 〈110〉 shock loading direction. Apart from that, a general characteristics such as a leftward shift of first coordination shell peak, thermal broadening, and decrease in intensity of the peak as well as the gradual loss of coordination (typified by flattening nature of the higher coordination shell peaks with increase in piston velocities) could be observed for shock loading in both the directions. No specific 'kink' and/or 'wavy' nature could be observed for the shocks along 〈111〉 direction, signifying almost no probability of transition to a body-centered phase. It should also be noted that at this point, even though there are presence of certain signatures for body-centered phase from the RDF plots for the shocks along 〈100〉 and 〈110〉 directions, but these RDF profiles are not able to comment on whether the body-centered phase is cubic, tetrahedral, orthorhombic or some other lattice structures. It may also happen that due to the presence of low volume fraction of the body-centered phase the effect/signature has not been reflected in the RDF plots.
To get deeper insights, x-ray diffraction patterns have been simulated from the atomistic trajectory data. At ambient temperature and pressure conditions, XRD plots (refer Fig. 5(a)) for Cu shows primary peak corresponding to (111) plane at a interplanar spacing (d hkl ) of 2.107 nm (i.e. 2θ value of 42.92°, where lattice parameter a = 3.615 Å; volume/atom = 11.68 Å 3 , is in good agreement with ICDD copper file No. 040836). The initial structure is an FCC since h, k, l shows either all odd or all even values 40 . Lattice deformation at piston velocity of 0.5 km/s in Cu(100) is demonstrated through a slight shift in the interplanar spacing of the (111) plane to 2.014 nm (i.e. 2θ of 45.05°). Amongst the nine peaks for the piston velocity of 0.5 km/s, five depict FCC phase with a = 3.512 Å and others (peaks 3, 5, 7, 8) illustrate the HCP phase, nucleated as a result of shock-induced plastic deformation. At piston velocities of 1.5 km/s and above we have observed even values of (h + k + l) thereby demonstrating a body-centered lattice structure 40 . For 1.5, 1.8, 2.0, and 2.5 km/s piston velocity the (111) peak is no longer observed and the (200) peak splits into two parts corresponding to (200) and (002), demonstrating a typical characteristic of tetragonal lattice formation. It should be pointed out that (002) plane spacing in a tetragonal lattice structure differs from the other two planes ((200) and (020)); whereas all the three planes have same spacing for a cubic structure 40 . Consistent with the (200) peak split higher order peak splits has also been observed for MD simulations. XRD analysis of ab − initio MD simulation results (shown in Fig. 5(b)) demonstrate a rigorous proof of shock-induced structural phase transition from FCC to BCT lattice structure (characterized by splitting of the peak at (200) into (002) and (200)  The mechanism as observed in this structural phase transformation, FCC-to-BCT, in Cu(100) is primarily an uniaxial lattice compression accompanied by lattice rotation through finite volume change. As shown in the schematic diagram (see Fig. 6), a lower symmetry tetragonal body-centered unit cell (with = c a 1 2 ) can be assumed (connected via red dashed line in Fig. 6) within the two consecutive high symmetry FCC lattice unit cells. As a consequence, adequate lattice compression of 1 2 , i.e. strain of ~30% along 〈100〉 crystallographic direction transforms this low symmetry BCT unit cell to BCC lattice. However, at a piston velocity of 1.5 km/s (i.e. shock pressure and temperature of ~80 GPa and ~1130 K, respectively) spatially averaged uniaxial lattice strain is obtained as ~21.75%, which is insufficient for complete transformation to BCC and eventually a body-centered phase with tetragonal symmetry (tetragonality is ~3.3%, = . 0 967 c a ) is identified. Successive increment of the applied strain by employing higher intensity shock, e.g. piston velocity of 1.8 km/s (shock pressure and temperature of ~100 GPa and 1520 K, respectively with spatially averaged uniaxial lattice strain is obtained as ~23%), the displacement of the atoms to attain higher symmetry BCC lattice is facilitated and that is why the measured tetragonality of the obtained BCT phase is decreased to ~1%. At 2.0 km/s shock speed, further reduction of tetragonality can be obtained with the increase in spatially averaged uniaxial lattice strain to ~25%. At the much higher intensity of shock, beyond piston velocity of 2.0 km/s, amorphous solid phase along with several small crystallites (with short range ordering up to ~10 Å) has been observed to exist prior to melting.
Signatures of body-centered phases (even values of (h + k + l)) has also been observed (see Fig. 7(a)) for the shocks along 〈110〉, at a piston velocity of and above 1.0 km/s. However, the nucleated body-centered phase has been observed to co-exist with the parent FCC matrix at the piston velocities of 1.0-1.8 km/s. From the x-ray diffraction profile as obtained from the EAM based MD simulations, at the piston velocity of 2.0 km/s the (200) plane has been observed to split into (002), (020) and (200) (with the interplanar spacing of 1.817, 1.814, and 1.804 nm, respectively, i.e. 2θ = 50.23, 50.33, and 50.63°, respectively) thereby demonstrating the formation of orthorhombic lattice 40 (as earlier postulated by Hirth et al. 27 ). At 2.5 km/s, four peaks can be observed near to a 2θ angle of 50°, i.e. at the range of the interplanar spacing of 1.818-1.786 nm, which exhibits an unknown/unidentified crystal lattice structure. The ab − initio molecular dynamics (AIMD) based MSST simulations (as shown in Fig. 7(b)) however illustrate only splitting of the (220) planes into (022) and (220) plane at piston velocities 1.8 km/s and above. Even values of (h + k + l) for the x-ray peaks obtained for piston velocities of 2.0 km/s and above, indicates towards the formation of a body-centered lattice as well. Even though peak splitting is observed for the (220) planes, unlike the primary (200) plane for the shocks at the piston velocity of 2.0 km/s and above. This feature thereby evidences that a body-centered phase is being formed under shock-loading, but the exact crystal lattice structure (i.e. whether tetrahedral or orthorhombic) cannot be ascertained easily which is thus open for further investigation. In the plots (refer Fig. 7(a)) obtained from MD simulations, 'bc' and 'fc' indicates planes corresponding to body-centered and face-centered lattice respectively. The two boxes indicated by a gray colored arrow represents the zoomed view near the (002) plane for the impact velocity of 2.0 and 2.5 km/s. For 2.0 km/s In conclusion, while performing shock simulation of Cu using different shock-intensities and different orientations of the lattices, we have stumbled upon a discovery of the existence of the body-centered phase of Cu, shock compressed along 〈100〉 direction. Using three different local lattice structure identification methods, these results have been derived based on classical MD simulations considering well accepted EAM-Mishin potential for Cu. These results have also been validated through computationally expensive ab − initio MD simulations at certain piston velocities. An extensive density functional theory (DFT) study could be undertaken to determine the dip in the energy associated with the phase transition in Cu under shock loading of different intensities and in different directions which will be addressed by the authors in some later publications. Laser ablated shock compression experiments with the facilities for in-situ real-time XRD diffraction analysis, can also validate this numerical study of the shock-induced phase transformation in Cu made in this article.

Methods
Detailed molecular dynamic simulations (NEMD and MSST) 33 have been carried out in this work in which single crystal Cu in different orientations has been subjected to shock loads of varying intensities (piston velocity from 0.8 to 3.0 km/s). The force potential used for Cu in MD simulations is the many-body EAM potential developed by Mishin et al. 31 . It should be noted that EAM-Mishin potential for Cu has been parameterized based on ab − initio and tight-binding methods and is well accepted by the research community as the standard for simulations of Cu both at ambient temperature and pressure conditions and also at high pressure and temperature conditions. The potential shows good agreement with experimental elastic constants, experimental phonon dispersion curves thereby demonstrating its global reliability. Apart from the calculation of lattice properties, the potential has also been extensively tested for various structural energies and transformation paths. It should be noted that EAM potential is well referenced in many papers on shock loading of Cu [41][42][43] Tables 1 and 2 for EAM based NEMD and MSST simulations.
Three different analysis methods have been used to identify the shock equilibrated crystal structure of the samples -adaptive Common Neighbor analysis (a-CNA) 32 (lattice structure identification in multiphase systems based on coordination number of atoms), Radial distribution function (RDF -which calculate the probability of finding out two neighbors at a certain distance apart) analysis and X-ray diffraction (XRD) analysis. Identification by any one method is not reliable since it may provide a misleading prediction, as has been discussed in the article.
Computationally expensive ab − initio molecular dynamics (AIMD) simulations (using AM05 exchange-correlation function) has been carried out (simulation details and sample sizes have been presented in Table 3) and through x-ray diffraction analysis the presence of the new phase is demonstrated.
In the NEMD simulations, the initial configurations of pure single crystal with desired orientations have been created and equilibrated by applying isothermal-isobaric, NPT ensemble integration scheme (for 100 ps with timestep size of 1 fs) along with Nose-Hoover thermostat algorithm at ambient temperature and pressure condition. Periodic boundary conditions (PBC) are maintained in all three orthogonal directions during equilibration period. The details of the samples, including, crystallographic orientation, sample dimensions and the total number of atoms, has been given in the Table 1. Please note, special care has been taken while choosing the simulation domain size for ensuring elimination of any artifact of transverse periodic boundary 44, 45 over the shock-induced plasticity mechanisms. The simulation domain size for this study has been chosen in such a way that it is adequate to contain the entire shock in transverse directions and also the shocks are steady in nature in each simulation with rise time of ~0.51-1.1 ps for the range of piston velocities of 3.0-1.0 km/s respectively for all above-mentioned crystal orientations. To generate a unidirectional planar shock wave in direct-shock/NEMD simulations, so-called standard 'momentum-mirror' method is applied to the equilibrated initial configurations conjointly with the micro-canonical ensemble, NVE, to achieve the conservation of energy during shock propagation. The size of the timestep during all shock simulation is 1 fs unless otherwise mentioned. All MD simulations have been accomplished by LAMMPS 46 open-source package.
In order to carry out simulations with larger timescales, MSST algorithm (multiscale shock technique 47 which is implemented under the LAMMPS 46 framework and has been tested for simulating a unidirectional planar shock propagation for various class of materials [48][49][50][51] with sufficient accuracy) has been used to perform simulations up to 10 ns after the passage of shock. The necessary sample details have been given in the Table 2. MSST is a simulation technique based on the Navier-Stokes equations for compressible flow and follows a Lagrangian point through the shock wave which is accomplished by time evolving equations of motion for the atoms, as well as volume of the computational cell to constrain the stress in the shock propagation direction to the Rayleigh line and the energy of the system to the Hugoniot energy condition. While performing the MSST simulations, care    Table 3. Necessary details of the initial configurations of the samples with desired orientations used for ab − initio molecular dynamics (AIMD). The mass density of all the samples are 8.814 g/cm 3 at ambient temperature and pressure.
has been taken to ensure that there is no significant drift in energy at any time instance with the chosen values of 'q' and 'tscale' (parameters required in MSST simulation) for a particular shock intensity. Although, it should be noted here that any choice of 'q' (chosen as 1 for our simulations) and 'tscale' does not affect the essential generality of the background physics i.e. typically does not affect the obtained results from the simulations.
Since ab − initio calculations are computationally expensive compared to MD simulations, thereby a comparatively smaller number of atoms are taken to simulate the system. It should be noted that the essential generality of the physics remains unaltered in the simulation with a limited number of atoms as used in ab − initio molecular dynamics (AIMD). The necessary sample details have been given in the Table 3. During equilibration of the samples with desired lattice orientation (a 3D periodic box), we use a time-step of 1.0 fs and the Nose-Hoover thermostat for NVT ensemble calculations at 300 K with exchange-correlation functional of Armiento-Mattsson 52, 53 , which describes crystalline solids with better accuracy 54 than the other functional of the same class (the density and gradient based functionals local density approximation (LDA) 55 , PBE 56 , BLYP 57, 58 and RPBE 59 ). The same MSST algorithm 33 has been utilized for shock compression through AIMD simulations. All the AIMD simulations are up to 15-20 ps in timescale and an averaged energy drift from the Hugoniot energy condition has been observed as less than ~1%. Spin-resolved Armiento-Mattsson exchange-correlation function 52,53 (which includes both density gradients similar to that of GGA (Generalised gradient based approximation) based methods and surface specific estimates for both exchange and correlation energies) has been used for this simulations. It should be noted that the above-mentioned functional have shown significant fidelity for different materials such as semiconductors, simple metals, transition metals, alkalihalides and oxides 60,61 . In these simulations, a plane wave cutoff of 250 Rydberg has been estimated as sufficient to get tight convergence in stress tensor. A wave function convergence criterion of 10 −6 a.u is used for all simulations until otherwise mentioned. These AIMD simulations are accomplished by the program QUICKSTEP 62 which is implemented under the program suit of CP2K 63 open source software package.
To calculate the x-ray diffraction pattern of the shocked samples from the atomistic trajectories, the diffraction algorithm 64 (as implemented in LAMMPS) 46 has been utilized in which a 3D mesh of reciprocal lattice points is built on a rectangular grid with certain reciprocal space mesh size. The diffraction intensity at each reciprocal lattice point has been computed by performing spherical integration and corresponding normalization of intensity data by assuming all orientations are equally probably. During the calculation of diffraction intensity of all the samples shocked at different piston velocities, the 3D reciprocal grid spacing has been chosen as c1/X, c2/Y, and c3/Z; X, Y, and Z are the dimensions of the simulation box. It should be noted that since this is a simulation study of the x-ray diffraction pattern, the calculated diffraction profiles has been presented as a function of the inverse of interplanar spacing (d hkl ). The wavelength has been chosen as 1.541838 nm. The use of inverse spacing of the simulation dimensions helps to ensure that points examined within the reciprocal space mesh may be close to the point of Bragg reflection. Series of calculations with the values of c1,c2,c3 as {0.0625,0.0625,0.0625}, {0.125,0.125,0.125}, {0.5,0.5,0.5}, {1,1,1}, {2,2,2} and {3,3,3} were done, but the computed diffraction patterns didn't showed significant dependency of the mesh resolution over the diffraction pattern. All the results shown in this paper were computed by using 1, 1 and 1 as c1, c2, and c3 respectively. This {1,1,1} values were evaluated as the optimized and accurate for all shock intensities. To examine the effect of finite volume over the simulated x-ray diffraction profiles, we have performed NEMD shock simulation taking larger transverse Y and Z direction (the simulation box dimension is 180.75 × 72.3 × 72.3 nm 3 ) and have calculated diffraction pattern (see the supplementary document); no significant differences in diffraction peak characteristics has been observed. To investigate the dynamic evolution of the BCT phase over the simulated time span, x-ray diffraction pattern has also been calculated for various time instances of the atomistic trajectory and no significant differences could be observed in peak profiles (see the supplementary document), which indicates the stability of the BCT phase up to 10 ns timescale.