Evolution of elastic precursor and plastic shock wave in copper via molecular dynamics simulations

Large-scale molecular dynamics (MD) simulations are performed to investigate shock propagation in single crystal copper. It is shown that the P-V plastic Hugoniot is unique regardless of the sample's orientation, its microstructure, or its length. However, the P-V pathway to the final state is not, and depends on many factors. Specifically, it is shown that the pressure in the elastic precursor (the Hugoniot elastic limit (HEL)) decreases as the shock wave propagates in a micron-sized sample. The attenuation of the HEL in sufficiently-long samples is the main source of disagreement between previous MD simulations and experiment: while single crystal experiments showed that the plastic shock speed is orientation-independent, the simulated plastic shock speed was observed to be orientation-dependent in relatively short single-crystal samples. Such orientation dependence gradually disappears for relatively long, micrometer-sized, samples for all three low-index crystallographic directions ⟨100⟩, ⟨110⟩, and ⟨111⟩, and the plastic shock velocities for all three directions approach the one measured in experiment. The MD simulations also demonstrate the existence of subsonic plastic shock waves generated by relatively weak supporting pressures.


Introduction
Previous molecular dynamics (MD) simulations of shock propagation in copper reported the dependence of the plastic shock speed u pl s on the crystallographic direction of shock propagation within the split-shock-wave regime, which involves a fast elastic precursor followed by a slower plastic shock wave [1,2]. However, recent experiments performed using single crystal copper samples observed no such orientation dependence [3]. One possible source of disagreement between theory and experiment is drastically different length scales: in experiment mm size samples were used, whereas in simulations the longitudinal dimension of the samples was less than 100 nm. This work aims to resolve the contradiction by extending simulation samples to experimental length scales. In this regard, the important questions to answer are whether the plastic branch of the P-V shock Hugoniot is unique regardless of the orientation of the sample, its microstructure, or its length; and whether the shock (u s ) vs piston velocity (u p ) Hugoniot is dependent on these parameters. If this is the case, then it is quite possible that the orientation dependence of the plastic shock speed u pl s should disappear in MD simulations using samples with longitudinal dimensions approaching those employed in experiment. By simulating shock wave propagation in micrometer-size samples, we investigate the attenuation of the HEL with the propagation distance. In addition, a transient regime of propagation of a subsonic plastic shock wave is observed.  [4], is also used to obtain the Hugoniot data in the overdriven regime of shock propagation. Two types of piston, one with a flat, and another with a rough surface, are used, the latter being applied to facilitate the appearance of plasticity in the case of weak shocks, see Section 4. The rough piston surface has a sinusoidal shape in the lateral y and z directions, spanning two wavelengths across each direction with amplitude 1 nm. The results seem to be insensitive to the choice of the parameters describing the corrugation of the piston surface.

