The collective density–density and hydrostatic pressure–pressure correlations of glass-forming liquids are spatiotemporally mapped out using molecular dynamics simulations. It is shown that the sharp rise of structural relaxation time below the Arrhenius temperature coincides with the emergence of slow, nonhydrodynamic collective dynamics on mesoscopic scales. The observed long-range, nonhydrodynamic mode is independent of wave numbers and closely coupled to the local structural dynamics. Below the Arrhenius temperature, it dominates the slow collective dynamics on length scales immediately beyond the first structural peak in contrast to the well-known behavior at high temperatures. These results highlight a key connection between the qualitative change in mesoscopic two-point collective dynamics and the dynamic crossover phenomenon.

## I. INTRODUCTION

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.^{50} To overcome the shortcoming associated with the traditional practice of analyzing intermediate scattering functions at discrete wave numbers, we adopt a “two-dimensional” approach^{51} 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.

## II. SURVEY OF EXISTING STUDIES ON MESOSCOPIC TWO-POINT COLLECTIVE DYNAMICS OF GF LIQUIDS

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,^{71} 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.

## III. APPROACH OF THE CURRENT WORK

*S*(

**Q**,

*t*)] and hydrostatic pressure–pressure [

*C*

_{pp}(

**Q**,

*t*)] correlation functions in the reciprocal space,

*N*is the total number of particles,

*p*

_{i}(

*t*) is the hydrostatic pressure element

^{127}defined by the diagonal components of the atomic pressure tensor for particle

*i*at time

*t*, and

**r**

_{i}(

*t*) is the particle position.

*S*(

**Q**) and

*C*

_{pp}(

**Q**) denote, respectively, the static density–density and pressure–pressure correlations:

*S*(

**Q**) ≡

*S*(

**Q**, 0) and

*C*

_{pp}(

**Q**) ≡

*C*

_{pp}(

**Q**, 0). Using

**Q**s that commensurate with the periodic boundary conditions in the simulations,

^{128}we systematically map out the normalized coherent intermediate scattering function

*S*(

*Q*,

*t*)/

*S*(

*Q*) and press-pressure correlation function

*C*

_{pp}(

*Q*,

*t*)/

*C*

_{pp}(

*Q*) on a dense grid of correlation times and wave numbers. Such a choice of wave numbers ensures that the minimum-image convention is implicitly implemented in the calculations of collective structures and dynamics in the reciprocal space.

^{128}Additionally, all of the simulation systems are sufficiently large so that no obvious finite-size effects have been observed. It is worth noting that the reciprocal space is the natural choice for analyzing the mesoscopic collective dynamics in liquids. While real-space analysis has its own appeal, it is easy to see that the hydrodynamic sound modes produce rather complicated functional forms in real space.

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.^{129} 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*).^{55} They serve as a baseline for understanding the collective dynamics of GF liquids.

*S*(

*Q*) and pressure correlation

*C*

_{pp}(

*Q*) exhibit oscillations (e.g.,

*Q*⪆ 3 for the coarse-grained polymer) and the mesoscopic regime where

*S*(

*Q*) and

*C*

_{pp}(

*Q*) are essentially flat (

*Q*⪅ 3). The most prominent feature of local relaxation at high wave numbers is the de Gennes narrowing phenomenon,

^{130}where the apparent relaxation time of

*S*(

*Q*,

*t*)/

*S*(

*Q*) displays a maximum around the first structural peak of

*S*(

*Q*). On the other hand, according to the conventional wisdom, in the absence of internal or slow viscoelastic relaxations, the density and pressure fluctuations in the long-time and wavelength limit are governed by ordinary hydrodynamics.

^{56}For example, the classical linear hydrodynamic theory

^{9,52–58}predicts that

*S*(

*Q*,

*t*)/

*S*(

*Q*) and

*C*

_{pp}(

*Q*,

*t*)/

*C*

_{pp}(

*Q*) of a single-component fluid are described by the following equations:

*γ*is the specific heat ratio,

*χ*is the thermal diffusivity,

*c*is the adiabatic sound velocity,

*b*is the kinematic longitudinal viscosity, and Γ is the acoustic attenuation coefficient. The first term on the right-hand side of Eq. (3) is the “heat mode” from entropy fluctuations at constant pressure, which is not present in

*C*

_{pp}(

*Q*,

*t*)/

*C*

_{pp}(

*Q*). The sinusoidal and cosinusoidal terms constitute the “sound mode” associated with the Brillouin doublet. All the thermodynamic and transport coefficients in Eqs. (3) and (4) can be evaluated directly from molecular dynamics simulations.

^{51,107}The details of the linear hydrodynamic theory calculations will be described in Sec. V.

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.

