Dynamic compression, including adiabatic shock compression (ASC) and quasi-isentropic compression (QIC), is of significant interest in inertial confinement fusion,1–6 laboratory astrophysics,7,8 and other high-pressure environments.9,10 In general, QIC is a process that lies between ASC and the conceptual nondissipative isentropic compression (ISE), but is closer to ISE in most cases.11–19 By contrast, ASC produces massive amounts of heat and entropy, easily leading to melting.12,19 The heating caused by plastic deformation usually echoes the increase in entropy in these types of compression and significantly affects the states of solids along the compression path. On the other hand, the microstructure that is created by plastic deformation and carries plasticity exhibits substantial dissimilarities between QICs with different strain rates and ASC.20,21 However, despite its importance, there is still only limited understanding of the relationship between plasticity and entropy production in dynamic compression.

The temperature increase ΔT in QIC can be decomposed into two components: the isentropic temperature increase ΔT_ise and the excess temperature increase ΔT_excess. The former is caused by a reduction in atomic spacing and is reversible upon decompression,22,23 and it should be path-insensitive. The latter, which echoes the entropy increase, is irreversible upon decompression and can be path-sensitive (e.g., rate-sensitive).24,25 It is generally supposed that ΔT_excess is produced in solids by thermomechanical conversion from plastic work, which is associated with evolution of strength and microstructure. It is also generally accepted that the yield strength of a given metal may undergo an abrupt increase when the strain rate increases beyond 10⁶/s.26–33 Such an increase in strength can be described by semiphenomenological models such as the Zerilli–Armstrong model,34 the Preston–Tonks–Wallace model,35 and the Steinberg–Guinan model.36,37 Characterization of microstructure evolution in dynamic compression have benefited very much from recent progress in dynamic x-ray diffraction (DXRD) technique and atomistic simulations. Phenomena including the transition from dislocation-dominated to twin-dominated plasticity and variations in dislocation and pattern as the strain rate increases in QIC have been observed both numerically and experimentally.20,21,38,39 On the other side, a temperature rise can have a significant impact on the properties of solids, including changes in the elastic modulus,28,40 the initiation of phase transitions under extreme conditions,27,41,42 and inevitably changes in strength.26,28,30,43 In addition, ΔT_excess under compression might cause an additional thermoelastic volumetric change,32,44,45 leading to another coupling in solids between the isotropic and deviatoric responses, described by the equation of state and the strength model, respectively. The understanding of temperature increase is crucial to obtaining an understanding of these couplings and predicting the variation of related properties.18 Unfortunately, it is not feasible to infer ΔT_excess from temperature data, owing to both experimental and theoretical challenges.

To estimate ΔT_excess and quantify entropy production, one method is to use elastoplastic hydrodynamic simulations with a continuum model (the Taylor–Quinney factor β_int)31,32,46–48 or discrete dislocation dynamics,49 which assumes that solids convert plastic work to heat in a constant ratio.46 Although the typical range of β_int is 0.9–1.0 in elastoplastic hydrodynamic simulations of dynamic compression, β_int has been shown to vary with factors such as strain, strain rate, and internal material structure in deformations in the absence of hydrodynamic pressure.50–54 In particular, Lazicki et al.27 have recently used DXRD to determine solid carbon phases in QIC and have conjectured that β_int for solid carbon may vary between 0.5 and 0.9. A similar technique has also been utilized to assist the development of strength models.55 Several atomistic simulations have shown that β_int for plastic deformations under ambient conditions varies with crystal structure.56,57 Despite this limited progress, there is still a lack of experimental knowledge of β_int at extremely high strain rates (>10⁵/s), hindering understanding QIC in terms of entropy production.

From the perspective of statistical physics, entropy production is related to heat flow between the atomic vibrational subsystem and the configurational subsystem.58–64 The entropy increase contains contributions from both of these subsystems. The former is related to plastic heating, reflecting the portion of plastic work converted into heat, β_intW_p. The latter should be correlated with the pattern of microstructures, which also store a portion of plastic work as what is known as cold energy, (1 − β_int)W_p.57,61,65 The interplay between plastic heating and microstructure evolution is a very unclear dynamic process. In general, the free energy of a system has the form U − TS, in which U is internal energy and S is entropy. Plastic heating and cold energy both have contributions from changes in both U and TS. Berdichevsky66–68 has pointed out that the lack of a plasticity model with satisfactory predictive capability may stem from the fact that our understanding of entropy for microstructures, which should serve as a key thermodynamic parameter, is very limited. Rao et al.69 and Xiao et al.70 have implemented the idea mentioned above to develop a two-temperature strength model for glassy systems. However, there is a lack of any direct link between microstructure characterization and thermodynamic entropy.

Our work reported here aims to provide thermodynamic insights into plastic deformation in QIC. For simplicity, a material that is isotropic, namely, metallic glass (MG), is taken as the model. We establish a partition of plastic work into heating and stored energy, quantify entropy production, derive its relationship with plastic work, and establish a partition into vibrational and configurational parts. The remainder of the article is organized as follows. Section II discusses the definition of a quasi-isentrope, on the basis of which simulation methods are developed and quasi-isentropic compression is demonstrated. Section III provides a detailed description of our analytical framework. Calculation of entropy production and its relationship with plasticity in nonequilibrium molecular dynamics (NEMD)-QIC are presented in Sec. IV. Section V presents conclusions.

