Excitations follow (or lead?) density scaling in propylene carbonate Free

For many decades, there has been evidence that discrete molecular rearrangements, so-called “hops,” contribute significantly to relaxation in viscous liquids.^1,2 At reduced temperatures, it is thought that molecules become trapped in a cage formed by their neighbors so that relaxation requires collective rearrangements of neighborhood particles³ that are thought to be facilitated by structural excitations.^4,5

The case for the importance of excitations in influencing liquid dynamics is supported through their ability to account for super-Arrhenius relaxation,⁶ the breakdown of the Stokes–Einstein relation,⁷ the appearance of dynamic heterogeneity,^8,9 and the behavior of specific heat at the calorimetric glass transition.¹⁰ If excitations are a causal factor in liquid dynamics, one would also expect a connection to the ubiquitous scale-invariance observed in liquid dynamics,^11–25 but this connection has not been explored.

As pointed out by Roland et al.,¹⁵ expectation of scale-invariance arises when the potential that controls system behavior is well-approximated by a pairwise additive inverse power law (IPL).^16,19 Dyre and co-workers^17,19,20,26 have expanded these ideas, demonstrating pressure and temperature invariance of dynamic, structural, and thermodynamic properties of model and real liquids when length and energy are properly scaled.²⁰ These properties are invariant along “isomorphs” of constant ρ^γ/T, where γ = n/3 and n is the exponent of the IPL pair potential. In systems where only the IPL potential is important, the isomorph scaling exponent γ is related to instantaneous fluctuations in energy and the first pressure virial.^17–20 

When the dominating interactions in a liquid do not conform to IPL potentials, the entire gamut of isomorph behavior can be lost. This is observed in systems with multiple interaction length scales, dominant electrostatic forces, or those with distinct directionality such as hydrogen bonding.¹⁹ For example, molecular dynamics (MD) simulations of the SPC water model have shown that there is essentially no virial–energy correlation and that isomorphic state points do not exist.^20,26 Similar results are found in experimental studies of heterogeneous hydrogen bonding (e.g., polyalcohol) systems.^27,28

When only minor motional modes of a system (such as intramolecular or long-range intermolecular degrees of freedom) do not conform to IPL potentials, some aspects of isomorph behavior are lost, but many other aspects may be retained. These systems exhibit what is sometimes referred to as “pseudoisomorphs,”²⁹ and properties associated with long timescales such as α relaxation may exhibit density scaling in that they are invariant at constant values of $ργTS/T$ where γ_TS is an empirical density scaling exponent. Density scaling has been experimentally verified for a wide variety of systems, including van der Waals, and other nonassociating liquids, polymers, ionic liquids, and some associating and hydrogen bonding liquids.^{12–15,21–25} Several theoretical and experimental works^19,30,31 have demonstrated a connection between γ_TS and the volume dependence of vibrational frequency for spatially extended collective modes of the liquid through the Gr$ü$neisen parameter.

In these pseudoisomorph systems, pressure–energy correlations are broken since the energy can partition to the non-IPL modes without strong coupling to changes in volume or pressure. The pressure of these fluids does not depend strongly on forces arising from angle vibrations and torsional rotations of at least some intramolecular modes.³² Olsen et al.³³ found in model diatomic and oligomeric systems with harmonic intramolecular bonds that vibrational eigenmodes representing center-of-mass (COM) motion scale properly and could be used to find the scaling exponent but that the highest frequency eigenmodes did not scale. The number of non-scaling eigenmodes was equal to the number of intramolecular bonds, so all of the intramolecular modes for the diatomic model and at least 1/3 of the highest frequency intramolecular modes for the oligomeric system did not scale.

Thus, while low-frequency vibrational modes and slow dynamic processes follow density scaling, there appears to be some cutoff in timescale or length scale below which dynamic processes in the liquid no longer scale. T$ö$lle et al. found that relaxation processes in ortho-terphenyl on timescales of 1 ns and longer follow density scaling.¹¹ Hansen et al.²³ and Ribeiro et al.³⁴ demonstrated invariance down to 100 and 10 ps, respectively. Puosi et al.²⁴ found density scaling at timescales of the “fast β” (β_fast) relaxation, roughly 1 ps. In that work, Puosi et al. quantified the mean-squared displacement $⟨u2⟩$ of particles in course-grained polymer systems, finding excellent correspondence between γ_TS obtained for $⟨u2⟩$ and α relaxation, but only when a small term quadratic in $(Tρ−γTS)$ was added to fit the latter.

