Mesoscopic two-point collective dynamics of glass-forming liquids

Glass-forming (GF) liquids exhibit dramatic viscous slowing down upon cooling, departing from the familiar Arrhenius temperature dependence.^1–5 This work is concerned with a simple and yet incompletely understood question about this seemingly enigmatic phenomenon: how do the collective motions of liquids^6–9—as characterized by the two-point density correlation functions, such as the dynamic structure factor S(Q, ω) and intermediate scattering function S(Q, t)—change during the liquid-to-glass transition, particularly around the crossover temperature where activated, cooperative dynamics start to emerge?^10–13 While S(Q, ω) and S(Q, t), in principle, can be probed by inelastic neutron scattering experiments, measurements of coherent dynamics of GF liquids in the relevant spatiotemporal regime are still technically challenging, and the existing data^14–20 are insufficient for fully addressing the aforementioned question. On the other hand, the extensive molecular simulations^21–27 during the past several decades have not led to a clear picture either partially due to the lack of a proper angle to approach the data analysis. Additionally, the majority of the experimental and computational studies of GF liquids have so far focused on analyzing the dynamics at relatively large wave numbers and overlooked the collective motions on mesoscopic scales.

From a theoretical perspective, cooperativity is a defining characteristic of molecular motions during the glass transition.^{10–13,28–33} It is widely believed that the cooperativity in GF liquids is best captured by high-order dynamical correlations^33–42 and four-point correlators,^{33,36–39,41,42} in particular. In this context, we point out that a number of scattering experiments and computer simulations have revealed qualitative changes in two-point dynamical correlations in GF liquids.^43–49 Some argue from the results of polarized quasielastic light scattering experiments that the dynamic crossover is accompanied by the appearance of slow long-range density fluctuations.^43,44,46 Nevertheless, because of the scarcity of systematic data, many issues regarding the fundamental features of two-body collective dynamics during the liquid-to-glass transition remain unresolved.

To establish the basic phenomenology of collective dynamics of GF liquids in terms of the classical two-point correlation functions, we spatiotemporally map out the normalized coherent intermediate scattering function S(Q, t)/S(Q) and hydrostatic pressure correlation function C_pp(Q, t)/C_pp(Q) of atomic as well as polymeric glass-formers over a large swath of wave numbers and correlation times using molecular dynamics simulations. For reasons that shall become apparent in the following, S(Q, t)/S(Q) is the most appropriate choice of two-point density correlation functions for the current purpose as opposed to the dynamic structure factor or van Hove function.⁵⁰ To overcome the shortcoming associated with the traditional practice of analyzing intermediate scattering functions at discrete wave numbers, we adopt a “two-dimensional” approach⁵¹ by directly visualizing S(Q, t)/S(Q) and C_pp(Q, t)/C_pp(Q) on a dense grid of correlation times and wave numbers. Our calculations encompass both small length scales at large wave numbers and mesoscopic scales beyond the first structural peak, where long-range collective density and pressure fluctuations take place. Using the linear hydrodynamic theory^9,52–58 as a reference, we show that the dynamic crossover below the Arrhenius temperature T_A is characterized by the emergence of collective, “nonhydrodynamic” density–density and hydrostatic pressure–pressure correlations at long wavelengths. The observed nonhydrodynamic mode is independent of wave numbers and closely coupled to the local dynamics at the first peak of the static structure factor S(Q). These findings point to a qualitative, mechanistic change in mesoscopic collective density and pressure fluctuations at the Arrhenius temperature and provide a new angle from the perspective of two-point collective dynamics to approach the dynamic crossover phenomenon.

This paper is organized as follows: To put our current investigation in perspective, we first present a brief overview of the relevant existing studies concerning mesoscopic two-point collective liquid dynamics in Sec. II in the context of the glass transition phenomenon. In Sect. III, we outline the general strategy of the current work for analyzing the collective dynamics of glass-forming liquids. The technical details of our molecular dynamics simulations and hydrodynamic theory calculations are given in Secs. IV and V. Section VI describes our observations from the simulations with regard to the mesoscopic two-point collective dynamics. Several implications of these results are discussed in Sec. VII.

While the origin of the glass transition phenomenon is still a subject under intense debate, it is generally believed at least according to the traditional views^10–13 that the dramatic viscous slowing down upon cooling involves some kind of cooperative, local structural rearrangements. As a result, significant attention has been directed to understanding the dynamics of GF liquids at such small length scales, which in the language of two-point density correlations are broadly defined by the “first” structural peak. On the other hand, glass transition is also a macroscopic phenomenon in the sense that the sharp change in structural relaxation time can be readily probed by bulk techniques, such as rheometry and dielectric spectroscopy. However, the collective dynamics of GF liquids on mesoscopic scales—length scales that are larger than the molecular distance but still fall short of the continuum limit^51,59,60—have not been systematically explored, experimentally or computationally.