A typical MG, Cu₅₀Zr₅₀, is chosen as the research model for three reasons. First, the isotropic nature of an MG allows simple mechanical analysis of, for example, strain decomposition. Second, size effects in molecular dynamics simulations of MGs are relatively smaller than in crystalline materials; see, for example, cases 1 and 2 in Ref. 71 and the results provided in the supplementary material. These small size effects in MGs are due to the fact that the mechanical interactions of an MG’s structural units decay as 1/r^∼2 rather than the 1/r decay exhibited by dislocations in crystalline materials.72 Third, the free energy and entropy during heating have been well studied in Ref. 73, providing a benchmark for our calculations. A small MG sample (S-Sample), containing 4000 atoms (dimensions 4 × 4 × 4 nm³), is prepared by rapid quenching from 2000 to 0.1 K in 200 ns. Subsequently, the sample is heated for 100 ps to 600 K, followed by relaxation at 600 K for 50 ns, after which it is cooled gradually to the target temperature. The glass transition temperature T_g of the sample is ∼620 K, as reported in Ref. 73. The atomic interaction in the MG is described by the embedded atom method (EAM) potential.74 All simulations are carried out using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) package.75,76 Timesteps are 1 fs for all simulations in this article. In contrast to the conventional usage in mechanics, positive values denote compressive stress and strain in this work, whereas negative values indicate tensile stress and strain.

Nonequilibrium molecular dynamics (NEMD) simulation are is employed to create quasi-isentropic compression of the MG. In particular, as well as ramp compression,21,24 another two artificial compressions are also simulated and analyzed, as depicted in Fig. 1(a), to exclude other potential dissipative mechanisms and to verify the validity of the proposed framework. The first compression method is uniaxially adiabatically compression (UC-QIC) of the S-Sample with controlled strain rate and no pressure and temperature control (NVE). UC-QIC results in a uniaxial strain state (ɛ_n ≡ ɛ_zz = Δρ/ρ > 0 and ɛ_t ≡ ɛ_xx = ɛ_yy = 0, where x, y and z are Cartesian coordinates, n represents the compressional direction, and t is orthogonal to n). Periodic conditions are applied to the x, y, and z directions. The computational cost of UC-QIC is relatively low, and therefore the effects of both size and strain rate in uniaxial compression are tested by UC-QIC cases. The results show that both the size effect and the effect of strain rate >10⁸/s on yield strength and temperature rise in UC-QIC are negligibly small (see Fig. S7, supplementary material).

The second compression method in this study is based on NEMD simulation of ramp wave propagation (NEMD-QIC) in a long sample (L-Sample), also leading to a uniaxial strain state. NEMD-QIC simulates experimental QIC. The L-Sample, containing 3 × 10⁷ atoms, has dimensions 20 × 20 × 1200 nm³, obtained by repeating the S-Sample 5 × 5 × 300 times. A flyer, whose particle velocity increases from zero to 1200 m/s in a rise time t_rise, moves from the lower end toward the upper end of the L-Sample, generating a ramp wave propagating toward the upper end within the sample along the z direction. Periodic conditions are applied to the x and y directions. NEMD-QIC is characterized by t_rise, during which time the pressure and density increase smoothly to a peak, as described in Refs. 21 and 24. t_rise is related to the strain rate according to γ̇∝1/trise. The rise time and propagation distance (L ∝ Ct_rise, where C is the wave velocity) achieved in NEMD-QIC are at least an order of magnitude lower than those in experiments.

The third compression method is “pseudo” isentropic compression (pISE) achieved through isotropic adiabatic compression of the S-Sample without pressure and temperature control (εn=εt=13Δρ/ρ>0), with periodic conditions along the x, y, and z directions. Note that NEMD-QIC results in spatial thermodynamic gradients (∇ρ, ∇T, ∇P, …) along the wave propagation direction n. By contrast, these gradients are not present in UC-QIC and pISE.

The three compression methods can lead to several potential candidate dissipative mechanisms,31,32 with plasticity possibly being dominant. The first mechanism is linked to the uniaxial strain state, which results in an increase in deviatoric strain (γ ≡ ɛ_zz − ɛ_yy > 0) and subsequently yields plastic deformation in solids. The second is based on the familiar concept of artificial bulk viscosity77 acting in the presence of a pressure gradient. The third is heat conduction due to the presence of a temperature gradient. All of these mechanisms can produce excess temperature and entropy increase. In contrast to NEMD-QIC, for UC-QIC, only plastic deformation is a candidate dissipative mechanism, owing the absence of gradients. Note that pISE provides a nondissipative compression benchmark. The term “pseudo” implies a small deviation from a perfect isentrope. Figures 1(b) and 1(c) demonstrate that the P–ρ–T paths in NEMD-QIC and UC-QIC almost completely overlap with that in pISE in early-stage compression, indicating that the behavior in this regime is consistent with an ideal isentrope. Dissipation may begin to appear in NEMD-QIC and UC-QIC at ρ/ρ₀ ≈ 1.1, where, as we will show later, yielding occurs, resulting in a significant deviation from pISE as depicted in Fig. 1(b).

We must stress that a comparison of the results for NEMD-QIC and UC-QIC in Fig. 1(b) may lead to an illusion that the pressure and temperature gradient contribute to the temperature increase in NEMD-QIC. However, we argue that the temperature increase induced by heat conduction in ramp compression at low pressure is negligibly small. The thermal conductivity κ of our MG sample is calculated by the Green–Kubo method,78 which relates the ensemble average of the autocorrelation of the heat flux J to κ through κ=V/3kBT2∫0∞J0Jtdt. The autocorrelation of J is computed for the small sample under the microcanonical ensemble. The results show an increase in κ from ∼0.6 W/mK at 0 GPa to ∼1.5 W/mK at 40 GPa (T = 300 K). The heat flow q caused by the temperature gradient is estimated as q = κ ∇T ≈ κT_max/Ct_rise, where T_max is the maximum temperature increase achieved in a Lagrange element and C is, for example, 5000 m/s. Therefore, the temperature increase caused by heat conduction in a time duration t_rise can be estimated: ΔT_c ≈ κ(T_max/Ct_rise)(1/ρCC_V). Under the assumption that the maximum temperature increase is not too high (e.g., 200 K), ΔT_c < 0.1 K is negligibly small for t_rise > 10 ps. Obviously, ΔT_c as a fraction of the total temperature increase is ΔT_c/T_max < 0.001 for t_rise > 10 ps and ΔT_c/T_max < 0.000 01 for t_rise > 1 ns. If the density is not high enough to lead to an order-of-magnitude increase in κ,79 then heat conduction in ramp compression can be ignored. In this study, the rise times vary from 10 to 400 ps, for which gradients may not be manifested.