We have recently shown that the Angstrom-length scale and picosecond timescale motion that has traditionally been classified as β_fast relaxation and measured through $⟨u2⟩$ contains signatures of two separate processes.^35,36 The faster of the two corresponds to elastic deformations of equilibrium local atomic structure, corresponding to so-called “inherent state” (IS) dynamics on a potential energy landscape (PEL).³⁷ The slower of the two processes represents^36,38 rearrangements of local atomic structure corresponding to interbasin (IB) barrier crossings on a PEL that lead to a new IS and the instantaneous population of IB crossing events is equivalent to the excitation population.³⁹ We have further found evidence that intramolecular modes may figure prominently in excitations,³⁸ leading us to question whether excitations and all the components of β_fast should be expected to have the same γ_TS as the α process or whether involvement of intramolecular modes in these short-timescale, local motions breaks this scaling.

To investigate to what extent excitations obey density scaling, we have conducted molecular dynamics simulations of propylene carbonate (PC) over a range of pressures and temperatures. PC was previously shown to follow density scaling^31,40 and provides an interesting case study for this investigation due to its strong electrostatic interactions (ɛ ∼ 65), molecular anisotropy, and the considerable flexibility of its ring, all of which affect its liquid state properties.^41–43 

Our simulations utilize an ab initio, polarizable force field capable of predicting properties that are in excellent agreement with experimental measurements near STP. By searching its phase diagram, we identify two P–T curves in PC that exhibit largely invariant structural and dynamical properties. We then quantify excitation populations at these phase points, utilizing a previously benchmarked approach.^35,36 We conclude that excitations fundamentally involve the atomic-scale intramolecular structure of PC and are not well characterized by center-of-mass motion alone. Furthermore, we demonstrate that excitation populations in this liquid obey density scaling as they are constant across the P–T curves exhibiting invariant structure and dynamics. We also show that the motions associated with excitation are the fastest dynamic processes that can be considered to follow density scaling.

Molecular dynamics (MD) simulations were conducted for propylene carbonate (PC) utilizing the OpenMM simulation package.⁴⁴ We utilize the atomistic, polarizable symmetry adapted perturbation theory-force field (SAPT-FF), which is an ab initio force field developed on the basis of symmetry adapted perturbation theory (SAPT) calculations.⁴⁵ We supplement these existing force field parameters with atomic charges and dihedral potentials that are specifically parameterized for PC utilizing additional density functional theory (DFT) calculations as described in the supplementary material. Polarization is incorporated via a Drude oscillator model utilizing a dual-Langevin thermostat as implemented in OpenMM, with a 1 ps⁻¹ friction coefficient used for both thermostats.⁴⁶ The particle–mesh Ewald method (PME)⁴⁷ is used for electrostatics and van der Waals (VDWs) interactions are computed up to a 1.4 nm cutoff. Each simulation consisted of 400 propylene carbonate molecules initially constructed using Packmol.⁴⁸ Equilibration consisted of 1 ns NPT simulations conducted with a Monte Carlo Barostat, followed by subsequent 5 ns NVT simulations at the equilibrated density. All simulations utilized a 1 fs time step. Because OpenMM utilizes a Monte Carlo Barostat, the virial is not directly accessible. We have thus utilized Gromacs to conduct an additional MD simulation in the NVT ensemble at 300 K to compute energy–virial correlation, shown in the supplementary material, Fig. S1.⁴⁹ The simulation details for the Gromacs simulation are largely identical to the OpenMM simulations, except that in the former Drude oscillator positions are treated strictly adiabatically (minimized at every timestep). See the supplementary material for more details about the simulation settings.

Finding isomorph families for Lennard-Jones (LJ) fluids and other simple systems for which pairwise IPL potentials dominate is accomplished by analyzing the virial–energy correlation function.²⁰ For IPL potentials, the virial and potential energy are perfectly correlated, so that the normalized correlation function is unity.²⁶ This is not possible for systems with significant intramolecular degrees of freedom. For organic liquids with molecular constituents of increasing size and complexity, it is expected that isomorphic behavior will not be found based on analysis of virial–energy correlation.^23,24,29 This is because intramolecular rather than intermolecular interactions dominate the potential energy of the system, with the former growing as 3N − 6 as the size “N” of the molecules gets larger. The intramolecular interactions do not follow inverse power laws; additionally, fluctuations in the intramolecular energy are not expected to correlate with virial fluctuations. The latter is inferred from the fact that the pressure of fluids does not depend on forces arising from angle vibrations and torsional rotations.³² For molecules of large size and complexity, near zero correlation of the virial and energy is expected due to the dominance of the intramolecular interactions. Experimentally, Hansen et al. have shown that density scaling breaks down for dynamic signatures indicative of intramolecular motion.²³ 