In the long-time and long-wavelength limit, small density and pressure fluctuations in liquids are governed by linear hydrodynamics.^53–58,61 Extending the hydrodynamic description of collective liquid dynamics to mesoscopic scales—at finite wavelengths and frequencies—has been a preoccupation of many researchers in the past fifty years or so. In this pursuit, the contributions of nonhydrodynamic modes to the collective density fluctuations have long been recognized.^{9,52,56–59,61–71} Most notably, it has been shown by various means that the influence of structural relaxation on collective dynamics should be considered, particularly at low temperatures.^{9,52,56,57,62,68–72} In fact, precisely because structural relaxation dominates the long-time collective density–density correlations in the low-temperature regime (on length scales comparable to the wavelength of visible light), polarized dynamic light scattering has been widely used to extract the structural relaxation time of GF liquids.^44,72–81 On the other hand, extensive experimental,^82–95 theory, and simulation studies^{7,9,21,48,49,56,58,59,63–71,96–112} have been carried out in the broad context of generalized hydrodynamics.

Despite the herculean effort in the past several decades, there are still some holes in our understanding of the mesoscopic two-point collective dynamics of GF liquids. First, traditional experimental and computational “generalized hydrodynamics” investigations have directed most attention to relatively simple fluids, such as noble gases,^7,82,83 liquid metals,^{71,84–86,89–95,101,102,108} and generic Lennard-Jones liquids,^{21,100,105,107,111} which are either incapable of forming glasses or not in the supercooled state. Second, while polarized dynamic light scattering has been widely employed to study GF liquids, the accessible scattering angles are generally limited due to practical constraints. On the other hand, mapping out the mesoscopic collective dynamics of GF liquids over a wide range of time and length scales by inelastic neutron or x-ray scattering techniques is often technically challenging as well. Consequently, the basic phenomenology of mesoscopic collective dynamics of GF liquids has not been firmly established. The present molecular dynamics study aims to fill this void.

In addition, two general types of models have been extensively used in molecular dynamics simulations of GF liquids: one is based on binary mixtures^{22,23,113–116} and the other is based on single-component polymer melts.^117,118 Although the classical linear hydrodynamic theory of binary mixtures^54–57 is considered standard textbook knowledge, the similarities and differences between the collective dynamics of one- and two-component liquids have not been thoroughly characterized or fully appreciated.^{110,119–121} This issue will be discussed in Sec. VII.

Finally, despite the remarkable development of the generalized hydrodynamic theory^{9,21,56,58,64,65,104,106,107} and its application to GF liquids,^48,49,68 the connection between qualitative changes in collective dynamics and the dynamic crossover phenomenon has not been clearly demonstrated. The present work helps to bridge this gap in our knowledge by providing detailed molecular dynamics simulation data on two types of model GF liquids.

As a side note, a few words should be said about the term “nonhydrodynamic mode(s),” which we shall use extensively throughout this paper. “Nonhydrodynamic mode/behavior” generally refers to Q-independent, nonpropagating collective dynamics in liquids as opposed to the ordinary, propagating hydrodynamics modes.^{44,56,70,71,122–124} It is also sometimes called the “kinetic mode.”^{67,70,122,124,125} Different physical processes can lead to the emergence of nonhydrodynamic modes,⁷¹ but in our current investigation, the most relevant mechanism is the Mountain structural relaxation.^{44,53,56,57,68,72,76,126} Therefore, the terms nonhydrodynamic mode, Q-independent mode, and Mountain structural mode are used almost synonymously in this paper.

We choose to analyze the hydrostatic pressure–pressure space—time correlation C_pp(Q, t) for two reasons. First, it has been shown that the stress–stress correlation functions are spatially anisotropic, in general,^47–49,129 but the hydrostatic pressure–pressure correlation is an exception.¹²⁹ In other words, analysis of the isotropic function C_pp(Q, t) is significantly simpler than the anisotropic stress–stress correlations, such as the shear stress correlation function.^48,49 Second, the classical linear hydrodynamic theory has well-known predictions for the behavior of the normalized hydrostatic pressure–pressure correlation function C_pp(Q, t)/C_pp(Q).⁵⁵ They serve as a baseline for understanding the collective dynamics of GF liquids.

As stated in the preceding discussions, a major objective of this work is to establish the basic phenomenology of mesoscopic two-point collective dynamics of GF liquids in the vicinity of the Arrhenius crossover temperature T_A. To this end, the normalized density–density and hydrostatic pressure–pressure correlations are systematically mapped out using molecular dynamics simulations, and the predictions from the classical linear hydrodynamic theory are presented as references to guide the data interpretation. Our current focus is to highlight the qualitative change in mesoscopic collective dynamics around T_A. Quantitative analyses with the generalized hydrodynamic theories therefore are not pursued here. Nevertheless, some qualitative discussions are offered in Sec. VII.