Before proceeding any further, we first present a conceptual discussion on the deviation from isentropic behavior. Figure 2 illustrates the increase in temperature with density under compression. It must be emphasized that, unlike plastic deformation with no pressure, a “perfect” reference isentrope can be established in dynamic compression from the initial state, along which there is no irreversible nonaffine deformation or microstructural change, an example being the conditions expected in pISE. This idealized isentropic path is represented by the solid blue line in Fig. 2. In dynamic compression, isentropic compression consumes the least energy. By contrast, in the case of QIC, plastic deformation is inevitable, leading to a waste of energy in dissipation. The quasi-isentrope yielded by experiments is shown by the red solid line in Fig. 2. A state C along the quasi-isentrope can be described by [TC,ρC,γaC≠0,…], where T is the temperature, ρ is the density, and γ_a measures the inevitable microstructural change. For a given state C, it is possible to define a corresponding reference isentropic state A with identical density, by [TA,ρA=ρC,γaA=0,…]. Before significant dissipation occurs, such as in the case of elastic compression, the temperature at state A, T(A), can be assumed to be equal to the temperature at state C, T(C). Along the perfect isentrope, the temperature is approximated by integrating over a small duration using the expression dTise|AB=−TAγgAdlnρA/ρB, where γ_g is the state-dependent Grüneisen parameter. However, we can only observe state C, with detailed atomic trajectories, in simulations but very limited observables in experiments. Hence, an isentropic path starting from state C is defined as dTise|CE=−T(C)γg(C)dln(ρA/ρB). To reach state D on the quasi-isentropic path, isochoric heating is required to compensate for the dissipation that occurs during compression from state E. By integrating dT|ED=dT|CD−dTise|CE, the amount of dissipation can be quantified.

The recovery of the reference isentrope is facilitated by the accumulated dissipation. In hydrodynamic simulations, the determination of the differential form dT|CD for temperature change along QIC can be expressed as44dT|CD=−Tγg(ρ,T)dlnρ0ρ+TdSCV(ρ,T),(1)where C_V is the isochoric heat capacity. The first term on the right-hand side of Eq. (1) is identical to dTise|CE. The second term, on the other hand, pertains to the dissipation contribution. It is worth noting that the black line, which results from integrating dTise|CE=−T(C)γg(C)dln(ρA/ρB), represents an instant isentrope rather than the reference isentrope. The difference between the temperature increase along the instant isentrope and that along the reference isentrope is given by ΔTerr=ΔT|CE−ΔT|AB=∫TAγgA−TCγgCdlnρA/ρB, which may increase as the compression proceeds. In general, ΔT|BE>ΔT|AC is expected, because of the difference between states A and C, raising subtle arguments about the reference isentrope. This difference is negligible when the temperature is as small as <600 K in our simulations.

The reference isentrope should be a fundamental property of a solid material that remains unchanged under varying compression methods, as long as there is no path-dependent phase transition. In our study, we quantify dissipation by rewriting Eq. (1) asdS=CV(ρ,T)dTT+γgρ,TCVρ,Tdlnρ0ρ,(2)and the entropy can be directly obtained by integrating this equation along the compressional path. The excess temperature increase appears as ΔT_excess = ∫TdS/C_V. The recovery of the reference isentrope is achieved by ΔT_ise = ΔT − ΔT_excess.

Our framework as illustrated in Fig. 3 aims to mine the data, from the perspective that it is the plasticity that contributes to dissipation and microstructural change. First, the fraction of the plastic work converted into heat, as the Taylor–Quinney factor, is calculated. As discussed in Sec. II C, we already have the plastic heat Q=∫ρ0ρTdS=∫ρ0ρCVdTexcess on hand according to Eqs. (1) and (2). To quantify plastic work, the Taylor strain condition is utilized to decompose uniaxial strain ɛ_n into an elastic contribution εne and a plastic contribution εnp:31,44dεne=dρ3ρ+dσn−dσt3G(ρ,T,γa),dεnp=dεn−dεne,(3)in which only the state-dependent shear modulus G is employed. Density, normal stress σ_n = σ_xx, and transverse stress σ_t = 0.5(σ_yy + σ_zz) are from molecular dynamics (MD) data. The strain in Eq. (3) is directly computed from MD simulations as dɛ_n = dρ/ρ. Plastic heating is not considered. According to the Lindemann-like rule,80 the linear thermal expansion caused by additional heat produced by plastic work is negligible as estimated by 0.75 × 10⁻²ΔT/T_g (<0.005). Plastic work dWp=(σn−σt)dεnp and β_int = Q/W_p then appear as functions of εnp.