We have computed the virial–energy correlation from an MD simulation of PC at 300 K, 1 atm; the results are shown in the supplementary material, Fig. S1. We plot the data as normalized fluctuations of each quantity, in analogy to previous work.^24,26,29 As expected from the dominance of the intramolecular energetics of PC and similar to results for other flexible molecules,^24,29 we find essentially no correlation between the virial and energy fluctuations.

We note that if it were possible to explicitly separate the intramolecular and intermolecular energy contributions, then one could analyze virial correlations with just the latter. However, such a separation is generally not tractable within a molecular dynamics code. The reason being that electrostatic interactions, which contribute to intramolecular energy for 1–4 and longer atom separations, are computed with an Ewald/PME sum and not easily separated on a pairwise interaction basis. Furthermore, the electrostatic energy is formally many-body for explicitly polarizable force fields, as utilized in this work.

Because of these considerations, in this work we utilize a brute-force approach to search for isomorph curves. We searched for two isomorph curves in the P–T phase diagram, shown in Fig. 1. We initially ran a 1 ns NPT simulation of 400 PC molecules at temperature/pressure values of 300 K/1 atm; simulation settings are the same as noted above. We then performed a 20 ns NVT simulation. We calculated the diffusion coefficient (D_T) and the inverse rotational correlation time (1/τ_rot) for this state point. D_T was calculated using the Einstein relation given by

and calculated τ_rot by integrating the rotational correlation function as follows:

with t_max equal to 1 ns and z_CO equal to the carbonyl bond vector. As isomorph properties are only apparent in reduced units, the calculated values of D_T and τ_rot were converted to reduced unit values $D̃T$ and $τ̃rot$⁠.²⁰ We performed a series of 1 ns NPT/20 ns NVT simulations at temperatures of 325, 350, 375, and 400 K and a range of pressures. $D̃T$ and $1/τ̃rot$ were computed at each state point. A point was classified as an isomorph point if the calculated values agreed with $D̃T$ and $1/τ̃rot$ at 300 K/1 atm (within statistical uncertainty).

We followed a very similar procedure for the second isomorph. Here, we initially ran a simulation at a temperature/pressure value of 250 K/1 atm. As the dynamics are noticeably slower at this state point compared to 300 K/1 atm, we ran 50 ns NVT simulations after the initial 1 ns NPT simulation; both $D̃T$ and $1/τ̃rot$ were then computed at this P–T point from the NVT simulation. We then followed a similar procedure as above, running simulations at varying pressures at the temperature points 275, 300, 325, and 350 K. The temperature/pressure points making up both isomorphs along with their calculated values of $D̃T$ and $1/τ̃rot$ are listed in Table I.

Isomorph and density scaling behavior present in reduced units of time, length, and energy as described by Gnan et al.²⁰ The reduced units of length and time used here are ρ^−1/3 and $m/kBT/ρ1/3$⁠, respectively. Here, we take the mass “m” as the molecular mass of PC, but note that most scaling arguments are derived for atomic fluids. We compute scaled diffusion coefficients as $D̃=ρ1/3m/kBTD$⁠. We compute scaled time constants τ using the reduced time unit. Dynamic scattering functions are given in real length/time units when comparing to experiment, and reduced length/time units otherwise, and this convention will be explicitly stated in context of discussion.

We note that experimental characterization of density scaling has primarily utilized real rather than reduced units. Formally, the latter should always be used, but in practice use of real units is largely inconsequential when relaxation times are plotted over many orders of magnitude. For example, Pawlus et al.⁴⁰ demonstrated density scaling for propylene carbonate by characterizing relaxation times over ten orders of magnitude, at different phase points. While relaxation times were plotted in real units, the conversion to reduced units differs by order unity for the investigated phase space and is thus not observable over plots that span many decades in time. We show a plot comparing reduced vs real units for D_T vs τ in Fig. S2 for the PC pseudoisomorph points.