We performed extensive molecular dynamics simulations on two types of representative GF liquids: coarse-grained (CG) bead-spring polymers^131,132 and Kob–Andersen (KA) Lennard-Jones (LJ) binary mixtures.^{22,23,133,134} This section describes the details of the simulations, which were carried out with the GPU-accelerated LAMMPS package.^136–137 All the models and results are presented in the reduced Lennard-Jones (LJ) units (e.g., σ for length, ɛ for energy, $τ=σm/ε$ for time, and ɛ/σ³ for stress), and the Boltzmann constant k_B is set to unity.

A coarse-grained (CG) bead-spring model^131,132 was employed in the polymer melt simulations. In this model, all beads interact with the truncated and shifted LJ potential^138, $ULJ(r)=4ε[(σr)12−(σr)6]−4ε[(σrc)12−(σrc)6]$⁠, r < r_c, where r_c = 2.5. The neighboring beads along the polymer chain are connected through the FENE potential^139, $UFENE(r)=−12kR02⁡ln[1−(r/R0)2]$⁠, with R₀ = 1.5 and k = 30. The polymer stiffness is controlled by a bending potential of the form¹³¹  U_bend(θ) = −a_θ sin²(b_θθ), 0 < θ<θ_c, where b_θ = 1.5, θ_c = π/b_θ, and the angle θ = cos⁻¹(b_j · b_j+1/|b_j‖b_j+1|), with b_j being the bond vector. Two different values of a_θ = 0 and a_θ = 2 were used to simulate fully flexible and semiflexible polymers, respectively.

In each case, the simulation contained M = 8000 linear chains with degree of polymerization z = 25. The initial configurations were generated by the hybrid MD/Monte Carlo method.¹⁴⁰ Specifically, thin polymer chains were first randomly placed in the simulation box, and the “size” of each bead in a polymer chain was then increased by using adaptive time steps and time-dependent potentials. Once the initial configurations were prepared, the system was equilibrated in the NPT ensemble at T = 2.0 and p = 0 for 5 × 10³ τ and subsequently quenched under a constant pressure of p = 0 to desired temperatures with a cooling rate of 10⁻⁶, which was slow enough to ensure full equilibration.¹³¹ The simulation temperatures were T = 1.5, 1.2, 1.0, 0.9, 0.8, 0.75, 0.7, 0.65, 0.6, 0.55, 0.5, and 0.48 for the fully flexible chains (a_θ = 0) and T = 1.5, 1.2, 1.0, 0.9, 0.8, 0.75, 0.7, 0.65, 0.6, 0.55, and 0.5 for the semiflexible chains (a_θ = 2). At each target temperature, the polymers were further relaxed in the NVT ensemble for more than 500× longer than the structural relaxation time before the production run.¹³² The reported structural and dynamical properties of the polymers S(Q), C_pp(Q), S(Q, t)/S(Q), and C_pp(Q, t)/C_pp(Q) were all analyzed in the NVT ensemble. The time step used in these molecular dynamics simulations was Δt = 0.002 τ. We have benchmarked our simulations with those of Xu et al.¹³² to make sure their results could be properly reproduced.

To investigate the behavior of atomic glass-formers, we adopted the Kob–Andersen (KA) model of LJ binary mixtures. Each simulated system consisted of a total of N = 60 000 particles with 80% large A particles and 20% B small particles.²² All the particles had the same mass and interacted through an LJ potential. Our choice of the interaction parameters followed the standard choice: σ_AA = 1.0, σ_AB = 0.8, σ_BB = 0.88, ɛ_AA = 1.0, ɛ_AB = 1.5, and ɛ_BB = 0.5. The interparticle interactions were truncated at 2.5 σ_AA.^133,134 The KA mixture was examined at two different pressures of p = 0 and p = 30. The simulation temperatures were T = 0.8, 0.75, 0.7, 0.65, 0.6, 0.55, 0.52, 0.5, 0.47, 0.45, 0.42, and 0.4 for p = 0 and T = 3.0, 2.5, 2.0, 1.7, 1.5, 1.4, 1.3, 1.2, 1.1, 1.05, and 1.0 for p = 30. It is worth noting that since the simulations were carried out under isobaric rather than isochoric conditions, our results should not be directly compared to those of Kob and Andersen.^22,23 Nevertheless, we also performed benchmark simulations to make sure that we were able to reproduce the results of Kob and Andersen under the same thermodynamic conditions. Before the production run, each system was first equilibrated at the target pressure and temperature for 5 × 10⁵ τ in the NPT ensemble and further relaxed in the NVT ensemble for more than 500× longer than the structural relaxation time. Similarly, the reported structural and dynamical properties of the polymers S(Q), C_pp(Q), S(Q, t)/S(Q), and C_pp(Q, t)/C_pp(Q) were analyzed in the NVT ensemble. The time step in the molecular dynamics simulations of the KA mixtures was Δt = 0.0025 τ.