## IV. SIMULATION METHODS

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, $\tau =\sigma m/\epsilon $ for time, and *ɛ*/*σ*^{3} for stress), and the Boltzmann constant *k*_{B} is set to unity.

### A. Coarse-grained polymers

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\epsilon [(\sigma r)12\u2212(\sigma r)6]\u22124\epsilon [(\sigma rc)12\u2212(\sigma 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)=\u221212kR02\u2061ln[1\u2212(r/R0)2]$, with *R*_{0} = 1.5 and *k* = 30. The polymer stiffness is controlled by a bending potential of the form^{131} *U*_{bend}(*θ*) = −*a*_{θ} sin^{2}(*b*_{θ}*θ*), 0 < *θ*<*θ*_{c}, where *b*_{θ} = 1.5, *θ*_{c} = *π*/*b*_{θ}, and the angle *θ* = cos^{−1}(**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.^{140} 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^{3} *τ* and subsequently quenched under a constant pressure of *p* = 0 to desired temperatures with a cooling rate of 10^{−6}, which was slow enough to ensure full equilibration.^{131} 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.^{132} 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.*^{132} to make sure their results could be properly reproduced.

### B. Kob–Andersen binary mixtures

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.^{22} 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^{5} *τ* 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 *τ*.

### C. Hydrostatic pressure–pressure correlation

*i*at position

**r**

_{i}and time

*t*can be evaluated as

^{141}

^{,}

**r**

_{ij}(

*t*) ≡

**r**

_{j}(

*t*) −

**r**

_{i}(

*t*),

**F**

_{ij}(

*t*) is the force applied on particle

*i*by particle

*j*, and $rij\alpha (t)$ and $Fij\beta (t)$ stand for the

*α*and

*β*(

*α*,

*β*=

*x*,

*y*,

*z*) components of

**r**

_{ij}(

*t*) and

**F**

_{ij}(

*t*), respectively, $Wi\alpha \beta (t)\u226112\u2211j\u2260irij\alpha (t)Fij\beta (t)$. The kinetic contribution to stress is small

^{127,142}and therefore not included in our calculations. Additionally, only interchain LJ interactions are considered for the CG polymers. The hydrostatic pressure “element”

^{127}

*p*

_{i}(

*t*) is defined by the diagonal components of $Wi\alpha \beta (t)$, $pi(t)=(Wixx(t)+Wiyy(t)+Wizz(t))/3$. The spatial Fourier transform of $\sigma i\alpha \beta (r,t)$ is $\sigma i\alpha \beta (Q,t)\u2261\u222b\sigma i\alpha \beta (r,t)eiQ\u22c5rdr$. The hydrostatic pressure–pressure correlation function

*C*

_{pp}(

**Q**,

*t*) is thus defined as

## V. HYDRODYNAMIC THEORY CALCULATIONS

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.

### A. Single-component systems

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 $\u27e8\delta H2\u27e9NVT=NkBT2cV$, whereas the specific heat at constant pressure *c*_{p} can be calculated from the NPT simulation as $\u27e8\delta (H+PV)2\u27e9NPT=NkBT2cp$.^{128} Here, *H* is the total energy, *N* is the total number of beads, and *V* is the volume of the system.

*χ*can be obtained from the thermal conductivity

*κ*as

*χ*=

*κ*/(

*ρc*

_{p}/

*m*), where

*ρ*is the bead number density and

*m*is the bead mass. The thermal conductivity can be evaluated from the energy current autocorrelation function $\u27e8j\alpha \epsilon (t)j\alpha \epsilon (0)\u27e9$ in the NVE simulation by using the Green–Kubo formula,

^{143,144}

The adiabatic sound velocity *c* can be computed from the specific heat ratio *γ* and the isothermal compressibility *χ*_{T} as $c=\gamma /(\rho m\chi T)$. In addition, *χ*_{T} is related to the zero-angle scattering as *χ*_{T} = lim_{Q→0}*S*(*Q*)/(*ρk*_{B}*T*).

*b*+ (

*γ*− 1)

*χ*]/2 and kinetic longitudinal viscosity

*b*=

*η*

_{L}/(

*ρm*) requires information about the longitudinal viscosity

*η*

_{L}, which can be evaluated using the Green–Kubo relation,

^{128}

*δp*

^{αα}is the fluctuation of the diagonal component of the pressure tensor in the NVE simulation,

*δp*

^{αα}(

*t*) =

*p*

^{αα}(

*t*) − ⟨

*p*

^{αα}⟩, with

*α*=

*x*,

*y*, or

*z*.

### B. Binary mixtures

^{54–56}predicts a complicated functional form for the collective density–density correlation

*S*(

*Q*,

*t*)/

