Quasi-isentropic ramp loading experiments can achieve Mbar-scale uniaxial compression while keeping samples at relatively low temperatures below the melting curve. Over the decades, it has become a very important technique in addition to shock compression experiments for the study of materials, especially solids under high compression.1–9 In a typical ramp loading experiment, a smoothly increasing pressure source is applied on one side of a multistep planar target, launching compression waves. Once the compression waves reach the rear side of each step, the free surface begins to accelerate, and the velocity profiles are measured by interferometry. The stress–density relation at the loading side is then inferred from these free surface data by performing inverse calculations of the flow field. Such data analysis methods were developed decades ago,10 and among the most widely used approaches are backward integration11,12 and characteristics-based Lagrangian analysis.13–16 Although these back-calculation methods have different mathematical implementations, the underlying physical basis is essentially the same, namely, solution of the inverse problem of one-dimensional isentropic hydrodynamics. These methods have been validated in many previous studies using data generated by forward hydrodynamic simulations,17,18 demonstrating that solving the inverse problem alone does not bring any technical difficulty by itself.

However, the problem remains that such a hydrodynamic approach of inverse calculation assumes the longitudinal stress P to be related to the density ρ by a well-defined function P(ρ) throughout the entire flow field. Although the function P(ρ) itself may not necessarily be the true isentrope, such a one-to-one correspondence between P and ρ is required to satisfy the formulation of the “isentropic” hydrodynamic equations19∂∂tρ0ρ=∂u∂h,ρ0∂u∂t=−∂∂hP(ρ),(1)where ρ₀ is the density before compression, u is the fluid velocity, h stands for the Lagrangian coordinate, and P is the longitudinal (not hydrostatic) stress.

However, for real-world materials, especially solids, the stress–strain behaviors may be strongly affected by structural phase transitions and/or elastic–plastic response. This makes the stress depend not solely on the density at a given moment,20 and thus in principle leads to violation of the basic assumption of these hydrodynamic-based methods for back-calculation.

Whether such lack of rigor will have serious consequences on the overall validity of the back-calculation process remains a problem. On the one hand, this technical issue is often overlooked in magnetically and laser-driven ramp loading studies. On the other hand, it has also been proposed from the perspective of continuum elastic–plastic modeling21,22 that the results of iterative Lagrangian analysis may be “compromised by local elastic release” in some situations, depending on the specific ramp loading waveform.

To investigate this problem in detail, we perform a large-scale molecular dynamics (MD) simulation of a typical ramp loading process. MD contains detailed physics of solids, not relying on any approximation of the constitutive relation. Microscopic details concerning ramp wave propagation have previously been approached using MD simulations,23–26 with most attention having been paid to analyzing the state of matter achieved in the ramp loading sample under MD simulations. However, in this study, we would like to move one step further by focusing on the back-calculation process. The free-surface velocity profiles generated by MD are used for back-calculation of the stress and density values.

The results turn out to exhibit significant deviations from what is directly measured from the MD simulations. By presenting a careful comparison between the one-dimensional flow field given by MD and hydrodynamic modeling, it is found that this kind of deviation can ultimately be attributed to the inherent irreversible behaviors of solids that are closely associated with structural phase transitions and elastic–plastic response—an important consequence of these effects being that rarefaction and compression waves may propagate at different sound velocities after entering the irreversible region, leading to failure of the isentropic hydrodynamic equations.

An embedded atom method (EAM) potential for high-pressure iron27 is used here. This interatomic potential has been carefully calibrated against experiments and density functional theory calculations,27 and so the insights drawn from its results should be relevant to real-world situations. Iron modeled by this potential displays considerable strength and a structural phase transition [from body-centered cubic (bcc) to a mixed structure of hexagonal close-packed (hcp) and face-centered cubic (fcc)] at roughly 15–25 GPa under quasi-isentropic compression, which provide the nontrivial stress–strain features in which we are interested.

Planar samples with four different thicknesses subjected to the same ramp drive are used to model the multistep target. For each of the thicknesses, we use a rectangular simulation box with periodic boundary conditions in the x and y directions. Compression is applied in the z direction by setting the lower z boundary as a reflective piston that accelerates following a designed path. The upper z boundary is a free surface [Fig. 1(a)]. The total simulation time is 100 ps, with step size of 1 fs. Simulation trajectories are recorded every 0.2 ps. At each sampling moment, the sample is divided into 400 equally spaced chunks along the z axis, and quantities such as density, particle velocity, virial stress tensor, and thermal pressure are computed and saved for each chunk. The mean particle velocity in the z direction of the last chunk is taken as the free-surface velocity.