This section details the calculations using the classical linear hydrodynamic theory.^53–58 The theory is used as a qualitative reference to identify the “nonhydrodynamic” density–density and pressure–pressure correlations. All the thermodynamic and transport coefficients required for computing the density–density and hydrostatic pressure–pressure correlations can be obtained, at least in principle, from molecular dynamics simulations. Various phenomenological extensions of the classical linear hydrodynamic theory are well documented in the literature.^9,56 However, we defer in-depth discussions of these approaches to a future work.

The hydrodynamic theory of single-component systems applies to our CG polymer simulations.^53,56 In the long-time, long-wavelength limit, the theory predicts that the normalized coherent intermediate scattering function S(Q, t)/S(Q) and hydrostatic pressure–pressure correlation function C_pp(Q, t)/C_pp(Q) are described by Eqs. (3) and (4), respectively. Since Γ = [b + (γ − 1)χ]/2, only four independent parameters, γ, χ, c, and b, are required for evaluating Eqs. (3) and (4).

To compute the specific heat ratio γ = c_p/c_V, simulations in NVT and NPT ensembles are required. Using the well-known fluctuation–dissipation relations, the specific heat at constant volume c_V can be estimated from the NVT simulation as $⟨δH2⟩NVT=NkBT2cV$⁠, whereas the specific heat at constant pressure c_p can be calculated from the NPT simulation as $⟨δ(H+PV)2⟩NPT=NkBT2cp$⁠.¹²⁸ Here, H is the total energy, N is the total number of beads, and V is the volume of the system.

The adiabatic sound velocity c can be computed from the specific heat ratio γ and the isothermal compressibility χ_T as $c=γ/(ρmχT)$⁠. In addition, χ_T is related to the zero-angle scattering as χ_T = lim_Q→0S(Q)/(ρk_BT).

In addition to γ, χ, c, and b, the mutual diffusion coefficient D_m is also needed for evaluating Eq. (9). Computing D_m from molecular dynamics simulations is not a straightforward task. For our current purpose, however, it suffices to use a simple estimation with reasonable accuracy,^145,146 D_m = x_BD_AA + x_AD_BB, where D_AA and D_BB are, respectively, the self-diffusion coefficients of particles A and B and x_A and x_B are the mole fractions.

This section presents the results of molecular dynamics simulations for the semiflexible (a_θ = 2) and fully flexible (a_θ = 0) coarse-grained polymers as well as the KA mixtures at p = 0 and p = 30. In Secs. IV A–IV C, the static structure factor S(Q) and hydrostatic pressure–pressure correlation function C_pp(Q) are first shown, followed by the spatiotemporal maps of S(Q, t)/S(Q) and C_pp(Q, t)/C_pp(Q) and further analyses of the correlation functions. Additional information regarding the self-dynamics, collective diffusion coefficient, and influence of threshold value on relaxation time is provided in the  Appendixes A and  B.

Figures 1(a) and 1(b) display the static density and pressure correlations of the semiflexible polymer (a_θ = 2) at various temperatures. The “first” peak of the static structure factor S(Q) is located at Q_max ≈ 7, where the structural relaxation time is calculated in our analysis. The flat region of S(Q) at low Qs defines the “mesoscopic” scale. The hydrostatic pressure–pressure correlation C_pp(Q) is qualitatively similar to S(Q) with an oscillatory region at high Qs and a flat region at low Qs. With the decrease in temperature, the first peaks of both S(Q) and C_pp(Q) become more pronounced, whereas the plateaus of S(Q) and C_pp(Q) at low Qs become lower. These observations are line with the previous reports of supercooled CG polymer melts in the literature.^147, Figures 1(c) and 1(d) present the collective dynamics of the semiflexible polymer at the first structural peak (Q = 7). Both the normalized coherent intermediate scattering function S(Q, t)/S(Q) and hydrostatic pressure–pressure correlation function C_pp(Q, t)/C_pp(Q) show a fast decay at high temperatures and exhibit a two-step relaxation as the temperature is lowered. Such a behavior of S(Q, t)/S(Q) resembles that of the self-intermediate scattering function F_s(Q, t) and is well documented.²⁴ In the high-Q region, S(Q, t)/S(Q) and F_s(Q, t) can be related by the Vineyard–de Gennes–Sköld convolution approximation,^130,148,149 and their resemblance is therefore not surprising. The characteristic “structural” relaxation time at Q = 7 is estimated as the correlation time when S(Q, t)/S(Q) or C_pp(Q, t)/C_pp(Q) decays to a constant value, e.g., 0.1. The results for the S(Q), C_pp(Q), and collective dynamics at Q = 7 of the flexible chain melts (a_θ = 0) are qualitatively similar and are shown in Fig. 2.