In computer simulations, excitation populations are generally quantified by searching for particles that undergo sufficiently large spatial displacements that persist in time, usually on the order of a picosecond.^5,50 In such an analysis, one must average out fast vibrations that correspond to motion within the same inherent structure, and one must justify the spatial and temporal cutoff values applied to discriminate hops.^5,50 Instead, we quantify excitation populations through analysis of the self-intermediate scattering function F(q, t)^35,36 by modeling it as a double-Gaussian in q as follows:

This approach has the advantage of being applicable to both experimental and computational characterization, and it requires no a priori assumptions about length scales or timescales of the excitations since the length scale and timescale signatures of excitations and their associated hops come directly from fitting to the model. Here, Φ(t) signifies the fraction of particles (H atoms or molecular centers of mass) that have undergone an excursion of length scale σ_IB up to time t. Φ(t ≈ 1 ps) is the instantaneous population of particles at an IB barrier and thus the instantaneous excitation population. Moreover, at t ≈ 1 ps, σ_IS and σ_IB are the characteristic length scales for inherent state elastic deformation and for hopping, respectively. In Eq. (3), we have neglected a diffusive term, so, at later times, σ_IS and σ_IB begin to grow slowly as they acquire diffusive character. This form provides excellent fits for experimental and simulation data as is shown in Figs. 3–5. The time-dependence of fit parameters σ_IS, σ_IB, and Φ is shown.

We first briefly discuss the accuracy of the PC force field utilized for MD simulations as this dictates the accuracy of all the simulation results presented in the manuscript. Compared to experiment, errors in the predicted density and enthalpy of vaporization are 3% and 9%, respectively, at 300 K, 1 atm; such accuracy is comparable to that of previous property predictions of organic liquids with the SAPT-FF force field.⁵¹ Additionally, the computed diffusion coefficients and rotational correlation times are within 10% of the experimental values over the temperature range of 250–350 K at 1.0 atm (Fig. S8). Because there are no empirical parameters in the model, similar accuracy is expected throughout the studied temperature–pressure phase diagram.

Utilizing a brute-force screening approach (Sec. II B), we were able to locate two curves on the temperature–pressure phase diagram that display largely invariant structural and dynamical properties (computed in reduced units). The curves that we located are shown in Fig. 1, in which we have mapped the pressure to the equilibrium density. Also shown in Fig. 1, for perspective, is the glass transition curve.⁵² One pseudoisomorphic curve starts at 300 K, 1 atm and extends to higher temperature/pressure, and the other starts at 250 K, 1 atm; by intention, both curves were chosen to be far away from the glass transition.

Our isomorphic curves display identical density scaling to that previously determined experimentally for propylene carbonate. Pawlus et al.⁴⁰ found that relaxation times and dc conductivity in propylene carbonate solutions measured across the phase diagram collapsed to a common curve when plotted against V^−3.7/T.⁴⁰ This scaling exponent of γ_TS = 3.7 agrees with our computed isomorph curves as shown in Fig. 1. All points on a given isomorph curve were fit to the form $ργTS/T=constant$⁠. A best fit of our data gave a value of γ_TS = 3.75 ± 0.09, which is consistent with the findings of Pawlus.

In Table I, we give calculated diffusion coefficients $(D̃T)$ and inverse rotational relaxation times $(1/τ̃rot)$⁠, both in reduced units, for the two pseudoisomorphic curves. These computed dynamical properties are invariant along each pseudoisomorph to within statistical uncertainty. In addition to invariance to these dynamic measures, we expect structural properties of the fluid to be invariant along pseudoisomorphs. To analyze the liquid structure, we compute the radial distribution functions (RDFs) between propylene carbonate molecules in the liquid as a function of center-of-mass positions. The RDFs are shown in Fig. 2(a) (top panel) and in Fig. S3 of the supplementary material as computed for all phase points on each pseudoisomorph curve. The RDFs for all points on each pseudoisomorph collapse to a single curve when the length scale is plotted in reduced units. In the inset, we highlight the similarity across phase points with a high resolution scale of the first RDF peak. For comparison, we show RDFs computed from 1 atm isobar phase points spanning the temperature range 300–350 K in Fig. 2(b). It is clear from the insets in Figs. 2(a) and 2(b) that there is structural deviation (albeit small) between the isobar points, in contrast with the identical structure among pseudoisomorph phase points.

In the bottom panel of Fig. 2, we show the self part of the intermediate scattering function F(q, t), defined as