Three different ramp loading cases are simulated in this study: liquid iron, single-crystalline iron, and nanocrystalline iron with random grain orientations. The sizes of the simulation boxes for each case are listed in Table I.

The simulation samples are carefully generated and equilibrated before ramp compression is applied. For the first two cases, the samples are generated in the bcc structure and then relaxed using a constant number of particles, pressure, and temperature (NPT) ensemble at zero pressure and the desired temperatures (1600 K for liquid and 300 K for solid). The residual stress after relaxation is of the order of 0.01 GPa. For the final case, the nanocrystalline samples are generated by Voronoi tessellation using the Atomsk code,28 and the same NPT relaxation is applied. All MD simulations in this study are performed with the LAMMPS code.29

We have also carried out forward hydrodynamic modeling for the ramp loading cases, against which the flow field extracted from MD simulations are compared. For this model, the isentropic hydrodynamic Eq. (1) are solved numerically by the forward Euler method, with the time step set small enough to guarantee convergence. These solutions of the isentropic hydrodynamic equations represent the primary model of ramp loading.

For each case, the P(ρ) relations used in solving the hydrodynamic equations are obtained by independent MD simulations of “wave-free” compression, in which a smaller sample with periodic boundary conditions in all three directions is deformed uniaxially and homogeneously at a fixed strain rate. The time scale for sound waves to transit through the simulation box is much faster than that of deformation, and so the compression is effectively homogeneous, free of “waves.” The MD samples for wave-free compression were prepared in initial conditions identical to those of the corresponding ramp-loading cases—being relaxed at the same temperature and density, and having the same grain size distribution for the final nanocrystalline case. The longitudinal stress vs density curves P(ρ) obtained from these wave-free simulations are also used as a reference for the subsequent results.

As one of the most widely used back-calculation methods, iterative Lagrangian analysis15 is implemented here. We do not repeat the data analysis with other methods, because they are essentially equivalent, all being based on the inverse calculation of isentropic flow field. The workflow of the iterative Lagrangian analysis is as follows:1.Starting from the present assumption of the P(ρ) relation, positive and negative characteristic lines are calculated backwardly from the free surface using the measured velocity profile. This gives the present reconstruction of the flow field, from which we find the corrected arrival times of each compression wave at the free surface as if they were not affected by rarefaction waves.2.The velocities of the compression waves can thus be fitted using step thicknesses and the corrected arrival times. This gives the Lagrangian sound velocity as a function of the particle velocity, c_L(u_p).3.The characteristic lines can be recalculated using these updated sound velocity data. Steps 1 and 2 are repeated until convergence.4.After convergence, the stress–strain relation is obtained by integrating the c_L(u_p) relation:ρ=ρ01−∫0udu′cL(u′)−1,Pzz=∫0uρ0cL(u′)du′.

Finally, we compare the back-calculation result against the real stress and density values measured directly from the MD trajectories.

We first present an example of liquid iron to demonstrate that our implementations of MD simulation and back-calculation work well when the nontrivial solid-dominated irreversible features are not present.

The piston acceleration history for ramp loading follows the standard design strategy that tries to delay shock formation as much as possible for a fixed ramp time.30 The reference P(ρ) relation used in the design is obtained by a preliminary wave-free simulation with strain rate 10¹⁰ s⁻¹. Here, the piston is designed to accelerate for 60 ps, from 0 to 2.5 km/s, generating ∼100 GPa pressure according to the reference P(ρ) relation. All compression waves are designed to converge at the same point, 196.4 nm away from the initial piston position. Hence, no shock is expected to form if the sample thickness is less than this value.

With this design, the typical one-dimensional Lagrangian flow field would be as in Fig. 1(b). No shock is generated. Shown here are the contours of particle/fluid velocity u_p extracted from MD simulation of a 108 nm thick sample (red). As expected, they are in good agreement with the hydrodynamic model for this liquid sample (blue), except for the very first compression waves at extremely low pressure.

MD simulations of this kind are carried out for a total of four different thicknesses ranging from 108 to 144 nm, representing a target with four steps. Figure 1(c) shows the free-surface velocity profiles in these cases, from which the c_L(u_p) relation [Fig. 1(d)] and then the P(ρ) relation [Fig. 1(e)] are back-calculated by self-consistent Lagrangian analysis.