As explained in the Introduction and Sec. II, a key objective of this work is to establish the basic phenomenology of collective dynamics of GF liquids above and below the Arrhenius crossover temperature. For this purpose, the normalized density–density and hydrostatic pressure–pressure correlation functions are systematically mapped out over a wide range of correlation times and wave numbers. Figure 3 presents S(Q, t)/S(Q) and C_pp(Q, t)/C_pp(Q) from the molecular dynamics simulations of semiflexible polymers along with the predictions of the hydrodynamic theory [Eqs. (3) and (4)]. At high temperature (T = 1.0), the major features of density–density and hydrostatic pressure–pressure correlations on mesoscopic scales can be qualitatively understood with the hydrodynamic theory. The characteristic slope of −1 of the sound mode can be clearly seen in both S(Q, t)/S(Q) and C_pp(Q, t)/C_pp(Q). By contrast, below the Arrhenius crossover temperature T_A ≈ 0.75, excess correlations develop on mesoscopic scales, which is not envisioned by the classical linear hydrodynamic theory.^53,56 These “nonhydrodynamic” modes completely dominate the long-time relaxation behavior at small wave numbers. The observation of Q-independent collective dynamics in the CG polymer is consistent with the previous study of Aichele and Baschnagel.¹⁵⁰ It is worth noting that the collective intermediate scattering function S(Q, t)/S(Q) should not be confused with the single-chain dynamic structure factor, which is highly dependent on the wave number for Rouse chains.^24,152–153

The spatiotemporal maps of S(Q, t)/S(Q) and C_pp(Q, t)/C_pp(Q) for the fully flexible polymer (a_θ = 0) are qualitatively similar (Fig. 4). A de Gennes narrowing region with slow collective dynamics is found for S(Q, t)/S(Q) around the first structural peak Q ≈ 7 in the simulations. On the other hand, slow dynamics is also observed on mesoscopic scales [Fig. 4(a)]. At high temperatures (T > T_A), the mesoscopic collective dynamics can be qualitatively understood with the classical hydrodynamic theory [Fig. 4(b)]. However, as the temperature is lowered (T < T_A), the mesoscopic collective dynamics from the simulations become significantly different from the predictions by the hydrodynamic theory: nonhydrodynamic modes dominate the long-time density fluctuations, which are coupled to the local structural relaxation. The mesoscopic collective dynamics of hydrostatic pressure fluctuations displays a similar trend, as shown by the spatiotemporal maps of C_pp(Q, t)/C_pp(Q) in Figs. 4(c) and 4(d). It is helpful to point out that long-ranged shear stress correlation also emerges in glass-forming liquids at low temperatures,^48,49 which clearly resembles the behavior of the hydrostatic pressure–pressure correlation.

The slow, nonhydrodynamic density–density correlations are to the first approximation nothing but the so-called Mountain “structural” modes^{44,53,56,57,68,72,76,126} previously observed in polarized quasielastic light scattering experiments. However, light scattering experiments are limited to small momentum transfers and thus incapable of accessing the critical spatiotemporal regime at relatively high wave numbers. Systematic exploration of the nonhydrodynamic modes with coherent inelastic neutron scattering is not an easy task either as a result of many technical challenges. There are currently no comprehensive experimental data on mesoscopic collective dynamics slightly beyond the first structural peak in the vicinity of the Arrhenius crossover temperature. In this regard, computer simulations afford a unique opportunity to reveal how the nonhydrodynamic modes emerge in the dynamic crossover regime, where the experimental tools fall short.

To quantitatively study the nonhydrodynamic modes and confirm their relation to the local structural relaxation, we extract the Q-dependent, characteristic relaxation time τ(Q) from density–density and pressure–pressure correlations (Figs. 5 and 6). Similarly, instead of fitting the correlation functions using a phenomenological equation, the relaxation time τ(Q) is defined as S(Q, t = τ)/S(Q) = 0.1 or C_pp(Q, t = τ)/C_pp(Q) = 0.1. Examples of determining τ of the nonhydrodynamic long tail are given in Figs. 5(a), 5(b), 6(a), and 6(b). Consistent with the observations in Figs. 3 and 4, these plots show that the nonhydrodynamic modes become pronounced only at low temperatures (T < T_A). More remarkably, the characteristic relaxation times τ(Q) of the nonhydrodynamic modes exhibit the same temperature dependence as the structural relaxation time below T_A [Figs. 5(c) and 6(c)], suggesting a close coupling of local and mesoscopic dynamics in the dynamic crossover region. It should be noted that the density and pressure relaxation times differ by a prefactor, which partially stems from the “arbitrariness” in our method of defining τ. Figures 5(d) and 6(d) indicate that the nonhydrodynamic modes are indeed approximately independent of the wave number (τ ∼ Q⁰), lending additional support to their “structural” origin.