where q is the momentum transfer wavevector and ⟨⋯⟩ denotes an ensemble average. In Fig. 2(c), we show F_s(q, t) computed for molecular center-of-mass (COM) positions, and in Fig. S3 we show F_s(q, t) calculated for hydrogen atom (H-atom) positions [i.e., the sum in Eq. (4) runs over COM or H-atom indices, respectively]. Here and for the remainder of the manuscript, all reference to the scattering function concerns the self part (and not the distinct contribution), and the quantity in Eq. (4) will be referred to simply as $Fq,t$ (“s” subscript dropped). The q values shown describe dynamics on length scales ranging from 70% to 2.5 times the average intermolecular spacing (q_max,COM = 5.76). At each of these length scales, COM dynamics among all the phase points investigated are identical to within uncertainty. These combined dynamic and structural analyses indicate that our identified phase point curves (Fig. 1) indeed behave as pseudoisomorphs. For comparison, in Fig. 2(d), we show F_s(q, t) computed from 1 atm isobar phase points spanning the temperature range 300–350 K (computed for COM positions) As expected, there is significant temperature dependence of dynamic relaxation for the isobar points; in contrast, the pseudoisomorph phase points exhibit essentially identical COM dynamics as indicated by F_s(q, t) in Fig. 2(c).

Based on our simulation analysis, we find that propylene carbonate exhibits density scaling such that structural and dynamic properties are functions of ρ^3.7/T, in agreement with previous experimental findings.⁴⁰ 

Interbasin barrier crossing events in atomic model systems are local cooperative particle reorganizations^3,53 and represent fundamental relaxation events. We have recently shown that motion associated with IB crossing and IS basin exploration can be quantified from the scattering function.³⁶ Specifically, we showed that experimental S(q, ω) [the time Fourier transform of F(q, t)] from quasi-elastic neutron scattering (QENS) on propylene carbonate is optimally fitted with three Lorentzians, having q-dependent amplitudes consistent with a diffusive process and two localized processes. The latter two processes had time-dependent length scale characteristics corresponding to localized processes, consistent with inherent state (IS) vibrational motion (σ_v) and interbasin (IB) hopping motion (σ_h) on a PEL. In order to characterize IB and IS motions, in this study we focus on F(q, t) calculated from Eq. (4). For comparison with QENS data from PC that were previously acquired by one of us,³⁶ we sum over H atoms because incoherent scattering from hydrogen atoms dominates the QENS signal in simple organics.

In Fig. 3, we compare the scattering function computed from our simulations to that previously measured by QENS.³⁶ There is generally good agreement between the simulated and experimental scattering functions, attesting to the accuracy of the employed force field. Since we do not include coherent scattering in Eq. (4), we expect deviations of a few percent at low q values, which are not visible in Fig. 3. However, there are prominent discrepancies between simulated and experimental $Fq,t$ at q > q_max = 1.1 Å⁻² and at short times and high q values.

These discrepancies are not surprising, considering the approximations of the intramolecular interaction potential (harmonic bonds, angles, etc.) used in the simulations (supplementary material). For example, the largest plotted wavevector q² = 5 Å⁻² corresponds to less than 3 Å real space distance, which is smaller than the molecular diameter of PC. Thus, $Fq,t$ at this wavevector is expected to be quite sensitive to the torsional potentials and even methyl group rotation, the latter of which is not expected to be accurately modeled by the force field. In contrast, the scattering at q < q_max (distances longer than intermolecular separation) is dictated largely by intermolecular forces and the intermolecular potential utilized in the MD simulations is expected to be quite accurate, as confirmed by the good agreement with experiment for $Fq,t$ in this regime.

We obtain parameters for excitations, IB and IS motion by fitting F(q, t) with Eq. (3). A straight line in a ln(F(q, t)) vs q² plot signifies a single relaxation process that is Gaussian in q. The bilinear form of $lnF(q,t)$ vs q² illustrates the two-step relaxation process. From a data-fitting perspective, Φ(t) serves as a time-dependent non-Gaussian parameter. From a particle dynamics perspective, Φ(t) represents the fraction of particles that have “hopped” (undergone a non-reversing IB crossing) up to time t. Bearing in mind that hops constitute significantly larger excursions than IS cage distortions, it is clear that until a molecule executes a hop, the q-dependence of its scattering signature will be characteristic of only small length scale σ_IS. Once a hop occurs, that signature will change to larger length scale for all subsequent times.