At the same time, the real density and longitudinal stress in the MD simulations are directly measured for all four cases, and for all temporal points before the maximum stress (i.e., excluding the final release phase, of which the free-surface velocity data are independent). The spatial sampling points cover the entire simulation sample. In Figs. 1(e) and 1(f), the shaded bands show the mean, standard deviation, and minimum/maximum values of the statistics. Clearly, the back-calculation results have satisfactorily reproduced the real stress–density trajectory. The remaining error is only about 0.5 GPa, comparable to the inherent fluctuation of the recorded MD trajectories. Therefore, for the case of liquid iron, MD simulation, free-surface velocity extraction, and Lagrangian analysis all work well, as expected.

Next, we move on to ramp compression for solids. A properly designed piston acceleration history that prevents shock formation is needed in the first place. However, such a design will not be as straightforward as in the previous case, owing to the unique properties of the stress–density relation here in solids. So, we first discuss briefly how the design is carried out.

Figure 2(a) shows the relation between longitudinal stress and density obtained by the preliminary wave-free compression at a strain rate ε̇=2×1010 s⁻¹. The sample starts as a bcc solid at 300 K. Compression is applied along the [111] crystalline direction, which exhibits much smaller rate effects compared with the [100] and [110] directions. As shown in Fig. 2(a), deformation remains elastic below 15 GPa; above 15 GPa, there is a region dominated by plastic deformation and structural phase transition, where the density increases without a significant change in pressure. The Lagrangian sound velocity c_L thus shows a dip between 8.5 and 9.5 g/cm³, which would cause splitting of elastic and plastic compression wave fronts as in Ref. 31. After smoothing, it is used for the subsequent designs.

We perform the design of the piston accelerating history by the method of characteristic lines. With the standard design strategy, i.e., requiring all compression waves after the c_L minimum to converge at the same point, the resulting characteristic lines in thick samples would be the black lines in Fig. 2(c). However, if the reflected rarefaction waves generated at the free surface are considered, because of the nonmonotonic behavior of the sound velocity as shown in Fig. 2(b), the rarefaction waves prior to the c_L minimum would distort the upcoming compression waves significantly, making them intersect each other much earlier than expected [Fig. 2(c), red curves]. In other words, a design that delays shock formation as much as possible in the absence of rarefaction waves will actually lead to early formation of a shock if the reflected rarefaction waves are taken into consideration.

To achieve shock-free ramp compression for this case with nonmonotonic c_L(ρ), we perform Monte Carlo optimization of the piston acceleration history. The piston is required to accelerate from zero to 1.3 km/s within 60 ps, which is expected to generate about 70 GPa pressure. The specific way in which the piston accelerates is adjusted by the Metropolis algorithm, and the characteristic lines are calculated forwardly using the reference c_L(ρ) relation [Fig. 2(b)]. The goal is to minimize self-intersection among the compression or rarefaction waves while avoiding reverberation, i.e., the rarefaction wave should not return to the piston surface before the acceleration phase finishes. Here, two cases with thicknesses 180 and 240 nm are considered simultaneously in this optimization, because ramp loading needs several different step thicknesses. The optimized result is satisfactory: self-intersections of characteristic lines have been effectively removed [Fig. 2(d)].

MD simulations of shock-free ramp loading can then be carried out using such a well-designed piston trajectory. Here, ramp loading is again applied in the [111] crystalline direction. Four different thicknesses of 176, 192, 208, and 224 nm are simulated, and the free-surface velocity profiles are shown in Fig. 3(a). The c_L(u_p) relation obtained by Lagrangian analysis using these data is given in Fig. 3(b). Notably, the back-calculated sound velocity shows significant deviations from the reference curve for u_p > 0.4 km/s, which corresponds to the onset of plasticity and phase transition.

These deviations propagate to the back-calculated stress–density curve. In Fig. 3(g), the longitudinal stress and density values measured directly from the MD trajectory are presented along with the back-calculation result. The latter is nearly 20% higher than the former. Note that such a deviation cannot be explained by the fluctuations in the MD simulation itself, since the it greatly exceeds the min/max range of the fluctuations.

Although different parts of the ramp-compressed sample experience different strain rates and in principle may exhibit different stress–density relations, such rate-dependent effects are actually small in the simulations here—owing to the carefully designed piston acceleration history that has prevented shock-like accumulation of compression waves. Here, the min/max range of the stress–density relation across the entire flow field is ∼±2 GPa, while the back-calculation error is around +8 GPa. Therefore, to a precision of ±2 GPa, a well-defined P(ρ) relation still exists inside the ramp-compressed sample, but Lagrangian analysis failed to reveal it.