We note that the heat mode [the first term in Eq. (3)] does not contribute to the hydrostatic pressure–pressure correlation [Eq. (4)], and this should be at least partially explain some of the difference in the initial decays of S(Q, t)/S(Q) and C_pp(Q, t)/C_pp(Q) (Figs. 3–6). Nevertheless, the initial decays of S(Q, t)/S(Q) and C_pp(Q, t)/C_pp(Q) do exhibit completely opposite temperature dependence, which is not anticipated by the linear hydrodynamic theory. Generally speaking, the short-time relaxation of the “hydrodynamic” density–density or pressure–pressure correlation becomes faster with the decrease in temperature due to the increases in thermal diffusivity, kinematic longitudinal viscosity, and sound velocity. This is indeed observed in S(Q, t)/S(Q) but not in C_pp(Q, t)/C_pp(Q). At this point, the origin of the abnormal temperature dependence of the initial decay of the hydrostatic pressure–pressure correlation is unknown. The contribution of the hydrodynamic modes to the short-time density–density correlation leads to the nonmonotonic behavior of mesoscopic relaxation time in Fig. 5(c) at high temperatures. At low temperatures, however, the nonhydrodynamic mode dominates the mesoscopic dynamics at long times.

To demonstrate the universality of our observations, we turn to the case of KA mixtures. The static correlation functions S(Q) and C_pp(Q) of the KA mixture at p = 0 are shown in Figs. 7(a) and 7(b). Similar to the polymeric system, the first structural peak of the KA mixture is also located at Q_max ≈ 7. The relaxation behaviors of S(Q, t)/S(Q) and C_pp(Q, t)/C_pp(Q) at Q = 7 are shown in Figs. 7(c) and 7(d). Figures 8(a) and 8(b) present the static correlation functions S(Q) and C_pp(Q) at p = 30. Compared to the results at p = 0, the oscillatory features of C_pp(Q) at Qs become better defined at high pressure. Interestingly, S(Q, t)/S(Q) and C_pp(Q, t)/C_pp(Q) at the first structural peak (Q = 7) closely resemble each other in this case unlike the behavior observed in low pressure systems (Figs. 1, 2, and 7). Since the atomic virial stress depends on the separation of particle pairs [Eq. (5)], such a connection between S(Q, t)/S(Q) and C_pp(Q, t)/C_pp(Q) can be intuitively appreciated.

Representative spatiotemporal maps of S(Q, t)/S(Q) and C_pp(Q, t)/C_pp(Q) at p = 0 are shown in Fig. 9, whereas the results for the p = 30 system are provided in Fig. 10. Overall, the basic phenomenology of two-point collective density–density and hydrostatic pressure–pressure correlations in KA mixtures is incredibly similar to that of the CG polymers but with a slight twist: the hydrodynamic density–density correlations in binary mixtures have an additional relaxation mechanism due to mutual diffusion.^54,55,109 With the approximations described in Sec. V, the S(Q, t)/S(Q) predicted by the hydrodynamic theory can be simplified to Eq. (9). Compared to Eq. (3) for the single-component system, the new term for a binary mixture is exp(−D_mQ²t).¹⁵⁴ Because the mutual diffusion is also coupled to the structural relaxation, analysis of nonhydrodynamic modes from S(Q, t)/S(Q) becomes unattainable (Figs. 9 and 10). Fortunately, the hydrostatic pressure–pressure correlations are free from the influence of such a mutual diffusion term: the hydrodynamic theory⁵⁵ predicts that C_pp(Q, t)/C_pp(Q) of a binary mixture is of the same form as Eq. (4) for a single-component system. This allows us to probe the nonhydrodynamic modes from the hydrostatic pressure–pressure correlations.

The results of our analysis are shown in Figs. 9–12. The behavior of C_pp(Q, t)/C_pp(Q) in the KA mixture at p = 0 and p = 30 is highly similar to that of the CG polymers: at T > T_A, the pressure correlations on mesoscopic scales can be qualitatively understood with the classical linear hydrodynamic theory; on the other hand, nonhydrodynamic modes appear at T < T_A. These modes are coupled to the local structural relaxation and display the same temperature dependence, as indicated by the analysis of the characteristic relaxation time in Figs. 11(a) and 12(a). The appearance of nonhydrodynamic tail also coincides with the emergence of the slow structural relaxation process [Figs. 11(b) and 12(b)].

Our extensive molecular dynamics simulations with the CG polymer and KA models suggest that the qualitative change in two-point mesoscopic collective dynamics at the Arrhenius temperature is a universal feature of GF liquids: at high temperatures, the collective density and pressure fluctuations on mesoscopic scales are relatively “simple” and can be described approximately by ordinary hydrodynamics^56–58 although deviations from the classical theory is well known and has been extensively discussed and documented.^{9,51,52,56,58,62,65,67,71,105,124} Below the Arrhenius temperature, the “nonhydrodynamic” Q-independent mode dominates the slow relaxations on intermediate time and length scales.