Second, we try to establish a connection between dissipation and microstructural change. As well as the part of the plastic work converted into dissipation, another part is stored as internal energy for microstructural change (without TS term): for example, defect energy in crystalline metals. However, the microstructural change also contributes to entropy increase by its disorder. We notice that the entropy S obtained by integrating Eq. (2) along the compressional path can be decomposed into vibrational and configurational components.59,73S directly quantifies the dissipation. Following the Boltzmann relation, we have S=kBlnΩ∼kBlnΩvΩc=kBlnΩv+kBlnΩc=Svib+Sconf, in which Ω and Ω are the number of vibrational and configurational states, respectively, and S_vib and S_conf are the vibrational and configurational entropies, respectively. To achieve the decomposition, we further notice that the Helmholtz free energy F containing only a vibrational contribution can be sampled by the Frenkel–Lad method,73 the details of which will be discussed in Sec. III B. Thus, S_vib = (U − F)/T (where U is internal energy) is obtained. Then, the decomposition of ΔS can be approximated by ΔS_conf = ΔS − ΔS_vib. It should be emphasized that there exists a correlation between ΔS_conf and microstructural change. Unfortunately, accurate quantification of this correlation remains challenging.67,68 Borrowing an idea from statistical physics, here we utilize a local structure analysis method, Voronoi polyhedron analysis,81 to reveal atomic structural changes. The mathematical quantification of microstructure is thus made via local fivefold symmetry (LFFS)82 and Shannon entropy83 so as to link with ΔS_conf.

Specifically in the case of QIC, the above framework has the capability to accurately determine the increase in entropy, its quantitative relationship with plastic work, and its correlation with changes in microstructure. Notably, it can uncover the configurational contribution to entropy, thereby enabling a thermodynamic characterization of microstructural change, a long-sought goal in materials science.

In the above framework, the shear modulus G, heat capacity C_V, Grüneisen parameter γ_g, free energy, and other thermophysical properties are assumed to be state-dependent. This assumption is also widely utilized in elastoplastic hydrodynamic simulations when rough estimates are made.18,32 For example, ργ_g = const is usually used.84 In our study, state-dependent properties are calculated directly from atomistic simulation of a small sample with only 4000 atoms instead of by employing phenomenological laws. Irreversible plastic deformation may impact the physical properties of an MG. Six samples are systematically prepared by subjecting the initial S-Sample to different levels of plastic deformation through shearing, followed by an un-shearing to a zero-stress state and subsequent NPT relaxation. The boundary conditions are consistent with those utilized in UC-QIC and pISE. As a result of the preparation procedure, residual strain is induced in the samples, which is quantified through atomic shear strain calculations.85 Taking the undeformed S-Sample as the reference configuration, the average atomic shear strains of these samples are determined by the Open Visualization Tool (OVITO)86 to be γ_a = [0.0, 0.0284, 0.0366, 0.110, 0.137, 0.166], respectively. γ_a = 0 corresponds to the undeformed sample. Each sample is subsequently compressed and heated to map the P–T subspace from 0 to 50 GPa and from 300 to 1200 K in pressure steps of 5 GPa and temperature steps of 50 K. The mechanical modulus is obtained by linear fitting of the strain–stress curve under small compression and shear deformations. The isochoric heat capacity C_V = −∂U/∂T and the Grüneisen parameter γ_g = (V/C_V)dP/dT are obtained by fitting U–T and P–T curves, respectively, under isochoric heating (δT = 10 K) from a given (P, T) state. The Helmholtz free energy F and the vibrational entropy S_vib = (U − F)/T are sampled by the Frenkel–Lad method.73,87 Then, functions of the form H(ρ,T,γa)=h0[1+h1ρ/ρ0+h2(ρ/ρ0)2](1+h3T+h4T2)(1+h5γa+h6γa2) for all properties are prepared for use. H(ρ, T, γ_a) assumes a linear deviation from ambient conditions with small density, temperature, and structural changes. h_i (i = 0, …, 6) are fitting parameters and are listed in Table I. The detailed calculation of F is discussed below.

The Frenkel–Lad method samples the Helmholtz free energy containing only the vibrational entropy contribution of the MG.73 The basis of this method is calculation of the free energy difference between a state of interest and a reference state for which an exact solution exists. In this study, an Einstein crystal in which atoms do not interact with each other via an atomic potential, is chosen as the reference state. The free energy difference between a state of interest and this reference state is calculated by the construction of a reversible path connecting a sequence of equilibrium states between the two thermodynamic states. By the introduction of a coupling parameter λ ranging from 0 to 1, the Hamiltonian of the system can change from the Hamiltonian of the Einstein crystal reference state (EC) to the Hamiltonian of the real system (R):Uλ=1−λUEC+λUR,(4)where U is internal energy. For the system in the NVT ensemble, the free energy F is a function of λ:Fλ=−kBTlnQλ=−kBTln∫e−U(λ)/kBTdr.(5)

The derivative of the free energy with respect to λ can be obtained asdF(λ)dλ=1Q(λ)∫e−U(λ)/kBT∂U∂λdr=∂U∂λλ.(6)

The free energy difference is obtained by integration of this differential equation:ΔF=∫01∂U∂λλdλ=∫01UR−UECdλ.(7)

The Hamiltonian of the Einstein crystal reference system has the harmonic formHEC=∑i=1Npi22m+12mω2(ri−ri0)2.(8)

The free energy of the Einstein crystal can be solved for analytically asFECN,V,T=3NkBTlnℏωkBT,(9)where ω is the oscillator frequency, and the harmonic potential with spring constant (mω²) is usually utilized in simulations to describe the Hamiltonians of a reference Einstein crystal. For the MG of interest here, corrections for re-centering of mass need to be considered:FECN,V,TkBT=∑i=Cu,Zrxi⁡lnβkΛi22π3/2−1Nlnβk2π∑i=Cu,Zrxiμi23/2−1NlnV+2.5NlnN−1NlnNNmol,(10)where x_i is the mole fraction, μi=mi/∑i=Cu,Zrmi, Λi=βh2/2πmi are the wavelengths of copper (i = 1) and zirconium (i = 2) atoms, and k is the spring constant. With F_EC and ΔF, the Helmholtz free energy of the MG is obtained as F = F_EC + ΔF, and also the vibrational entropy S_vib = (U − F)/T.