To identify the factors that may be related to this kind of deviation, we present two additional analyses as negative controls. First, nonmonotonicity of the sound velocity alone should not affect the validity of Lagrangian analysis in the context of isentropic hydrodynamics. To prove this, we generate simulation data using the hydrodynamic model with the reference stress–density relation of ε̇=1010 s⁻¹, which is strictly a one-to-one correspondence between P and ρ without any kind of irreversibility. For these u_fs data [Fig. 3(c)], Lagrangian analysis yields exactly the correct c_L(u_p) relation [Fig. 3(d)]. Therefore, in the framework of isentropic hydrodynamics, nonmonotonic c_L does not give rise to any problem.

The second possibility we would like to exclude concerns MD simulations: can the speed of infinitesimal compression waves in MD simulations be different from that obtained by dP/dρ just by itself? We note that sound velocities in MD simulations can be directly measured for simple compression waves [Fig. 3(e)]. Before being affected by rarefaction waves, the simple compression waves are straight lines, the linear-fitting slopes of which are exactly the corresponding Lagrangian wave velocities. This fitting can be performed for each contour line of u_p. Thus, the true Lagrangian sound velocity c_L as a function of particle velocity u_p for simple compression waves is measured directly. This again matches well with the reference c_L(u_p) curve [Fig. 3(f)] and gives the correct stress–density relation after integration [Fig. 3(g), the magenta line]. So, before the rarefaction waves are encountered, MD simulation still matches the hydrodynamic model well.

Therefore, these deviations must be related, first, to the difference between MD simulation and the hydrodynamic model, and, second, to the wave interaction region, i.e., the region where flow field is reversely constructed using the characteristic line formulation of the hydrodynamic equations.

Since the deviations are related to the difference between MD and hydrodynamic model in the wave interaction region, we present detailed comparisons between them in the sense of flow field. The 208 nm case is taken as an example. In Figs. 4(a) and 4(b), the time derivatives of the longitudinal stress at each Lagrangian position, ∂P/∂t, are shown as heat maps for visualization of the compression and rarefaction waves. Overlapped are the first rarefaction wave and the last compression wave given by Lagrangian analysis. For the case of hydrodynamic model, the back-calculated characteristic lines match well the underlying flow field features [Fig. 4(a)].

However, for MD, the back-calculated trajectory of the first rarefaction wave shows visible deviations after u_p ≳ 0.4 km/s [Fig. 4(b)]. Note that u_p ∼ 0.4 km/s corresponds exactly to where the deviation in the back-calculated c_L(u_p) relation begins to appear in Fig. 3(b). The first rarefaction wave in the MD simulation seems to be faster than that expected from back-calculation. It is also faster than that in the hydrodynamic model, which is evident on comparing the positions marked by the white arrowheads in Figs. 4(a) and 4(b).

Further examination of the MD flow field indicates that the “abnormally” fast rarefaction wave goes beyond the scope of the isentropic hydrodynamic model. For the one-dimensional flow field governed by Eq. (1), the sound velocity is fully determined by the local value of the density by dP(ρ)/dρ. Both the infinitesimal compression and rarefaction waves propagate at this sound velocity when passing through this spatial–temporal point. This is the immediate consequence of a well-defined one-to-one P(ρ) relation. However, this basic property is not satisfied, as shown in Fig. 4(c). The local sound velocities of the first rarefaction wave are characterized by its tangents, while the simple compression waves next to it also have the local sound velocities as their slopes. By mirror-imaging graphically the linear fittings of the latter, it is clear that these two sets of slopes do not match each other, especially for u_p between 0.4 and 0.75 km/s. That is to say, at least in this zone of the MD flow field, rarefaction waves and their sibling compression waves travel at different sound velocities. This violates the basic assumption of Eq. (1) on which the commonly used back-calculation methods are based.

Not only is the first rarefaction wave faster than expected, but deviations also exist in the entire network of characteristic lines in the wave interaction region. Given the Lagrangian flow field, a set of phenomenological characteristic lines can be constructed from the contours of fluid velocity and the longitudinal stress (the method is explained in Appendix A). As shown in Fig. 4(d), there are notable differences between this set of phenomenological characteristic lines and those given by Lagrangian analysis in the region marked by the dashed red polygon.

Note that the region marked by the dashed red polygon in Fig. 4(d) is exactly where the sample is subjected to irreversible compression. Therefore the difference in compression and rarefaction wave velocities should be closely related to the inherent asymmetry between loading and local elastic unloading, for solids beyond the elastic limit.