Figure 4 shows F(q, t) calculated from hydrogen atoms, fits to Eq. (3), and fit parameters. As shown in the top panel, the model fits quite well, fully parameterizing F(q, t) for this q and t range. The middle and bottom panels in Fig. 4 show excellent agreement in the time-dependence of model fits for each of the 300 K, 1 atm pseudoisomorphs. In the middle panel, Φ(t) increases monotonically as expected,³⁶ with (1 − Φ) exhibiting an initial rapid drop during the period of quasi-ballistic motion when atoms are exploring the cage formed by their neighbors. As we have explained elsewhere,³⁹ just after this initial drop in (1 − Φ), Φ(t ≈ 1 ps) represents the instantaneous fraction of H atoms involved in large-length scale motion characteristic of structural excitations. After this initial period of cage exploration, (1 − Φ) decays almost exponentially as excitations visit new regions of the sample and facilitate hops of an increasingly large fraction of particles. In the bottom panel, the characteristic length scales σ_IS and σ_IB initially show signatures of localized motion during the period of cage exploration, reaching quasi-asymptotic values of ≈0.1 and ≈0.2 respectively, consistent with IS vibration and IB crossing.⁵³ The values of these fit parameters continue to grow, taking on a diffusive character at longer times in this model because, for simplicity’s sake, we do not explicitly include diffusion through a separate parameter.

The high degree of overlap among the fit parameters at each of the pseudoisomorphic phase points signifies excellent correspondence between F(q, t) for H atoms at each phase point and demonstrates that excitation populations as well as dynamic timescales and length scales associated with IS vibrations and IB crossing obey density scaling. We have previously shown that the IS and IB modes are related in a simple way to α relaxation and translational diffusion,^35,36,38 so the fact that they exhibit scaling is consistent with Hansen’s conclusion (based on longer timescale data) that motion coupled to α relaxation should exhibit density scaling.²³ 

While the dynamics encoded by the diffusion coefficient, rotational relaxation time (Table I), and F(q, t) for COM integrate over intramolecular motion, the atomic scattering function formally encompasses dynamical processes over all temporal and spatial scales relevant to IB and IS motion. The fast timescales and short length scales of IS and IB motion make them candidates for being associated with intramolecular modes. Comparing the top panels of Figs. 4 and 5 provides evidence that IS and IB modes in PC involve significant intramolecular motion. The top panel of Fig. 5 shows F_COM(q, t) calculated from molecular centers of mass for the 300 K, 1 atm phase point. Comparing with F_H(q, t) from H atoms at the same phase point in Fig. 4, it is quite clear that the COM motion is of lower amplitude overall since, for all q and t, F_H(q, t) ≤ F_COM(q, t). A similar relationship was found between COM and H-atom F(q, t) for glycerol.³⁸ Molecular rotation will cause some differential decay of F_H(q, t) compared to its COM counterpart, but a detailed comparison of the two scattering functions suggests another significant contribution. The two characteristic length scales for motion are similar for H-atom and COM fits, suggesting that the molecule moves in a quasi-rigid fashion, but Φ_H ≥ Φ_COM at all times. This pronounced muting of the non-Gaussian nature of the COM F(q, t) with respect to the H-atom F(q, t) directly demonstrates that IS and IB motion executed by PC are largely intramolecular in nature (although they could in principle involve parts of more than one molecule). Hopping of individual H atoms—parts of the molecule—precedes similar length scale motion of the entire molecule reflected in the COM.

We turn now to the question of which intramolecular modes are likely to contribute to IB and IS motions and why, if they contribute significantly, do they not spoil the density scaling demonstrated in Fig. 4. The origin of the effective IPL appears to be collective in nature since the strong energy–virial correlations are observed only when averaging over the ensemble in rigid particle models.¹⁹ Thus, we expect that if any intramolecular mode does contribute to density scaling, the contributing modes would most likely be low-frequency ones that involve significant portions of the molecule.

In the top left panel of Fig. 6, we overlay Raman spectra of PC from two sources. The red solid trace shows the lowest frequency intramolecular vibrational modes, obtained through spontaneous Raman scattering. These Raman modes at 191, 318, and 451 cm⁻¹ are all molecularly extended ring deformation modes. The black trace is the Fourier transform of the β_fast component from time-domain optical polarizability measured by Optical Kerr effect (OKE) spectroscopy.⁵⁴ The β_fast OKE signal contains a collective and librational part, which we demonstrated to have the same temperature-dependent frequencies and relative amplitudes as that of IB barrier crossing and IS motion measured by neutron scattering.⁵⁴ The blue dashed line represents the spectrum of IB motion and the red dashed line represents the IS spectrum. The black line is the sum of the IS and IB spectra. In the PEL picture, both the IS and IB modes are coupled to relaxation in the sense that the IS modes assist in basin exploration but the molecular reorganizations that make up relaxation events proceed directly along IB mode coordinates. The intramolecular vibrational modes with frequencies below 500 cm⁻¹ overlap spectrally with and thus may couple to the IS and IB modes.