Figure 4 shows several results for the calculated thermophysical properties according to Table I. To discuss the rationality of these results, a density-sensitive parameter Hρ=1/H∂H/∂ρ, temperature-sensitive parameter H_T = (1/H)∂H/∂T, and deformation-sensitive parameter Hγa=(1/H)∂H/∂γa of the thermophysical properties are defined. At the ambient states with small irreversible deformations, it is clear that all these properties are sensitive to temperature, all except heat capacity are sensitive to density, and only shear modulus is sensitive to irreversible deformation. These expected scenarios are shown in Fig. 4. Figures 4(a) and 4(c) show that shear modulus and vibrational entropy are sensitive to density and temperature. Figure 4(c) also indicates a clear separation of glass and liquid states on the map, and it suggests that the glass transition increases with density as expected. The vibrational entropy is not sensitive to plastic strain, because the sampling method is irrelevant to configuration, as confirmed by its small Hγa according to Table I. The deformation-insensitive isochoric heat capacity as shown in Fig. 4(d) also supports this scenario. In particular, Fig. 4(b) demonstrates that γ_g decreases with increasing density. This trend can be explained by the fact that t phonons also exist in the MG, despite its lack of long-range order.88 Additionally, ργ_g = const is not applicable for an MG over wide density and temperature ranges. This is reasonable, considering that ργ_g = const is only a phenomenological law inferred from ASC experiments that access limited regions of density and temperature in ρ–T space. Because γ_g plays a key role in determining the isentropic temperature increase according to Eqs. (1) and (2), our method for obtaining γ_g at high pressure and high temperature may be applicable for interpretation of QIC experiments. For instance, γ_g at high density and pressure can be accurately calculated by first-principles calculation or machine-learning-enabled MD with quantum accuracy.

With the help of the framework established above, the quantification of dissipation by excess temperature increase is illustrated in Fig. 5(b), enabling decisive exclusion of alternative dissipative mechanisms apart from plasticity. The analysis shows that the process of pISE results in a minimal ΔT_excess, indicating a near-isentropic behavior as designed. By comparing Figs. 5(b) and 5(c) in which the strength is calculated by Y = σ_n − 0.5(σ_yy + σ_zz), ΔT_excess within the elastic regime in both NEMD-QIC and UC-QIC is negligible, with NEMD-QIC displaying gradients ∇T and ∇P. As discussed in Sec. II B, heat conduction does not lead to significant ΔT_excess as t_rise > 10 ps, as demonstrated in the elastic regime [Fig. 5(b)]. Also, the contribution of the bulk viscosity to dissipation in QIC is small when the magnitude of ∇P is much lower than that at the shock front. This suggests that the contributions from both heat conduction (as discussed in Sec. II B) and bulk viscosity to dissipation in NEMD-QIC are inherently negligible in the elastic regime. After a short period elapsing from yielding [as can be seen by comparing Figs. 5(b) and 5(c)], ΔT_excess in NEMD-QIC and UC-QIC begins to increase significantly, while the thermal conductivity and bulk viscosity are not greatly changed by the yielding transition. These findings convincingly demonstrate that plastic deformation dominates dissipation in a MG subjected to QIC, as well as in crystalline solids, despite the distinct microscopic mechanisms involved.31,32

By subtracting ΔT_excess, the reference isentrope is recovered as a substantial and unchanged property of MG, as illustrated in Fig. 5(a), in which the ΔT_ise for NEMD-QIC, UC-QIC, and pISE overlap with each other. In our simulations, we argue that no phase transition occurs in MG. Unlike crystalline solids, MG structures lack long-range order (LRO) and exhibit short-range order changes upon compression, as revealed by Voronoi analysis in terms of local fivefold symmetry (LFFS) increase.74 The smooth changes in LFFS without disruption of disordered structures in MG suggest the absence of a phase transition. The reference isentrope recovered from all NEMD-QIC and UC-QIC simulations close to the pISE path is thus rational. In addition, it is clear that ΔT_excess in NEMD-QIC is higher than that in UC-QIC over the whole plastic deformation, as illustrated in Fig. 5(b). This discrepancy is likely due to the unique deformation protocol employed in UC-QIC, which includes periodic and uniform rescaling of atomic positions followed by a structural relaxation. This rescaling results in an unexpected portion of strain energy being directed toward the microstructural subsystem during the plastic deformation regime. A structural comparison between NEMD-QIC and UC-QIC demonstrates this significant energy flow (see Fig. S11, supplementary material). Thus, UC-QIC substantially deviates from a real QIC process represented by NEMD-QIC, where plastic deformation and structural change should be taken into account.

An unexpected finding is revealed in regard to the Taylor–Quinney factor β_int, which quantitatively measures the relation between plastic work and excess temperature. The results indicate that β_int is insensitive to strain rates but sensitive to strain. All β_int vs plastic strain curves exhibit the same behavior as the plastic strain increases, converging in the vicinity of one another, as depicted in Fig. 6(b). This outcome contradicts experimental work conducted by Nieto-Fuentes et al.,53 who found that β_int decreases with strain rates that are about six orders of magnitude lower than that in our simulations for a variety of metals. Additionally, an incubation period in which β_int remains near zero is observed shortly after yielding. During this period, the plastic work is primarily consumed by microstructural change rather than converted to heat, in parallel with Figs. 5(b) and 5(c). After the incubation period, β_int increases monotonically with plastic strain and follows an empirical master curveβint(εp)=0,εp<εpc,β0lnεpεpc,εp≥εpc,(11)where εpc is the ending of the incubation period, and β₀ is a fitting constant. On setting β₀ = 0.45 and εpc=0.03, Eq. (11) is found to fit all curves in Fig. 6(b) with satisfactory accuracy. The observed trend strongly suggests that β_int is sensitive to structural change in the early stage of plastic deformation. As compression proceeds to the late plastic deformation stage, extrapolation of Eq. (11) suggests that β_int reaches a saturation value of unity at ρ ≈ 1.8. This value is consistent with that typically used in hydrodynamics simulations.31,32