First, in our MD results, the rarefaction waves do indeed cause local unloading with small magnitudes. Figure 4(e) shows the longitudinal stress as a function of time for several different Lagrangian positions. After being compressed by the elastic and plastic compression waves successively, the sample is slightly unloaded upon encountering the first rarefaction waves [Fig. 4(e), arrowhead]. Second, it is during this very short period of unloading that the stress–density curves are observed to deviate significantly from the original path of loading [Fig. 4(f)].

The zone of unloading is small in magnitude, and thus the material should response elastically. The stress–strain relation in this region will be dominated by local elastic properties. On the other hand, however, for the original P(ρ) path followed by the material during loading, the slope is mainly determined by plastic deformation and the increasing fraction of the high-pressure phase. The irreversibility between plastic compression and local elastic decompression has been mentioned in this context in Refs. 21 and 22; however, we would like to emphasize here the effect of phase transition. For local elastic unloading, the fraction of the high pressure phase should remain roughly the same as a result of hysteresis. Thus, the change in phase ratio should contribute to the stress–density relation of loading, but not to that of local elastic unloading.

The behavior of the iron sample during local elastic unloading is illustrated in Fig. 4(g) by “wave-free” simulations. The unloading phase is exaggerated a bit for better visualization. A small sample of single-crystalline iron containing 2 × 10⁵ atoms is deformed along the [111] crystalline direction at ε̇=±1010 s⁻¹. To simulate the loading–unloading–reloading process, here, the sample is first compressed along the z axis to a certain extent and then stretched by 2% in volume in the z direction, and then the compression continues until the final pressure is reached. By exaggerating the unloading phase, 2% expansion in volume is still small enough to keep the unloading process elastic. Irreversibility in the stress–density curve is thus obvious. At the same time, the hysteresis of phase transition is also evident. The ratio between the low- and high-pressure phases is presented in Fig. 6 in Appendix B, which shows hysteresis in accordance with the longitudinal stress.

This kind of irreversibility affects the overall ramp loading flow field in a subtle way. On the one hand, the overall magnitude of the deviation in P(ρ) is actually rather small—less than 1 GPa as in Figs. 4(e) and 4(f)—and appears as small fluctuations in Fig. 3(g). However, the slope dP/dρ of the release curve differs considerably from the loading path, and it is the derivative dP/dρ that determines the wave velocity. Therefore, a general argument could be that two different sound velocities can actually arise from a seemingly one-to-one relation between stress and density in the MD simulations.

Our final example concerns ramp compression of polycrystalline iron (Fig. 5). Because the probability of defect nucleation is too small with the limited size of th MD sample, the onset of plasticity is delayed until 15 GPa in the previous single-crystalline case. Such an abnormally high elastic limit may not reflect the correct physics. For example, it makes plasticity and phase transition occur almost simultaneously.

Here, with nanocrystalline samples, MD simulation shows that plastic deformation begins at about 5 GPa owing to grain boundary sliding. The onset of plasticity thus separates from the phase transition—the stress–density curve obtained by wave-free simulation shows two “steps” at around 5 and 30 GPa, respectively [Fig. 5(d), dashed curve]. The stress and density inside the ramp loading sample display the same features, although the fluctuation in the MD simulation here is a bit larger than in the previous single-crystalline case because of the randomly initialized grain configurations.

For the back-calculation results, there exist similar kinds of errors as shown above, owing to the irreversibility introduced by plasticity and phase transition. The undulations in the c_L(u_p) relation are not correctly captured, and the resulting P(ρ) curve shows deviations of the order of 5 GPa—which exceed the inherent MD fluctuation range by many times the standard deviation. More importantly, these errors in the final back-calculated result represent a qualitatively different behavior of the stress–density curve. Two step-like features corresponding to plasticity and phase transition are greatly compromised—the second one has even vanished completely. Without being aware of this kind of back-calculation error, the conclusion one may draw from the resulting P(ρ) relation could be that under such a high strain rate, the structural phase transition would happen smoothly, which could be an incorrect interpretation. This observation may be particularly interesting, because one of the major goals of quasi-isentropic ramp loading experiments is to probe high-pressure phase transitions in solids.

According to our explanations, the ultimate cause of the deviation is related to half-way unloading of the irreversibly compressed solid sample, during which the stress–density relation deviates from that during loading. We would like to discuss whether this kind of unloading can be avoided, just like the undesired shocks, by using more sophisticated designs for the piston acceleration history and the target geometry.