In the right top and bottom panels of Fig. 6, we overlay power spectra of velocity autocorrelation functions for ring atoms and molecular centers of mass, respectively, at the three simulated pressures indicated. The right top panel shows that the pressure dependence of the intermolecular mode characteristic frequencies drops with increasing frequency, indicating that higher frequency peaks contribute less to density scaling. Consistent with this, we found that eliminating the much higher frequency bond-stretching modes by constraining all bond lengths led to no change in density scaling (supplementary material, Fig. S4). Also consistent with these findings, Raman modes found below 800 cm⁻¹ in organic molecules with more than one backbone atom typically involve molecular deformations and typically have larger Gr$ü$neisen parameters than do higher frequency modes involving bond stretches and H-X bending.^55–57 

The broad β_fast features below ≈150 cm⁻¹ in the right bottom panel are due to rotational and translational motion, but they arise only from translational motion in the right bottom panel. In this spectral region, the IS and IB responses are not well-separated, but their characteristic frequencies appear to have relatively strong pressure dependencies. We quantify the characteristic frequencies of the IB and IS motion by fitting the COM power spectra to functions commonly used to characterize OKE data⁵⁴ as described in the supplementary material. We find acceptable fits using functions corresponding to IB and IS motion and an additional Gaussian process at slightly higher frequency than the IS mode. In the left bottom panel of Fig. 6, we show the model fit to the 300 K, 1 atm COM power spectrum. Analogous fits for 1900 and 4600 atm are shown in supplementary material Fig. S5 and are of similar quality.

In the inset to the right bottom panel of Fig. 6, we show Gr$ü$neisen parameters for each of the modes found in the power spectra, calculated as

where i represents the ith vibrational mode. Among the modes characterized, we see a continuous progression of increasing γ_G with decreasing frequency. We emphasize that, although the IB and IS modes are distinct from intramolecular eigenmodes, the differences in F_H(q, t) and F_COM(q, t) suggest they involve intramolecular motion and thus likely couple to these modes.

Finally, we note that γ_IB ≈ γ_TS. This suggests that it may be generally possible to find γ_TS in a simulated liquid by simply finding γ_IB. It also suggests that the relationship between γ_TS and the Gr$ü$neisen parameter may be more direct than previously thought. Associations between γ_TS and γ_G have been made theoretically through the latter’s connection to configurational entropy³⁰ or through the kinetic energy derivative of the pressure in the harmonic limit,¹⁹ but the striking result in the inset to Fig. 6 suggests that the Gr$ü$neisen parameter may influence γ_TS directly through γ_IB.

The Gr$ü$neisen parameter formally characterizes the volume dependence of normal mode frequencies in a crystal, but, of course, there is no relaxation in an ideal crystal. Similarly, Debye–Waller (DW) factors strictly describe only harmonic vibrations but are used to describe liquid relaxation through the mean square displacement $⟨u2⟩$⁠. An extensive literature, going back three decades^35,58–61 has explored the relationship between α relaxation and $⟨u2⟩$ measured at ≈1 ps, the time of the β_fast process. In these approaches, it is assumed either that the amplitude of harmonic motion describes a localization volume for particles trapped in a cage of their neighbors or is proportional to the probability of a local relaxation event. Recently, we have made a more direct connection between motions that constitute $⟨u2⟩$ in amorphous systems and relaxation processes by demonstrating that $⟨u2⟩$ contains signatures of both IS basin exploration and elemental relaxation events—the IB barrier crossing. Thus, when properly parsed, dynamic signatures of the elemental relaxation process are obtained directly from the β_fast relaxation of $⟨u2⟩$ or F(q, t).^35,36