As an indirect measure of entropy production, quantification of ΔT_excess reveals an unexpected rate-insensitive to rate-sensitive transition in entropy production. This is demonstrated in Fig. 6(a), where, for t_rise > 100 ps, NEMD-QIC produces ∼100 K excess temperature at ρ/ρ₀ = 1.29, while the relative isentropic temperature is about 100 K. As t_rise decreases to 40 ps, the excess temperature increases to around 120 K. As a comparison, adiabatic shock compression produces ∼370 K, which is approached by NEMD-QIC as t_rise further decrease. Surprisingly, the presence of rate-insensitive entropy production challenges the conventional understanding obtained from low-strain-rate mechanical tests that plastic heating behavior is rate-sensitive.51,53 Our simulations suggest that rate-sensitive ΔT in a single-phase solid subjected to high-pressure compression can only be achieved at extremely high strain rates, even though the strength is rate-sensitive over a wide range of strain rates, while a trend of increase in both yield strength Y and W_P at ρ/ρ₀ = 1.29 with decreasing t_rise is evident.

Phenomenologically, with the quantification of β_int, the transition from rate-insensitive to rate-sensitive entropy production can be roughly explained by the breakdown of Y/G ≈ const instead of Y, where Y is the yield strength and G is the shear modulus at yielding. This can be argued by a simple scaling. It is noted that ΔT_excess should be proportional to β_int × W_P/ρ. This scaling is expressed asΔTexcess∝βint×WPY×Gyρy×ρyρ×YGy∝βint×WPY×Gyρy×YGy.(12)

For all simulated cases, the increase of W_P with strain rate can be attributed to increasing Y (W_P/Y ≈ const for all t_rise). Also, G/ρ ≈ const for all t_rise holds because of the linear dependence of the shear modulus on density at low pressures. Therefore, the following simpler scaling is proposed:ΔTexcess∝βint×YGy.(13)

The rate-insensitivity requires that β_int × Y/G be independent of t_rise. As previously discussed, it is found that there is minimal dependence of β_int on t_rise in all cases, as illustrated in Fig. 6(b). As a result, the magnitude of entropy production can be attributed to variations in Y/G. This relationship is illustrated in Fig. 6(a), where a transition from rate insensitivity to rate sensitivity of Y/G can be observed, as well as a corresponding change in entropy production. Thus, Eq. (13) may explain the trend of how entropy production in QIC varies with strain rate.

The observed strain rate dependence of Y/G in MGs is in line with previous studies. According to the cooperative shear model (CSM), at low strain rates, a weak rate dependence of Y/G in MGs is expected if yielding is viewed as a fold catastrophe.89 The CSM suggests that MGs contain inherent defects that become activated and contribute to plasticity upon yielding. Measurements of the yielding strength of an MG in a QIC experiment support this hypothesis, since they show a lack of sensitivity around γ̇∼105/s.90 It has been observed that in many crystalline metals, at higher strain rates, a transition of Y/G from rate insensitivity to rate sensitivity occurs at γ̇∼108/s.31,43,91 This transition may indicate a shift in the mechanism of yielding from activation of existing defects to nucleation of new plasticity carriers such as dislocations.21,91 On the basis of these findings in crystalline metals, it is hypothesized that a similar transition in the yielding mechanism may occur in MGs.

To further investigate this hypothesis, a concept from statistical physics known as a replica is employed. Specifically, the distance from the atomistic replica for each t_rise(X) to that for t_rise = 400 ps (Y only appears as a subscript when representing a replica) is defined as DX,Y2=(1/N)∑idX,Yi=(1/N)∑i(xXi−xYi)2+(yXi−yYi)2, where x and y are atom coordinates, and dX,Yi is the spatial distance between the same atom in different replicas. In an MG, each atom is trapped in a cage formed by its neighboring atoms. On deformation, certain cages may break, serving as plasticity carriers. dX,Yi=1.0 Å is adopted as a criterion below which the behavior of the cage centered on atom i is as the same as that at yielding of t_rise = 400 ps compression, and above which its behavior deviates from that at yielding of t_rise = 400 ps compression. Similarly, dX,Yi<1.0 Å should be expected if the cages centered on atom i both break or both survive at these two rates. It is supposed that our MG is prepared with born defects. The cages regarded as the born defect may break at yielding for all rates, and thus their center atom will have dX,Yi<1.0 Å. In the scenario of the CSM, such cage breakages associated with born defects may cooperate, leading to an activation volume of ∼100 atoms. Atoms with dX,Yi>1.0 Å can then be regarded as an extra broken cage by comparison with yielding in t_rise = 400 ps compression. The term “broken” here means that the cage changes in a different way to its behavior at t_rise = 400 ps, and so it cannot sustain shear stress after yielding. We observe that the fraction of extra cage breakage remains as small as 1% in the rate-insensitive regime [Fig. 6(c)], suggesting that the yielding behavior for these rates is almost the same as that for t_rise = 400 ps. Meanwhile, sharp increases in both atom fractions with dX,Yi>1.0 Å and the distance between yielding replicas, as well as Y/G, are observed at t_rise ≈ 150 ps. This suggests that there is an extensive occurrence of extra cage breakage at high strain rates. For instance, Fig. 6(d) shows that the extra broken cages are dispersed throughout the bulk of the yielding replica in the case of t_rise = 67 ps. Note that the extra breakages are not related to born defects and will not cooperate. Therefore, the activation volume of extra breakage should be as small as several atoms, which is much smaller than that for the CSM. It is well known that small activation leads to a higher strain-rate-sensitive factor m=3kBT/σV* (where V* is the activation volume). Thus, our observation provides strong evidence for a transition in the yielding from the activation of born plasticity carriers to a mechanism aided by the nucleation of new carriers.