*S*(

*Q*) of a binary mixture, which we shall not try to reproduce here. While all the thermodynamic and transport coefficients involved in this expression, in principle, can be evaluated from molecular dynamics simulations, the computational cost associated with such calculations is prohibitively high. However, with some reasonable approximations, the equation can be significantly simplified. For dense liquids, we can assume the ratio of adiabatic compressibility

*β*

_{s,c}and isothermal compressibility

*β*

_{T,μ}is approximately unity,

*β*

_{s,c}/

*β*

_{T,μ}≈ 1. Additionally, we may assume that the thermal diffusion ratio

*k*

_{T}and barodiffusion coefficient

*k*

_{p}are small. Under these assumptions, it is easy to show that the following simplified expression can obtained:

*D*

_{m}is the mutual diffusion coefficient and Γ ≈ [

*b*+ (

*γ*− 1)

*χ*]/2. On the other hand, the pressure–pressure correlation function

*C*

_{pp}(

*Q*,

*t*)/

*C*

_{pp}(

*Q*) of a binary mixture, as predicted by the hydrodynamic theory,

^{55}takes exactly the same form as Eq. (4). Similarly, in the current qualitative analysis, we ignore the modification of the decay rate for the hydrostatic pressure–pressure correlation with the assumption of small

*k*

_{T}and

*k*

_{p}.

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*_{B}*D*_{AA} + *x*_{A}*D*_{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.

## VI. RESULTS

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.

### A. Coarse-grained polymers

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 *Q*s defines the “mesoscopic” scale. The hydrostatic pressure–pressure correlation *C*_{pp}(*Q*) is qualitatively similar to *S*(*Q*) with an oscillatory region at high *Q*s and a flat region at low *Q*s. 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 *Q*s 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.^{24} 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.^{150} 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*^{0}), 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.

### B. Kob–Andersen binary mixtures

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 *Q*s 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*_{m}*Q*^{2}*t*).^{154} 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^{55} 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)].

## VII. DISCUSSION

### A. Relation to current theoretical understanding

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*^{0}), 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*^{0} 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.

*mechanistic*explanation for the qualitative change in mesoscopic dynamics around the crossover temperature. For example, we may consider a generic expression for the normalized dynamic structure factor [normalized Laplace transform of the intermediate scattering function

*S*(

*Q*,

*z*)/

*S*(

*Q*)] for a single-component fluid,

^{106}

^{,}

*D*(

*Q*,

*z*) is a generalized diffusion function. Equation (10) can be derived by a phenomenological generalization of the Fickian law of diffusion or the projection operator technique.

^{106}In the hydrodynamic limit,

*D*(

*Q*,

*z*) takes the following form if the entropy fluctuations are neglected (this assumption is not essential but significantly simplifies the discussion),

^{9,106}

*b*=

*η*

_{L}/(

*ρm*) is the kinematic longitudinal viscosity defined earlier. It follows that

*bk*≪ 2

*c*

_{T}, two poles are located at

*Q*-independent mode is found for a viscous liquid in the intermediate wave number region, where

*bQ*

^{2}≫

*z*,

*η*

_{L}and modulus

*M*

_{L}is used,

*η*

_{L}=

*M*

_{L}

*τ*

_{α}. The preceding analysis, in essence, is the same as the one presented by Schweizer and Saltzman,

^{69}which is a small-

*Q*approximation of Eq. (12). Since Eq. (12) simply describes “isothermal” ordinary hydrodynamics in the small-

*Q*limit,

^{9}this derivation does not explicitly introduce any structural relaxation process and is thus not fully consistent from both the physical and mathematical points of view. A natural and better approach to illustrate the nonhydrodynamic mode is to consider the viscoelastic model

^{52,56,57,62,73,103}by taking into account the relaxation of kinematic longitudinal viscosity,

*b*

_{∞}and

*b*

_{0}are kinematic longitudinal viscosities in the high- and low-frequency limits, respectively. In this simple model, the dynamic longitudinal modulus is assumed to follow the Maxwell behavior, $M(z)=z\eta L(z)\u221dz/(z+\tau \alpha \u22121)$. Instead of Eq. (12), we now have

*Q*-independent mode with a pole at $zM=\u2212\tau \alpha \u22121$.

^{68}using the mode-coupling theory has a similar mathematical structure. In their case, the normalized dynamic structure factor is of the form

*m*(

*z*) is the frequency dependent relaxation function. They assume that

*m*(

*z*) has the asymptotic behavior of

*m*(

*z*) ∝

*z*

^{−1}at high frequencies and

*m*(

*z*) ∝

*τ*

_{α}at low frequencies, which is exactly what one would expect for a Maxwell model. Therefore, perturbation analysis also yields a Mountain mode at $zM=\u2212\tau \alpha \u22121$ along with two Brillouin lines.