From a theoretical perspective, the dominance of nonhydrodynamic, structural mode at low temperatures should be expected and in some sense is nothing surprising. The Q-independent (Q⁰), structural mode has been derived by various methods for the case of internal relaxation or viscous liquids.^52,62,68,69 Since such a mode is controlled by the structural relaxation, which becomes extremely slow in the deeply supercooled state, the nonhydrodynamic, Q⁰ mode is expected to prevail at low temperatures. However, there are no quantitative predictions as to when and in what spatiotemporal region the nonhydrodynamic structural mode takes over the ordinary hydrodynamic modes. In other words, the dramatic change in mesoscopic collective dynamics around the Arrhenius crossover temperature was not known previously.

Evidently, the aforementioned phenomenological treatments of the nonhydrodynamic mode cannot by itself explain the connection between the qualitative change in mesoscopic dynamics and the dynamic crossover phenomenon as no “glass physics” is embedded in these viscoelastic models. Additionally, the calculations are performed in the small-Q, hydrodynamic limit. The present simulations show, however, that the viscoelastic Q⁰ mode arises immediately below the first structural peak, where the static structure factor S(Q) becomes flat. In other words, the “small-Q limit” envisioned by the theories, in fact, starts from intermediate Qs on length scales only slightly larger than the cage size.

It is worth pointing out that our results are consistent with the observation of the increase in the Landau–Placzek ratio^72,155 below the Arrhenius temperature T_A by quasielastic light scattering as the Mountain structural mode contributes more prominently to the central line rather than the Brillouin doublet.^56,72 Most interestingly, our findings seem to be in line with the picture of coupled local and nonlocal relaxations at T < T_A,^157–160 where local topological excitations are expected to interact and show collective dynamics. In this context, we note that the wave number dependent viscosity η(Q) of GF liquids has been empirically modeled in terms of growing correlation lengths.¹⁶¹ The potential link between these ideas and our result remains to be explored.

The present work also demonstrates some key differences and similarities between the mesoscopic collective dynamics of one- and two-component liquids, which can be rationalized by the linear hydrodynamic theory of binary mixtures.^54–57 Despite the fact the theory has been around for more than half a century, to the best of our knowledge, there is no clear demonstration of its rudimentary predictions with computer simulations. Theoretically, the mutual diffusion in binary mixtures appears as an additional relaxation process for the density–density correlation [Eq. (9)] but only modifies the linewidth of the Brillouin doublet in the hydrostatic pressure–pressure correlation. Such expected behavior is nicely illustrated by the simulations of CG polymer melts and KA mixtures (Figs. 3, 4, 9, and 10). In particular, our simulations show that in the presence of the mutual diffusion, the nonhydrodynamic, Q⁰ mode becomes obscured in the density–density correlations of KA mixtures. Figures 9 and 10 also suggest that there is a strong coupling of density fluctuations to thermal processes in binary mixtures. These results should be useful for understanding the Q-independent slow dynamics in atomic and molecular systems, in general.¹¹⁰ 

Coarse-grained polymer melts and Lennard-Jones binary mixtures are two types of commonly used models for molecular dynamics simulations of GF liquids with the former being a single-component system and the later being a two-component system. While the collective hydrostatic pressure–pressure correlations in these systems exhibit great similarity, the density–density correlations show some clear, qualitative differences (Figs. 3, 4, 9, and 10) as discussed above. Additionally, the slow diffusion processes in polarized dynamic light scattering experiments of binary GF liquids^119–121 are often discussed without any reference to the mutual diffusion predicted by the ordinary linear hydrodynamic theory.^54–57 Such an oversight may lead to incorrect interpretation of the data. For example, it is possible that the slow diffusive process observed in some ionic liquids¹¹⁹ is simply due to mutual diffusion rather than the presence of Fischer clusters.⁴⁶ We hope that the data presented in this work can help raise awareness about the basic phenomenology of collective dynamics in one- and two-component systems.

Recently, Handle, Rovigatti, and Sciortino reported a computational study of the Q-independent slow dynamics in atomic and molecular systems,¹¹⁰ including the TIP4P/2005 water model,¹⁶² the Beest–Kramer–van Santen (BKS) silica model,¹⁶³ and the KA binary mixture model.^22,23 In this subsection, we briefly comment on the connection of the present investigation to the work of Handle et al. so that our results can be properly understood in a broader context.