In analogy with the idea that $⟨u2⟩$ is directly related to relaxation through its reporting on IB barrier crossings, we propose that a causal connection exists between γ_G and γ_TS through γ_IB as follows: The normal modes of the liquid that couple to the IB process are by definition found at the Frenkel frequency,⁶² ν_F = 1/τ. At this frequency, the liquid supports longitudinal modes, but transverse modes are overdamped since the liquid can reorganize on the timescale of these local shear deformations, and this motion directly couples to relaxation. In an amorphous system, longitudinal and transverse modes will also be coupled since there is no well-defined extended lattice. Since the motion along these mode coordinates is coupled to elemental relaxation (IB crossing) events, the resonant frequency of the mode can be thought of as a relaxation attempt frequency. Thus, for a fixed barrier height, small changes in the frequency of these modes would lead to directly proportional changes in the rate of elemental relaxation events. We have previously demonstrated that for T > T_B (T_B being where the α and Johari–Goldstein β relaxations bifurcate), the timescales for IB crossing events and α relaxation have the same temperature dependence.³⁶ Thus, for this temperature range, the Gr$ü$neisen parameter of the Frenkel modes (γ_IB), which describes the pressure dependence of their frequency, also describes the pressure dependence of the α relaxation time (γ_TS) through the pressure dependence of the elemental relaxation time. As we see in Fig. 6, γ_IB ≈ γ_G.

Above, we showed that high-frequency modes associated strictly with bond stretching have no impact on density scaling. Beyond arguing that the pressure dependence of lower frequency intramolecular modes and their spectral overlap with IB modes suggest that the former may contribute to IB barrier crossings and to density scaling, we can only speculate as to how important they may be. As we have shown above, hops associated with IB transitions appear to occur at the atomic rather than molecular level. Thus, intramolecular modes that conserve molecular center of mass may well contribute to hops of constituent atoms. However, to quantify this would require a series of simulations in which we systematically constrain these various intramolecular modes. While standard MD simulations are done using Cartesian equations of motion, an alternative approach for exploring the role of intramolecular modes on liquid dynamics would be to perform molecular dynamics in general/internal coordinates (or hybrid) to enable explicit interrogation of specific internal modes.^63–65 This direction will be explored in the future work.

While evidence of density scaling in β_fast at ∼1 ps timescales has been demonstrated previously, the relationship involved a quadratic correction term.²⁴ Given that β_fast relaxation has contributions from distinct modes of motion, some of them potentially intramolecular, and given that intramolecular vibrational modes can destroy the virial–energy correlation that is often the formal basis for density scaling arguments, it was unclear whether all aspects of β_fast motion would couple or whether the necessity of a quadratic correction term signaled a partial decoupling. In this work, we have shown that the IB barrier crossing aspect of the β_fast, a fundamental structural excitation mode, seems to obey the same density scaling, whereas the higher frequency IS mode is associated with a nearly twofold smaller scaling factor. Because structural excitations correspond to fundamental hopping (intrabasin barrier crossing) processes on the liquid potential energy landscape (PEL), our results suggest a novel microscopic perspective for density scaling. Interpretation of density scaling based on excitations may prove a promising general perspective applicable to molecular liquids for which virial–energy correlation or single inverse power law arguments may not hold due to the intramolecular modes or due to multiple effective potentials.

For propylene carbonate, we have demonstrated that the characteristic length scales σ_IS and σ_IB have the same values in F(q, t) computed for COM and to H-atom motion but that Φ_(t) calculated for COM always trails that of H atoms. We interpret this as demonstrating that excitations must thus be interpreted as essentially atomic rather than molecular displacement events. This important conclusion has implications for generalizing excitation theories of atomic and model fluids^5,50 to real molecular liquids. The fact that the Gr$ü$neisen parameter for collective IB modes is approximately equal to γ_TS suggests that finding these modes in molecular liquids will be the key to connecting excitation approaches in these liquids to those of atomic model systems.

The supplementary material includes figures and brief discussion concerning the following: the lack of virial–energy correlation for PC, a brief comparison of the use of real and reduced units for finding pseudoisomorphs in PC, a description of the radial distribution functions and intermediate scattering functions for pseudoisomorph families, the pressure dependence of bond-stretching vibrations and the fitting of velocity autocorrelation power spectra to find the pressure dependence of low-frequency intermolecular collective motion and intramolecular molecule deformation vibrations, and, finally, force field construction and benchmarks.

J.P.S. and J.G.M. acknowledge support from the Air Force Office of Scientific Research under Award No. FA9550-22-1-0025.

The authors have no conflicts to disclose.

John P. Stoppleman: Data curation (equal); Formal analysis (equal); Writing – original draft (equal). Jesse G. McDaniel: Funding acquisition (equal); Project administration (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). Marcus T. Cicerone: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Project administration (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.