By fine-tuning the piston acceleration history only, it would be challenging to avoid unloading completely if the sample material shows significantly nonmonotonic c_L(ρ) during compression. For a given sample thickness, after entering the region with dc_L/dρ < 0, the acceleration of the piston needs to be slow enough that those compression waves with very low sound velocities are not overtaken by subsequent compression waves. Also, to avoid unloading of the sample, we must guarantee that the rarefaction waves (generated by earlier elastic compression waves) should arrive with even lower intensity. This means that in the previous elastic loading phase, the piston must accelerate more slowly. However, too slow an acceleration, for a given sample thickness, is inherently in conflict with the requirement to avoid reverberation—i.e., the rarefaction wave should not return to the piston surface before the loading finishes, as has been emphasized repeatedly in the literature.32 Therefore, avoiding unloading sets strict constraints on the maximum pressure that can be achieved before reverberation. For example, for the case of single-crystalline iron studied here, this upper limit may only slightly exceed the elastic limit, undermining the significance of ramp loading greatly.

However, on the other hand, the strength of the rarefaction waves causing unloading may be greatly reduced by covering the rear surface of the sample with a window layer. This is a natural result of impedance mismatching. Since the transparent window materials (e.g., LiF and quartz) certainly have much higher impedance than the vacuum, the rarefaction waves generated on the interface upon shock/compression wave transition would be much weaker than those in the case of a free surface. This point has also been mentioned in Ref. 21, where it was noted that no local elastic unloading is observed for ramp loading samples with LiF windows. Although introducing the window material brings in additional uncertainties concerning the reference equation of state and the issue of correcting the refractive index under compression, it would at least help alleviating the problem of back-calculation studied here by reducing the effects of unloading.

Admittedly, the MD simulations in this study are performed at much higher strain rates than those in typical magnetically or laser-driven ramp loading experiments. The typical strain rate in this study is around 10¹⁰ s⁻¹ owing to computational limitations, which is about 20 times greater even than that in laser-driven experiments.33 However, this strain rate is still much smaller than that at the shock front (∼10¹² s⁻¹), and therefore the compression may still be regarded as quasi-isentropic, as suggested in Ref. 34. We have to use such a high strain rate because we want to simulate the entire process of ramp wave propagation and free-surface acceleration. The sample size itself is proportional to the ramp loading time. Therefore, the required computational power is inversely proportional to the square of the strain rate (or even the fourth power if the cross section is scaled accordingly). Much lower strain rates can be achieved if we choose to simulate only a small fraction of the sample, i.e., to perform simulations on homogeneous deformation of a small-size system. However, the effects related to the propagation of compression and rarefaction waves, which turned out to have strong influence on the validity of back-calculation, cannot be exhibited unless the entire ramp loading process is simulated as we have done here. Our proposed mechanism for the deviations of the back-calculation methods is related to the intrinsic irreversibility related to elastic–plastic response and hysteresis for first-order phase transitions—and these effects should still exist at lower strain rates.

Another issue related to strain rates is the possible rate dependence of the material’s stress–strain response. Since the ramp wave sharpens by itself, different strain rates are experienced by different parts of the sample. Thus, if the stress–strain relation is strongly rate-dependent, a single-valued stress–strain curve should not exist at all. We note that the broadening of the stress–strain relation in Fig. 3(g) [and Fig. 5(d)] is mainly caused by this kind of rate effect. However, in this study, we have tried to minimize rate effects by choosing an appropriate material and loading direction, and using a carefully designed piston acceleration history. Therefore, the broadening remains small enough compared with the deviation caused by irreversibility in which we are interested.

In this paper, we have examined the data interpretation process of ramp loading with the help of MD simulation. Our results indicate that the commonly used backward calculation methods may generally exhibit deviations when applied to solid samples, showing abrupt elastic-to-plastic transition and/or structural phase transitions. These kinds of transitions make the stress–strain response no longer reversible between loading and unloading. As a result, compression and rarefaction waves in a ramp-loaded planar sample may travel at different sound velocities, which could essentially affect the back-calculation process. The evidence presented in this study could serve as a reminder about the complexities associated with the back-calculation methods themselves, especially when these hydrodynamics-based back-calculation methods are to be applied to situations with sharp elastic–plastic transition and/or structural phase transitions.

Acknowledgment. We thank Dr. Xiaoxi Duan, Dr. Zhebin Wang, and Dr. Yin Zhang for helpful discussions.