In the small wave number limit, $ΩQ2=cT2Q2$ with $cT=kBT/(mS(0))$ being the isothermal sound velocity so that Eq. (19) becomes identical to Eq. (17).⁶⁸ In particular, when the term γ₀Q² describes the diffusion of momenta, γ₀ should take on the meaning of kinematic longitudinal viscosity, which has the same dimension as diffusivity. In fact, the idealized mode-coupling theory derivation by Latz and Schmitz starts from a Dyson-type equation that is incredibly similar to Eq. (19) except that their short-time relaxation is described by a Q-dependent “bare” longitudinal viscosity rather than a simple constant γ₀. Furthermore, in the case of exponential relaxation, $mQs(t)=As⁡exp(−t/τα)$⁠, a Q-independent, Mountain mode can be found by the perturbation theory as explained in the preceding discussion. In other words, the mode-coupling equation employed by Handle et al. in essence is similar to that of Latz and Schmitz.⁶⁸ 

Second, the results of the present investigation are consistent with the findings of Handle et al. for TIP4P/2005 water, BKS silica, and KA mixture. It is worth noting that while Handle et al. draw attention to the “gel mode”^164,165 in a variety of soft matter systems,^167–173 the emergence of nonhydrodynamic modes in the current context generally does not require the formation of a transient network.⁷¹ In particular, as discussed in Sec. VII A, the Mountain mode is a result of viscoelastic relaxation, which is completely general. In the case of KA mixtures, the presence of mutual diffusion masks the nonhydrodynamic mode in density–density correlation, which is in agreement with the theoretical expectation^55,109 and confirmed by both the present study and the work of Handle et al. Because of the network nature of silica, the interdiffusion process is suppressed,¹¹⁰ and the observation of the Q-independent mode is in line with the generalized hydrodynamic theory. Similarly, interdiffusion is unimportant for water and the presence of the Q-independent mode¹¹⁰ is consistent with the aforementioned theories for single-component systems.

In summary, the two-point collective dynamics of two model glass-forming liquids—coarse-grained bead-spring polymers and Lennard-Jones binary mixtures—have been systematically mapped out in a wide range of correlation times and wave numbers using molecular dynamics simulations. It is shown that the emergence of the nonhydrodynamic, Q-independent mode coincides with the dynamic crossover phenomenon. Despite the extensive theoretical studies on nonhydrodynamic mode in liquid dynamics, the observed connection between the qualitative change in mesoscopic collective dynamics and viscous slowing down remains to be explained. By providing comprehensive data on two model systems, this work also draws attention to the similarity and difference between the mesoscopic collective dynamics of one- and two-component liquids.

The research is supported by the U.S. Department of Energy (DOE), Office of Science, Office of Basic Energy Sciences, Early Career Research Program Award Grant No. KC0402010, under Contract Grant No. DE-AC05-00OR22725. The computational work was carried out at the Oak Ridge National Laboratory’s Center for Nanophase Materials Sciences, which is a DOE Office of Science User Facility. Our investigation used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725. We gratefully acknowledge helpful discussions with Professor K. S. Schweizer.

The authors have no conflicts to disclose.

Zhiqiang Shen: Data curation (lead); Formal analysis (lead); Investigation (equal); Methodology (equal); Writing – original draft (lead). Jan-Michael Y. Carrillo: Data curation (equal); Writing – review & editing (equal). Bobby G. Sumpter: Supervision (equal); Writing – review & editing (equal). Yangyang Wang: Conceptualization (lead); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Supervision (lead); Writing – review & editing (lead).

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

At low temperatures, the self-intermediate scattering function of GF liquids develops a “plateau” at the first structural peak—a phenomenon that is well known in the literature.^23,175 For CG polymers, the characteristic slope of −1/4 for Rouse dynamics can be observed at long correlation times and low Qs (Fig. 13).^51,153 Comparing Figs. 3 and 4 with Fig. 13, it is easy to see that the total collective normalized intermediate scattering function S(Q, t)/S(Q) ≈ 0 on time and length scales of the Rouse dynamics in contrast to the behavior of the single-chain dynamic structure factor. In other words, on the time and length scales of “polymeric motions” (e.g., Rouse motions, reptation, etc.), the density–density correlation due to local fluctuations has completely decayed and the liquid therefore can be considered to be “incompressible.”¹⁵³ 

In the main text, we compute the relaxation time τ from the density–density and pressure–pressure correlations according to the criteria: S(Q, t = τ)/S(Q) = 0.1 and C_pp(Q, t = τ)/C_pp(Q) = 0.1 or 0.07. To demonstrate that the qualitative conclusions of this work do not critically depend on the “threshold value” in such calculations, we present in Fig. 16 the relaxation times computed according to S(Q, t = τ)/S(Q) = 0.2 and C_pp(Q, t = τ)/C_pp(Q) = 0.05 for the semiflexible polymer (a_θ = 2) and KA mixture (p = 0). Indeed, the results in Fig. 16 are similar to those in Figs. 5(c) and 11(a), confirming the close relation between local structural relaxation and mesoscopic relaxation at low temperatures (T < T_A).