From the above discussions, our results indicate the surprising extrapolation that there should be minimal variation in excess temperature when the rise time is varied from 10¹ to 10² ns in typical QIC experiments. This is due to the fact that both scaled strength and β_int are insensitive to strain rate within this range, which may also hold under experimental conditions. To effectively manipulate the thermodynamic path, it is necessary to change the initial temperature of the sample or to manipulate the internal structure to control the yield strength at low strain rates in QIC experiments. Additionally, our data suggest that there may be a significant increase in excess temperature as the rise time is decreased to the nanosecond regime or lower in these experiments.92

Direct calculation of the thermodynamic entropy increase provides a quantification of dissipation in parallel with the excess temperature. In the NEMD-QIC simulations, we observe that, as well as excess temperature, ΔS at ρ/ρ₀ = 1.29 exhibits a transition from rate-insensitive to rate-sensitive behavior, as depicted in Fig. 7(b). As t_rise > 150 ps, ΔS remains constant at ∼0.6 k_B/atom at ρ/ρ₀ = 1.29, and the production rate is about 2.0 × 10⁻²k_B/atom per GPa. However, as t_rise decreases to ∼30ps, which is comparable to the adiabatic shock compression wave front width, ΔS at ρ/ρ₀ = 1.29 approaches 0.7 k_B/atom, with a slight increase compared with that for rate-insensitive compression. For comparison, we also calculate ΔS for adiabatic shock compression to ρ/ρ₀ = 1.29, with a value of ∼1.6 k_B/atom. This is significantly higher than that for all QICs, in agreement with the thermodynamic theory. This direct quantification of entropy reveals that QIC is substantially different from shock compression, even when the ramp wave width is as small as 30 ps. It should also be stressed that Eq. (2) decomposes the entropy increase into an isothermal compression contribution and an isochoric heating contribution. In the elastic regime, these two contributions cancel each other out, resulting in an isentropic process, as expected. On yielding, the positive isochoric heating contribution becomes dominant over the negative isothermal compression contribution, leading to an overall increase in entropy.

An advantage of entropy quantification is that it can reveal how the dissipated energy is partitioned. Entropy is split into two parts, a vibrational part S_vib and a configurational part S_conf, as shown in Fig. 7(a). From the perspective of the potential energy landscape (PEL),93,94S_vib arises from intrabasin vibration, while S_conf arises from interbasin exploration, which is enhanced during plastic deformation. It is important to note that S_conf is a thermodynamic measure of structural change, as has been discussed extensively in the literature.59,63,67,68 We find that the magnitude of ΔS_conf in plastically compressed MG subjected to QICs is on the order of 10⁻¹k_B/atom per plastic strain of 0.1 [Fig. 7(a)]. This is comparable to that in the melting of an MG.73 This fact suggests that MGs also undergo rejuvenation under dynamic compression, even though the increase in density leads to free volume reduction and local order increase in terms of LFFS. Our result implies that the essential structural change in rejuvenation is the diversity of local atomic structures instead of the local symmetry, which will be demonstrated later.

It should be emphasized that the increase in ΔS_vib is coupled with ΔS_conf. In a highly harmonic approximation, ΔS_vib can be estimated asΔSvib=3NlnUKρUKρ0ωρ0̄ωρ̄,where N is the number of atoms, U_K is the kinetic–vibrational energy, and ω̄ is the average vibrational frequency. With little structural change in the elastic regime, the relation UK∝ω̄ is established and leads to ΔS_vib ∼ 0. Figure 7(a) shows that ΔS_vib is relatively small in the elastic regime of all QICs. The immediate increase of ΔS_vib on yielding is attributed to vibrational mode softening caused by structural change,95 which is reflected by ΔS_conf. As compression proceeds, UK∝ω̄ may be reestablished, and saturation of ΔS_vib is observed, as in Fig. 7(a). The saturation of ΔS_vib in a real QIC is also reasonable, because ΔT_ise may totally prevail over ΔT_excess at very high density. Additionally, the increase in ΔS_conf lags behind that in ΔS_vib, but contributes more at ρ/ρ₀ = 1.29 to the overall entropy increase in our low-pressure QICs. This is in agreement with the continuous evolution of the microstructure of MGs under QIC. Therefore, it can be concluded that plastic deformation-induced microstructural change is the fundamental cause of the entropy increase in MGs subjected to QIC in a low-pressure regime.

An analysis is performed to investigate the relationship between structural change, thermodynamics, and β_int. Microstructural changes in an MG can be described by LFFS and the Shannon entropy S_Shan defined on the LFFS distribution.83,96 As shown in Fig. 8, an increase in LFFS is observed under QIC with increasing density. Note that LFFS of an MG usually decreases under normal plastic deformation without hydrodynamic pressure. This decrease in LFFS can be explained by the more liquid-like structure associated with the continuous rejuvenation and an increase in accessible inherent structures.85,97 Along the rejuvenation direction toward liquid, the configurational entropy should increase.73Figures 7(a) and 8(a) clearly imply the opposite, i.e., that the configurational entropy and LFFS all increase with plastic strain under QIC. These striking results suggest that LFFS may not be able to represent the plasticity carriers in an MG. However, Figs. 7(a) and 8(b) show a continuous decrease in SShanCu using only a copper-centered polyhedron, and an increase in S_conf as compression proceeds from yielding. This trend was also observed by Han et al.,83 who quantified configurational entropy from PEL diversity sampling and defined S_Shan in the same way as we have done. Thus, the Shannon entropy, which characterizes structural diversity, rather than LFFS, which characterizes the average local structure, is more suitable to describe the microstructural changes during plastic deformation of MGs. We also observe in Fig. 8(b) that at ρ ≈ 1.05ρ₀, pressure may dominate the structural change, causing an abrupt decrease in SShanCu.