^{68}In passing, it is useful to note that the nonhydrodynamic mode for hydrostatic pressure–pressure can be demonstrated in a similar fashion with a viscoelastic extension of ordinary hydrodynamics. In the case of shear stress correlation function, the Maxwell model has been invoked to describe the behavior of GF liquids.

^{48,49}In particular, the emergence of long-ranged shear stress correlations in glass-forming liquids at low temperatures

^{48,49}is closely connected to the results of the current investigation.

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*^{0} 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 *Q*s 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.^{161} The potential link between these ideas and our result remains to be explored.

### B. Single-component systems vs binary mixtures

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*^{0} 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.^{110}

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^{119} is simply due to mutual diffusion rather than the presence of Fischer clusters.^{46} 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.

### C. Relation to the recent work of Handle and co-workers

Recently, Handle, Rovigatti, and Sciortino reported a computational study of the *Q*-independent slow dynamics in atomic and molecular systems,^{110} including the TIP4P/2005 water model,^{162} the Beest–Kramer–van Santen (BKS) silica model,^{163} 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.

*et al.*numerically solved a simplified mode-coupling theory equation to fit their simulation data in the

*time*domain,

*ϕ*

_{Q}(

*t*) ≡

*S*(

*Q*,

*t*)/

*S*(

*Q*), $\Omega Q2=vT2Q2/S(Q)$ with $vT=kBT/m$ being the thermal velocity,

*γ*

_{0}is a constant associated with the damping of fast microscopic dynamics, and $MQs(t)\u2261AsQ2\u2061exp[\u2212(t/\tau \alpha )\beta ]=mQs(t)Q2$ describes the structural relaxation with

*A*

_{s}being a numerical constant,

*τ*

_{α}being the structural relaxation time, and

*β*being the stretch exponent. Using Laplace transform, it is straightforward to find the solution of the above integro-differential equation in the

*frequency*domain,

^{58}

*ϕ*

_{Q}(

*t*), $MQs(t)$, and $mQs(t)$ respectively.

In the small wave number limit, $\Omega Q2=cT2Q2$ with $cT=kBT/(mS(0))$ being the isothermal sound velocity so that Eq. (19) becomes identical to Eq. (17).^{68} In particular, when the term *γ*_{0}*Q*^{2} describes the diffusion of momenta, *γ*_{0} 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 *γ*_{0}. Furthermore, in the case of exponential relaxation, $mQs(t)=As\u2061exp(\u2212t/\tau \alpha )$, 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.^{68}

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.^{71} 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,^{110} 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^{110} is consistent with the aforementioned theories for single-component systems.

## VIII. SUMMARY

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.

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

### APPENDIX A: SELF-DYNAMICS AND COLLECTIVE DIFFUSION COEFFICIENT

*F*

_{s}(

**Q**,

*t*) is defined as

*N*is the total number of particles in the system and

**r**

_{i}(

*t*) is the position of the

*i*th particle at time

*t*. Following our previous approach,

^{51,174}the self-intermediate scattering functions of the CG polymer melts and KA mixtures are spatiotemporally mapped out on dense grids of correlation times and wave numbers. Representative results are shown in Figs. 13 and 14.

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 *Q*s (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.”^{153}

*Q*s is dominated by the mutual diffusion of large and small particles. To provide further information on the diffusive behavior, the collective,

*Q*-dependent relaxation time

*τ*is evaluated according to the criterion

*S*(

*Q*,

*t*=

*τ*)/

*S*(

*Q*) = 0.1 and presented in Fig. 15. The collective diffusion coefficient can be determined from the low-

*Q*relaxation time as

### APPENDIX B: INFLUENCE OF THRESHOLD VALUE ON RELAXATION TIME

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}).

## REFERENCES

*Metastable Liquids: Concepts and Principles*

*Jamming and Rheology: Constrained Dynamics on Microscopic and Macroscopic Scales*

*Complex Dynamics of Glass-Forming Liquids: A Mode-Coupling Theory*

*Dynamical Heterogeneities in Glasses, Colloids, and Granular Media*

*Statistical Physics of Liquids at Freezing and Beyond*

*Dynamics of the Liquid State*

*Molecular Hydrodynamics*

*Dynamic Light Scattering: With Applications to Chemistry, Biology, and Physics*

*Theory of Simple Liquids: With Applications to Soft Matter*

*Photon Correlation and Light Beating Spectroscopy*

*Neutron Inelastic Scattering*

*Physical Acoustics*

*Computer Simulation of Liquids*

*The Theory of Polymer Dynamics*

(*γ* − 1)/*γ* ≈ *γ* − 1, because *γ* ≈ 1 for dense liquids.