Computational methods
The effect of the sample's microstructure is investigated by introducing vacancies and dislocations into perfect crystal sample. In addition, a new embedded atom method (EAM) potential for copper was specifically devised by us to describe properly the extreme conditions created by the shock wave.
3. Uniqueness of the P -V plastic Hugoniot Figure 1 shows the P -V plastic branch of the Hugoniot, and the comparison between our MD results and experiment by Mitchell and Nellis [5] using polycrystalline copper samples. The experimental results are represented by the line, which is a fit of experimental data points. As is seen, all MD points lie on the experimental line, which is a good validation of our EAM potential. Importantly, the MD results include shocks initiated by flat and rough pistons, propagating in samples of different crystallographic directions, varying microstructure, and different lengths. Therefore, as in experiment, the P -V plastic Hugoniot is unique in the sense that it does not depend on the sample's orientation, its microstructure, roughness of the piston, or the propagation distance. Instead, each final state on the plastic branch of the Hugoniot is determined solely by the piston velocity. However, the P -V pathway to this final state, passing through the elastic precursor state, is affected by many factors, see details in the next section.
4. Dependence of the P -V pathway on crystallographic orientation, sample microstructure, sample length, and roughness of the piston The dependence of the P -V Hugoniot pathway to the final plastic state, on the sample's crystallographic orientation, its microstructure, and the type of piston, are displayed in figures 2, 3 and 4. The orientation dependence is shown in figure 2, where a specific case of shock waves initiated by a piston with u p ∼ 1.2 km/s, and propagating in the 100 , 110 , and 111 crystallographic directions is shown. An orientation-dependent Hugoniot pathway consists of segments of Rayleigh lines connecting initially uncompressed and intermediate elasticallycompressed states, and then the elastic and final plastic states; the latter being the same for all three directions and experiment. Split elastic-plastic shock waves are observed in the 110 and 111 directions, whereas a single plastic shock wave propagates in the 100 sample. As seen from figure 2, the pressure in the elastic precursor, which is by definition the Hugoniot elastic limit (HEL), is different for the 110 and 111 crystallographic directions: P 110 HEL = 35 GPa and P 111 HEL = 45 GPa, for the particular case of a perfect sample with longitudinal dimension L = 400 nm. As the HEL is substantially reduced upon propagation over a large distance, it becomes very small, < 1 GPa, for very large longitudinal dimensions of the samples, resulting in the absence of orientation dependence for mm-sized samples used in experiments. In contrast, the orientation dependence exists for shorter samples ∼ 100 nm used in MD simulations, which is a consequence of the crystalline anisotropy of the fcc lattice.
The effect of the sample's microstructure on the P -V Hugoniot pathway to the final plastic state is shown in figure 3, where results for a perfect Cu 110 sample are compared to those obtained for samples containing vacancies, and pre-existing dislocations, with all three cases corresponding to the same piston velocity u p = 1.6 km/s. For all three samples, the splitshock-wave regime is observed, resulting in the same final plastic state. However, P HEL is drastically reduced: by about one third in the sample with vacancies at the intial concentration 10 −3 , and by one half in the sample with pre-existing dislocations at an initial concentration around 5 × 10 15 m −2 . The reduction of the HEL upon the introduction of point and extended defects brings the results of MD simulations closer to experiment, which is the right trend as the experimental samples are never perfect crystals, containing a substantial amount of vacancies, dislocations, stacking faults, or grain boundaries.
The effect of piston roughness is shown in figure 4. The MD simulations demonstrate that the piston quality does not affect the physical phenomena generated by strong shock waves. However, the piston roughness facilitates the onset of plasticity in the case of weak shock waves. In particular, for a flat piston moving with u p = 0.47 km/s in Cu 110 , only an elastic shock wave is observed. In contrast, a rough piston induces a substantial shear stress, resulting in dislocation nucleation and the development of plasticity in the vicinity of the piston [6]. In this split-shock-wave regime, the speed of the plastic wave, as reflected by the slope of the Rayleigh line, is very small. In fact, it is a subsonic shock wave propagating in the pre-compressed material with a velocity less than the local sound speed. Nonetheless, it is still a shock wave as it exhibits jumps in temperature, pressure, and density at the plastic front; see figure. 5. This unusual regime of shock-wave propagation is discussed in section 5. The important result obtained in our large-scale MD simulations is the decay of the pressure in the elastic precursor, or HEL, with propagation distance. figure 6 displays the evolution of the pressure in the elastic front (P el ≡ P HEL ) and the the jump in pressure at the plastic front ∆P pl = P pl − P el , (bottom panel), together with the elastic and plastic shock wave speeds as a function of the propagation distance for the case of a Cu 110 sample containing pre-existing dislocations and u p = 0.92 km/s. P HEL is reduced from 35 GPa to 10 GPa after the shock wave propagated about 1.5 µm. The attenuation of the HEL is due to the growth and multiplication of pre-existing dislocations, and the elastic-plastic collapse of the overcompressed elastic state [7]. The decrease of P HEL results in the increase of the jump in pressure at the plastic front ∆P pl , causing the acceleration of the plastic wave. At the end of the MD simulation, the plastic wave speed approaches the experimental value corresponding to the same piston velocity [3]. Concurrently, the elastic wave slows down, its speed approaching that of the plastic wave; see figure 6. Although the attenuation of the HEL with propagation time in metals is well-known among experimentalists [8][9][10], it is often ignored by MD simulators as short samples ≤ 100 nm do not allow observation of the HEL decay.
5. u s -u p Hugoniot and subsonic plastic shock waves Due to the dependence of the HEL on the propagation distance, the calculated u s -u p plastic Hugoniot reports shock velocities corresponding to the specific length of the sample used in  simulations. The MD data would change if samples of different lengths were used. In contrast, the experimental data by Mitchell et al. were obtained using mm-sized samples. Such large propagation distance yields extremely weak elastic precursor, the corresponding P HEL being very small compared to the pressure in the plastic wave. Therefore, some disagreement between theoretical predictions and experiment is expected. In fact, the MD results agree well with experiment for strong shock waves; but for relatively weak shock waves, the plastic shock speeds obtained from MD simulations are substantially smaller than corresponding experimental values; see figure. 7. It is expected that the MD points would approach those of experiment if the simulated sample length was extended to experimental scale.  Figure 7. u s -u p plastic Hugoniot for copper. Comparison of MD results (this work) with single crystal (Chau [3]), and polycrystalline (Mitchell [5]) experiments. Points corresponding to transient subsonic plastic shock waves lie below the horizontal line corresponding to the bulk sound speed in uncompressed solid.
As it was mentioned above, subsonic plastic shock waves are observed in our MD simulations; see figure 7. Although they might be seen in previous MD simulations [2], and experiment [11], such unusual behavior has not yet been analyzed in detail. Conventional wisdom implies that a jump in thermodynamic quantities corresponding to the shock wave front always propagates with a supersonic speed; however, such an assertion assumes steady-state flow.
Therefore, it is quite plausible that a subsonic regime might occur during the initial stage of shockwave propagation, which is characterized by the slow development of plasticity.
As the pressure in the elastic precursor drops substantially over the course of the shock wave propagation, the velocity of the plastic front, defined by the pressure jump at the plastic front, constantly increases, and at some point the plastic wave becomes supersonic, its velocity exceeding the sound speed in the elastically pre-compressed material. Further details are provided in the forthcoming publication [12].

Conclusions
Our large-scale MD simulations on shock-compressed copper unambiguously demonstrate that the P -V plastic Hugoniot is a unique function of the piston velocity, in accordance with experimental observations. The P -V pathway, however, is shown to depend on the crystallographic orientation, microstructure, and the length of the sample. The MD simulations demonstrate the decay of the pressure in the elastic precursor with propagation distance, in good agreement with experimental observations. This phenomenon explains the appearance of the unusual regime of subsonic plastic shock wave propagation.