It is quite interesting that in the rate-insensitive regime, our S_conf and SShanCu fall on almost the same line [Fig. 9(a)], suggesting their equivalence in describing the structural change of MG. The configurational entropy clearly has a thermodynamic sense, while the Shannon entropy is a mathematical description of the microstructure of our atomic model glass. Since the Shannon entropy can be considered in terms of the distribution of local environments of atoms, it is possible to use atomic force microscopy to measure the configurational entropy of an MG. A similar measurable microstructural feature, the fractal dimension of dislocation patterns, is found to increase with plastic deformation.98 Our work may help relate microstructural characterization with the thermodynamic entropy of the microstructure, thereby providing constraints for plasticity models.59,99,100

Another issue of interest is the relation between the configurational entropy and the Taylor–Quinney factor. Figure 9(b) shows that β_int increases monotonically with S_conf on yielding. This is consistent with Langer’s proposal58–61,63 that β_int depends on the interplay between the configurational subsystem and vibrational subsystem. As the compression proceeds after yielding, the configurational entropy increases with plastic strain, suggesting that more and more accessible inherent states and degrees of freedom are activated by deformation. Thus, it is speculated that plastic work becomes more easily converted into heat as the configurational entropy increases. Figure 9(b) also shows that β_int vs ΔS_conf falls onto almost the same master curve for compressions in the rate-insensitive regime. The following formula, without any physical preconditions, fits this master curve very well:βint=1−exp−ΔSconf−ΔSCab.(14)

The fitting parameters ΔS_C, a, and b are −0.044, 0.234, and 0.662, respectively. It is interesting that exp−ΔSconf−ΔSCab(b<1) is of Kohlrausch–Williams–Watts (KWW) form, which is usually used in describing relaxation of glasses.

The effective temperature T_eff is defined as the conjugate of the configurational entropy, T_eff ≡ ∂U_C/∂S_conf,59,100 where U_C is the configurational energy, or, less accurately, “defect” energy, while, obviously, T is the conjugate of the vibrational entropy, T ≡ ∂U_V/∂S_vib, where U_V is the internal energy minus U_C. We simply calculate U_C as W_p − Q = (1 − β_int)W_p. Figure 10(a) shows that U_C increases significantly after yielding and decreases slightly in late compression, while S_conf continues increasing. This suggests that the energy storage capability of the configurational subsystem of an MG is limited during QIC compression. Polynomial functions are used to fit the data in Figs. 7(a) and 10(a) (replacing relative density with plastic strain). With these polynomials, the effective temperature is calculated asTeff≡∂UC/∂εp∂Sconf/∂εp.

The trend of decreasing T_eff with plastic strain after yielding is depicted in Fig. 10(b). The extremely high T_eff on yielding can be considered to be a result of mechanical rejuvenation associated with yielding.69,70 After yielding, the fact that T_eff gradually decreases indicates that the physical aging effect starts to favor mechanical rejuvenation,24,26 which is expected in plastic deformation of MGs. T_eff becomes lower than the thermodynamic temperature at ɛ_p ≈ 0.05 and becomes negative at ɛ_p ≈ 0.09, at which point U_C starts to decrease, suggesting a unique nonequilibrium state. The entropy of the configurational subsystem becomes a nondecreasing function of energy at ɛ_p ≈ 0.09. T_eff should have a similar trend to that of T.59 The higher the value of T_eff, the more disordered should the configurational subsystem be. In this sense, the most disordered state is expected to occur at yielding, which causes an extreme mechanical rejuvenation in MGs, and the configurational subsystem evolves toward the direction of order after yielding, while the configurational entropy keeps increasing. Note that Fig. 8(a) indeed indicates that the structure of an MG under QIC evolves toward greater and greater order, characterized by increasing LFFS. This suggests that the heat capacity for the configurational subsystem is first negative and then positive under QIC. Unfortunately, to the best of our knowledge, there is no appropriate theory that can explain our results, even though we have successfully obtained the effective temperature as expected from Ref. 59.

We have presented a framework for the analysis of the thermodynamic path undergone by a metallic glass subjected to compression and to provide insight into plastic heating. Using this framework, it is possible to accurately recover the same isentrope of the initial sample from P–ρ–T paths obtained by both isotropic compression and NEMD simulations. The deviation from the isentrope is quantified using the excess temperature and thermodynamic entropy. We have demonstrated that this deviation stems from plasticity under QIC of metallic glass and quantifies the portion of plastic work being dissipated. Our results suggest that rate-insensitive entropy production under compression with strain rates achievable experimentally is of the order of 10⁻²k_B/atom per GPa. However, at higher strain rates, a transition from rate-insensitive to rate-sensitive entropy production is observed, which is attributed to the breakdown of the relation Y/G ≈ const. In particular, the entropy is partitioned into configurational and vibrational entropies. Thus, the framework can provide a thermodynamic variable to describe microstructural change. The configurational entropy is found to be proportional to the Shannon entropy form describing structural change. Overall, our framework provides a comprehensive understanding of compression, plastic deformation, structural change, entropy production, and their relationships in QIC. Furthermore, it greatly extends the capability of atomistic simulations to interpret QIC experiments.

Acknowledgment. We are grateful for insightful discussions with Professor Yang Yong at the Department of Mechanical Engineering, City University of Hong Kong. This work was supported by the NSAF under Grant No. U1830206, the National Key R&D Program of China under Grant No. 2017YFA0403200, the National Natural Science Foundation of China under Grant Nos. 11874424 and 12104507, and the Science and Technology Innovation Program of Hunan Province under Grant No. 2021RC4026.