Fast Vibrational Modes and Slow Heterogeneous Dynamics in Polymers and Viscous Liquids

Many systems, including polymers and molecular liquids, when adequately cooled and/or compressed, solidify into a disordered solid, i.e., a glass. The transition is not abrupt, featuring progressive decrease of the microscopic mobility and huge slowing down of the relaxation. A distinctive aspect of glass-forming materials is the microscopic dynamical heterogeneity (DH), i.e., the presence of regions with almost immobile particles coexisting with others where highly mobile ones are located. Following the first compelling evidence of a strong correlation between vibrational dynamics and ultraslow relaxation, we posed the question if the vibrational dynamics encodes predictive information on DH. Here, we review our results, drawn from molecular-dynamics numerical simulation of polymeric and molecular glass-formers, with a special focus on both the breakdown of the Stokes–Einstein relation between diffusion and viscosity, and the size of the regions with correlated displacements.


Introduction
When polymers, liquids, biomaterials, metals and molten salts are cooled or compressed, if the crystallization is avoided, they freeze into a microscopically disordered solid-like state, a glass [1][2][3]. On approaching the glass transition from states with high fluidity, the viscosity exhibits a huge increase of more than 10 orders of magnitude [1,2], along with the parallel decrease of the diffusivity [3,4]. Correspondingly, at microscopic level, solid-like behaviour becomes apparent, e.g., it is observed that a particle spends increasing time within the cage formed by the first neighbours where it rattles with amplitude u 2 1/2 on picosecond time scales [5]. This temporary trapping is rather persistent and the particle has average escape time, the structural relaxation time τ α , which increases from a few picoseconds in the low-viscosity liquid up to thousands of seconds close to the glass transition [6]. The quantity u 2 appears in the expression of the Debye-Waller (DW) factor, which, assuming harmonicity and isotropy of the thermal motion, takes the form exp −q 2 u 2 /3 , where q is the absolute value of the scattering vector [7]. Researchers investigating the cage motion in viscous liquids usually refer, as a metonym, to u 2 as the DW factor too, e.g., see the work in [8][9][10]. To keep maintain similarity with this literature, the same convention is adopted here.
The transition from a liquid to a glass is accompanied by the growth of transient domains which exhibit different mobility, e.g., see Figure 1. The phenomenon is usually dubbed "dynamical heterogeneity" (DH) and has been extensively studied, e.g., see the reviews in [4,6,11,12] and topical papers [13][14][15][16]. The size of the domains is relatively small involving approximately 10 molecule diameters [11], corresponding to a few nanometres [14]. On a more general ground, the size of DH domains is strictly related to the possible presence of characteristic length scales in glass-forming systems. Starting with the seminal paper by Adam and Gibbs, who invoked the presence of "cooperatively rearranging regions" in viscous liquids [17], increasing interest has been devoted to identifying possible growing length scales as mobility decreases [18,19]. A broad classification in terms of either static or dynamic length scales is usually used. Static (thermodynamic) length scales are determined by the free-energy landscape, whereas dynamic length scales are set by the rules governing the time evolution of the system and extracted from finite-time behaviour of time-dependent correlation functions and associated susceptibilities [6]. Even if growing static length scales have been reported by experiments [20] and simulations [21], there is still debate if they control the glass transition [22]. It is still not clear to what extent dynamic correlations are the consequence, or the primary origin of, slow dynamics occurring close to the glass transition [19]. . Bright yellow particles have squared displacements no less than 1. Notice that the two states have comparable mean square displacement (∼0.21, homogeneous state; ∼0.28, heterogeneous state) but rather different relaxation times τ α (∼9, homogeneous state; ∼1550, heterogeneous state). Homogeneous, i.e., position-independent, dynamics of the monomers is an aspect of systems with fast relaxation. Conversely, in the presence of heterogeneous dynamics, clusters of particles with extremely high mobility coexist with nearly immobile ones, slowing down the relaxation.
Even if rooted at nanometric length scales, DH exerts clear influence at macroscopic level. One widely studied phenomenon is the breakdown of the Stokes-Einstein (SE) relation involving the diffusion coefficient D and the shear viscosity η (the more debated analogous phenomenon involving the rotational diffusion, where the breakdown is revealed [23,24] or missing [11], will not be considered here). For a single particle with radius R moving in a homogeneous fluid with viscosity η at temperature T, the SE relation states that k B denotes the Boltzmann constant and ζ denotes a number depending on the boundary condition between the fluid and the particle [25]. Under a no-slip condition, ζ = 6. Roughly, the SE law states that the quantity k B T/Dη is a constant of the order of the size of the diffusing particle. Remarkably, despite its macroscopic derivation, SE also well accounts for the self-diffusion of many monoatomic and molecular liquids, provided the viscosity is low ( 10 Pa ·s) [26]. On the other hand, the finite diffusion coefficient of guest atoms in solid hosts, where viscous transport is missing strongly, suggest the SE failure close to the solidification occurring at the glass transition. In fact, a common feature of several fragile glass formers is the SE breakdown for increasing viscosity. The failure manifests itself as a partial decoupling between the diffusion and the viscosity, in the sense that D −1 increases less than η [4,11,[27][28][29][30]. The decoupling is well accounted for by the fractional SE D ∼ η −κ [31], where the non-universal exponent κ falls in the range 0.5 ≤ κ < 1 [15]. The usual interpretation of the SE breakdown relies on DH and the subsequent presence of a spatial distribution of characteristic relaxation times τ close to the glass transition [6,11,31]. The neat argument is that, although the viscosity is more sensitive to the longest relaxation times, the diffusivity is influenced by the shortest ones. As the shape of the distribution tends to widen on approaching the glass transition, the gap between D −1 and η increases as well, leading to the SE breakdown [4]. Diffusion, viscous transport, and structural relaxation involve time scales that are much longer that than the typical vibrational time t of the particle rattling in the cage of the first neighbours, typically a few picoseconds. The diffusion coefficient is expressed as D = 6δ 2 /τ D , where τ D is the minimum time ensuring that the particle random displacements at a pace τ D are statistically independent with finite mean square value δ 2 [25]. On the other hand, the viscous flow requires the relaxation of the shear stress fluctuation, which occurs in a Maxwell time τ M = η/G, where G is the intermediate-time shear modulus [32]. On approaching the glass transition, t τ D , τ M , τ α . Despite the huge difference in time scales, earlier [33] and later theoretical studies [5,[34][35][36][37][38][39][40], and experimental ones [41], addressed the rattling process within the cage to understand the slow dynamics, rising a growing interest on the DW factor [8,29,30,[42][43][44][45][46][47][48][49][50][51][52][53][54][55][56][57][58][59][60][61]. Within this context, most interest has been devoted to the correlations between DW factor and the structural relaxation time τ α , which are found to be strong and encompassed by a universal master curve [47]: An analytical expression of the master curve is derived in Section 2. Alternative forms of the master curve are reported by Douglas and coworkers [8][9][10]. Correlations between DW factor and the structural relaxation time τ α are found in polymers in bulk [30,[47][48][49] and thin films [61], binary atomic mixtures [48,55], colloidal gels [52], antiplasticised polymers [8,9], water [57] and water-like models [59,60]. The DW factor also provided an alternative interpretation of the so-called thermodynamic (or temperature/density) scaling [58]. The correlation between structural relaxation and DW factor has been inspected in the experimental data concerning several glass-formers in a wide range of fragility-the steepness index m defined by Angell [1] (20 ≤ m ≤ 191), including polymers, van der Waals and hydrogen-bonded liquids, metallic glasses, molten salts and the strongest inorganic glass-formers [47,50,51,[55][56][57][58].
The structural relaxation time τ α is an average quantity which is certainly affected by DH but not in a straightforward way. Nonetheless, given the scaling expressed by Equation (2), it is legitimate to wonder if DH and fast vibrational dynamics exhibit correlations. Working in that direction, we have found positive answers and the present paper collects and reviews a selected part of our results, with a focus on the breakdown of the SE law [29,30]. Even if strictly related, we will not discuss here a study concerning ultrathin molecular films with strong mobility gradients analogous to DH, where the same scaling observed in bulk, Equation (2), has been revealed [61].
Our approach relies on the increasing evidence that the master curve, Equation (2), is a manifestation of a more fundamental correlation between the vibrational dynamics and the slow relaxation. It may be presented in the following terms. Let us consider a generic space-and time-dependent correlation function C(x 1 , t 1 ; x 2 , t 2 ), where x denotes a configuration of the liquid at a given time t, i.e., the set x of all the positions of the elementary microscopic particles (monomers, atoms, molecules, etc.). For steady states, C(x 1 , t 1 ; x 2 , t 2 ) depends on the time difference t = t 2 − t 1 . Let us set t 1 = 0 and define C(x 0 ; x, t) ≡ C(x 0 , 0; x, t). If two states, labelled by a and b, have equal DW factor, the correlation function, when evaluated over the two states, has coincident time evolution at least between the typical vibrational time scale t and τ α [49]. Said otherwise, for t * t τ α , it holds [32,47,[49][50][51][53][54][55]: In selected systems, Equation (3) holds beyond τ α and extends up to the diffusive regime, e.g., unentangled polymers and atomic binary mixtures [29,30,49,55].
Our studies were prompted by the finding by Widmer-Cooper and Harrowell that DH are predicted by particle displacements at short times [44]. However, it must be stressed that our DW factor is evaluated within the vibrational time scale t and not the time scale in [44], which is approximately one order of magnitude longer, a choice leading to differences for states with low viscosity.
The review outlines a model of the slow heterogeneous relaxation and transport in terms of vibrational dynamics in Section 2. The model is presented for completeness, but it is not essential to the understanding of the simulation results discussed in the rest of the paper. Later, a broad introduction to relaxation and transport in polymeric melts, and the correlation with the vibrational fast dynamics is given in Sections 3 and 4, respectively. The signatures identifying the presence of heterogeneous dynamics are discussed in Section 5. The SE breakdown is presented in Section 6, with a final discussion on the length scale of the mutual influence between particle displacements in Section 7.

A model of the Slow Heterogeneous Relaxation and Transport in Terms of Vibrational Dynamics
An in-depth, microscopic understanding of the link between the fast and slow dynamics is still missing, even if the impact of anharmonicity has been noted [43,58,62]. Here, we present a model, extending first seminal ideas [34], where the key role is played by the DW factor u 2 , which is a single-particle quantity. Alternative pictures, in terms of the same quantity, are known [8][9][10]. Notice that, even if a single-particle quantity, the DW factor encodes information on collective dynamics and spatially extended cooperative phenomena [32,53,54,[62][63][64].
At the present level of development, the model delivers expressions of the diffusion coefficient and the structural relaxation time in terms of the DW factor. It also accounts for the nonexponential character of the relaxation, an aspect which will be not presented here. However, even if it incorporates some consequence of DH, i.e., the presence of a wide distribution of relaxation times p(τ), it does not cover any spatial aspect related to DH, which instead has been revealed by the simulations, as we will see in Sections 5.1 and 7, and accounted for by Equation (3).

Relaxation Time
A first basis to connect fast and slow degrees of freedom was developed by Hall and Wolynes who, assuming that atomic motion is restricted to cells, pictured the glass transition as a freezing in an aperiodic crystal structure [34]. As a result, the viscous flow is described in terms of activated jumps over energy barriers ∆E ∝ k B Ta 2 / u 2 , where a is the displacement to reach the transition state. The usual rate theory leads to the Hall-Wolynes equation: Equation (4) has the form of Equation (2). A very similar relation was derived by Buchenau and Zorn, in terms of soft vibrational modes [41]. Equation (4) is expected to fail when τ α becomes comparable to the typical rattling times of each atom in the cage, corresponding to picosecond timescales. This condition is quite mild, e.g., in selenium it occurs~100 K above the melting temperature [41].
One basic assumption of Equation (4) is that the distance to reach the transition state has a characteristic value a. Actually, this length scale is dispersed. To constrain the related distribution, p(a 2 ), it is assumed that the latter does not depend on the state parameters such as the temperature, the density or the interacting potential. This complies with the spirit of the work in [34], where the a distance is said to be mostly controlled by the geometrical packings. It is also known that, irrespective of the relaxation time, τ α , the average distance moved by the relaxing unit within τ α is approximately the same, i.e., a fraction of the molecular diameter [1]. Averaging Equation (4) over the distribution p(a 2 ) yields the structural relaxation time Note that Equation (6) assumes that the distribution of the relaxation times is mainly due to the distribution of the displacement to reach the transition state in the different local environments, whereas the average DW factor u 2 is taken as homogeneous across the sample. This viewpoint relies on the picture that relaxation is related to long wavelength soft modes [41,46]. Support has been provided by the strong correlation observed in glass-formers between u 2 and the elastic modulus under quasi-static mechanical equilibrium [32].
As a suitable choice, the distribution of the squared distances p(a 2 ) is taken as a truncated Gaussian form [47,48] where A is the normalization ensuring ∞ 0 p(a 2 ) da 2 = 1 and a 2 min is the minimum displacement to reach the transition state. Given the weak influence, and to get rid of an adjustable parameter, one takes a 2 min = 0 [47,48]. The motivations behind the Gaussian form of p(a 2 ) mainly rely on the Central Limit Theorem. In fact, a 2 (r 2 0 in the notation in [34]) is the cumulative squared displacement of the N mon particle that move [34].
Plugging Equation (7) into Equation (6) leads to the following generalized HW equation (GHW), τ 0 is a suitable constant. An analogous law is anticipated for the viscosity η, given the known near proportionality with τ α [3]. Equation (8) is the form of the master curve Equation (2) being adopted in our studies. Other variants useful in the comparison with numerical and experimental results are listed in Appendixs A and B. Obviously, if the distribution p(a 2 ) is narrow and centred at a 2 0 , Equation (8) must reduce to the expression derived by Hall and Wolynes, Equation (4), τ (HW) α (a 2 0 ). For the specific choice of p(a 2 ), given by Equation (10), Equation (8) shows that this happens if σ 2 a 2 /8 u 2 2 a 2 /2 u 2 , namely, the ratio R defined as is vanishingly small. Equation (9) depends on the magnitude of DW factor so that, being the parameters σ 2 a 2 and a 2 independent of the physical state, the presence of the distribution p(a 2 ) is negligible when the DW factor is large, thus leading to a very narrow distribution of relaxation times, a characteristic of homogeneous dynamics. This suggests to read the condition R = 1 as the crossover between homogeneous and heterogeneous dynamics, i.e.,

R
1 homogeneous dynamics R 1 heterogeneous dynamics (10) Finally, we notice that the distribution p(a 2 ) in Equation (7) with a 2 min = 0 may be recast via Equation (4), as a log-normal distribution of relaxation times p(ln τ) where B is the normalization ensuring p(ln τ) d ln τ = 1, τ = τ 0 exp(a 2 /2 u 2 ). An interesting feature of p(ln τ) is that its width ∼ σ a 2 / u 2 increases by decreasing the DW factor.

Diffusion Coefficient
The diffusion coefficient D may be expressed via the above model by the relation [30] The above equation assumes that displacements as large as a occurring in a time τ (HW) α (a 2 ) are statistically independent. Notice that, although Equation (6) signals that the structural relaxation time is affected by the larger a 2 values, i.e., the longest relaxation times of the distribution p(ln τ), the diffusivity, according to Equation (12), is influenced by the shorter ones.
The explicit expression of the diffusion coefficient and an approximated version are given in Appendix A.2.

Stokes-Einstein Product
The Stokes-Einstein (SE) relation, Equation (1), states that the quantity Dη/T is constant if the diffusing particle changes neither the size nor the boundary conditions with the liquid. As the numerical evaluation of the viscosity is a delicate point, proxies are often used [27,65]. As an example, as η ∝ Tτ α in unentangled polymers [66], it is more suitable to study the breakdown of the SE law by considering the SE product where M is the number of monomers. K SE is expected to be independent of the chain length, as D ∝ 1/M in unentangled polymers [66] and the monomer relaxation at τ α poorly senses the chain connectivity. The above equation with M = 1 may be also used for liquids where the elementary units are atoms or small molecules, as the temperature factor in the ratio Dη/T provides a change of approximately ∼20% in fragile glass-formers [3], much less than the observed increase of K SE on approaching the glass transition [4,11]. The explicit expression of the SE product K SE derived within the vibrational dynamics model and an approximated versionK SE are given in Appendix A.3.

Transport and Relaxation in Polymeric Melts
The correlation between diffusivity, slow relaxation and fast vibrational dynamics has been studied by Molecular-Dynamics (MD) simulations of a coarse-grained model of a melt of linear unentangled polymer. Details about the model are given in Section 9. Even if rather crude, the model was proven to capture the universal aspects of the correlation and allowed an effective comparison with the experiment [47].
To provide a microscopic picture of the transport, the mean square displacement (MSD) of the monomer r 2 (t) is usually considered: where x i (t) is the position of the i-th monomer at time t. In addition to MSD, with the purpose of characterizing the relaxation, the self part of the intermediate scattering function (ISF) is also considered [26]: In an isotropic liquid, ISF depends only on the modulus of the wavevector q = ||q|| and features the rearrangements of the spatial structure of the fluid over the length scale ∼2π/q, leading to a decaying profile in time starting from F s (q, 0) = 1. In our case, ISF was evaluated at q = q max , the maximum of the static structure factor (7.13 ≤ q max ≤ 7.55) corresponding to the length scale of the monomer size. F s (q max , t) vanishes when the monomer displacement in a time t largely exceeds the monomer diameter. The time needed to make F s (q max , t) small is a measure of the escape time of the monomer from the cage formed by the neighbours, also known as the structural relaxation time τ α , customarily defined by the relation F s (q max , τ α ) = e −1 . Figure 2 shows typical MSD and ISF curves of the polymeric monomers. At very short times (ballistic regime), MSD increases according to r 2 (t) ∼ = (3k B T/m)t 2 and ISF starts to decay. The repeated collisions with the other monomers slow the displacement of the tagged one, as evinced by the knee of MSD at t ∼ √ 12/Ω 0 ∼ 0.17, where Ω 0 is an effective collision frequency, i.e., it is the mean small oscillation frequency of the monomer in the potential well produced by the surrounding ones kept at their equilibrium positions [64,67]. At later times, a quasi-plateau region, also found in ISF, occurs when the temperature is lowered and/or the density increased. This signals the increased caging of the particle. Trapping is terminated after an average time τ α . For t τ α , MSD increases more steeply. The monomers of short chains (M 3) undergo diffusive motion r 2 (t) ∝ t δ with δ = 1. For longer chains, owing to the increased connectivity, the onset of the diffusion is preceded by a subdiffusive region (δ < 1, Rouse regime) [68]. At long time, the monomer displaces in a diffusive way with diffusion coefficient D = lim t→∞ r 2 (t) /6t.

Vibrational Caged Dynamics and Debye-Waller factor
In our model polymer, the term "vibrational dynamics" refers to the rattling of the trapped monomer within the cage formed by the closest monomers. It is crucial to provide a robust criterion to assess the presence of the cage, which is anticipated to lack in liquids with high molecular mobility and fast relaxation. Compelling evidence of the cage effect is provided by the time velocity correlation function, which, after a first large drop due to pair collisions, reverses the sign since the monomer rebounds from the cage wall [64]. As an alternative route to reveal the cage effect, we consider the slope of MSD in the log-log plot Representative plots of ∆(t) for the polymer system are given in the inset of Figure 2a. ∆(t) tends to 2 at short times, due to the ballistic motion, and reaches the plateau level 1 at long times, owing to the diffusive motion. In the absence of caging effect, ∆(t) decreases in a monotonous way on increasing time. Caging is indicated by the presence of a minimum of ∆(t) occurring, irrespective of the physical state in the present model polymer, at t = 1.0(4). In actual time units, t is~1-10 ps [69].
The presence of the minimum paves the way to a robust definition of the DW factor u 2 , the mean square rattling amplitude of the monomer during the trapping period. In fact, the minimum, corresponding to the inflection point in the log-log plot of r 2 (t) ), separates two regimes. At short times, t < t , the inertial effects dominate, whereas for t > t , early escapes from the cage become apparent. Therefore, a convenient definition of the DW factor as a mean localization length is just MSD at t :

Debye-Waller Scaling of the Slow Relaxation
The monomer dynamics depends in a complex way on the state parameters. Nonetheless, there is clear correlation between the DW factors and the long-time relaxation dynamics. First examples are shown in Figure 2 by considering MSD and ISF. Note that states with equal DW factor have coincident time evolution of both MSD and ISF at least between t and τ α [49]. In Section 5.1, it will be shown that these results are a manifestation of Equation (3).
It is seen that the coincidence of the MSD curves is lacking at times longer than τ α for states corresponding to polymer chains with different length. This effect is not a failure of the scaling at times exceeding τ α , but a mere consequence of the complex dependence of MSD on the chain length since it is affected by all the Rouse modes [66]. In fact, if the correlation function of the single Rouse mode with the slowest relaxation time is singled out, i.e., the one with characteristic relaxation time given by the average chain reorientation time τ ee [66], the scaling is still observed after proper account of the chain length dependence, see Figure 3. The finding proves that Equation (3) holds also at a time τ ee being much longer than τ α .
As a side product of the coincidence of the ISF curves in states with equal DW factor seen in Figure 2, one has that states with equal DW factor u 2 have equal structural relaxation time τ α too. This can be reformulated via the master curve Equation (2), which, according to the model detailed in Section 2.1, takes the form given by Equation (8), i.e., a simple parabolic law between log τ α and 1/ u 2 [47,48]. Figure 4 tests Equation (8), written in the form given by Equation (A1) for a wide variety of physical states of our model polymeric melt [48]. It is also shown that the scaling holds if one considers the end-end chain reorientation time τ ee , i.e., the time needed by the correlation function C ee (t) to drop to e −1 , see Figure 3; although, in this case, it is described by a different master curve. The results prove that Equation (3) holds also at times τ ee much longer than τ α . Data from [48].

Signatures of the Heterogeneous Dynamics
MSD and ISF well-expose the cage effect, whereas the possible DH influence on their shape is less apparent. Figure 1 shows that DH is characterized by the presence of clusters of monomers with rather different mobility [11,12]. We now present and discuss two quantities well tailored to provide quantitative insight into this aspect.

van Hove Function
One central quantity of the DH analysis is the self part of the van Hove function G s (r, t) [26]: where x i (t) is the position of the i-th monomer at time t, and δ[·] is the three-dimensional Dirac delta function. In isotropic liquids, the van Hove function depends on the modulus r of r. The interpretation of G s (r, t) is direct. The product G s (r, t) · 4πr 2 is the probability that the monomer is at a distance between r and r + dr from the initial position after a time t. The moments of G s (r, t) are of interest: For n = 2, one recovers the usual mean square displacement (MSD). If the monomer displacement is a Gaussian random variable, G s (r, t) reduces to the Gaussian form [26]: Equation (20) is the correct limit of G s (r, t) at very short (ballistic regime, r 2 (t) = 3k B T/µt 2 ) and very long times ( diffusion regime, r 2 (t) = 6Dt, where D is the monomer diffusion coefficient).
The spatial Fourier transform of the self part of the van Hove function yields the ISF function, Equation (15) [26]. Figure 5a presents the self-part of the van Hove function G s (r, t), evaluated at τ α for the set of states with different mobility and relaxation shown in Figure 2. It is seen that if the relaxation and the mobility are fast, the shape of G s (r, τ α ) decreases by increasing the displacement r from the initial position. On the other hand, the states belonging to the D and E set, the ones with slowest relaxation, exhibit a tendency toward a bimodal pattern, namely, in addiction to particles undergoing small displacements, a shoulder at r ∼ 1 (the monomer diameter) is observed. This signals the presence of particles exhibiting fast displacements by solid-state jump dynamics [27]. Said otherwise, the quasi-bimodal pattern of the van Hove function is clear signature of DH. Four other aspects are to be noted: • The self-part of the van Hove function is expressed by suitable correlation functions, see Appendix B. Then, the coincidence of G s (r, τ α ) in states with equal DW factor observed in Figure 5a (the sets of states labelled as A, . . . , E) is in harmony with Equation (3). • Equation (3) also holds if one inspects the spatial dependence of the correlation function, e.g., the van Hove function, at τ α . In particular, even in the presence of DH. • Given their relation with G s (r, t), the coincidence of both MSD and ISF observed in Figure 2 for the sets of states labelled as A, . . . , E is strictly linked to the one observed in Figure 5a. • The pattern of the D and E sets of states is not consistent with the Gaussian limit G g s (r, τ α ), Equation (20), predicting a progressive decay with r, i.e., the DH dynamics is not Gaussian;  (21), of the states of Figure 2. On increasing the structural relaxation time from A states to E states, the system tends to increase the fractions of monomers with either much lower or much higher mobility with respect to the fraction predicted by the Gaussian approximation. Data from [53].
To quantify the deviations of the self-part of the van Hove function G s (r, τ α ) from the Gaussian limit, one defines the quantity [27,70] Figure 5b plots the ratio N s (r, τ α ). It exhibits increasing positive deviations at both short and large r values, evidencing the excess of nearly immobile and highly mobile monomers with respect to purely Gaussian behaviour, respectively. The analysis, in terms of the ratio N s (r, τ α ), reveals the wide distribution of mobilities pictured in Figure 1, right.

Non-Gaussian Parameter
An effective quantity to expose the time evolution of the non-Gaussian character of DH dynamics is the non-Gaussian parameter (NGP) [26]: where r 2 (t) and r 4 (t) are the mean square and quartic displacements of the particle at time t, respectively. α 2 vanishes if the displacement is Gaussian, i.e., it follows from a series of independent elementary steps with finite mean and variance. Figure 6 plots the NGP time evolution, Equation (22), for the set of states A, · · · , E and additional states with very slow relaxation. It is seen that NGP vanishes at very short times, as the ballistic regime is Gaussian in nature. At intermediate times, a peak value α max 2 is observed increasing with the relaxation times [46,71,72]. The maximum occurs at a time slightly shorter than the structural relaxation time τ α , as in simpler molecular systems [27]. A snapshot of the microscopic mobilities in a lapse of time τ α , where DH is quite apparent, is plotted in Figure 1 (right). At later times, NGP decreases as the monomer dynamics enters the homogeneous diffusive regime, which is a Gaussian process [25]. It is seen that states belonging to the same set A, · · · , E, i.e., with equal DW factor, have identical NGP in the time window [t , τ α ] at least. This agrees with Equation (3), given the relation of NGP with the moments of the self part of the van Hove function G s (r, t), Equation (19), and the expression of the latter in terms of suitable correlation functions, see Appendix B [73]. Note also the exponential increase of α max 2 with the ratio R defined in Equation (9) [47,48]. This is in harmony with the inequalities in Equation (10), stating that DH is characterized by R > 1.

SE Breakdown in Unentangled Polymers
The SE law is usually derived by considering the diffusivity of macroscopic bodies displacing in homogeneous viscous liquids [25]. The diffusion in the presence of strong DH does not comply with the SE law [15]. We have studied the SE breakdown in melts of unentangled polymers [29]. In these systems, helpful features are found [66]: (i) the diffusion coefficient D is inversely proportional to the chain length M, and (ii) the viscosity η is proportional to the end-end reorientation time which, in turn, is proportional to the structural relaxation time, e.g., see Figure 3, showing that states with equal structural relaxation time also have equal end-end reorientation time. Then, as discussed in Section 2.3, the study of the validity of the SE law is more efficiently carried out in terms of the product DMτ α , which is anticipated to be state-independent if the SE law holds. Figure 7 shows that in states with homogeneous Gaussian dynamics, i.e., with small α max 2 values, the R values are comparable or less than the unit value and the product D M τ α is nearly constant, i.e., the SE law holds true. On the other hand, in the presence of significant DH, i.e., α max 2 > α max 2,c = 0.40(5), one finds R > R c = 1.9(1) and the product D M τ α tends to increase, i.e., the SE law fails [11,12,27]. The comparison between α max 2 and R substantiates the conclusion that the magnitude of the ratio R allows one to conclude whether DH is appreciable or not, as suggested by Equation (10). As the ratio R-apart from constants-depends only on DW, see Equation (9), the finding supports previous conclusions that the long-time DH is rooted in the fast dynamics [44].   It is seen that states with equal R (Figure 7a), i.e., states with equal DW factor according to Equation (9), exhibit nearly equal values of the product D M τ α . A similar result has been reported for atomic binary mixtures [55] and metallic alloys [30]. Recognising that the diffusivity D and the structural relaxation time τ α reflect the long time behaviour of MSD and ISF, respectively, Figure 7a reveals that Equation (3) is valid even in the diffusive regime which is entered in polymer melts at times fairly longer than τ ee , being τ ee τ α .

Quasi-Universal SE Breakdown of Fragile Glass-Formers
Having noted that the SE failure is tracked by the DW factor in unentangled polymers, we now pose the question if this finding exhibits universal features. To this aim, we consider the ratio K SE /K 0 with K SE defined in Equation (13) and K 0 a scaling factor to ensure the unit limit value at large DW factor.
In Figure 8, we plot the ratio K SE /K 0 as a function of u 2 /u 2 g . We complement the MD results on unentangled polymers already presented in Figure 7 with other MD data, considering the diffusion of Cu and Zr atoms in metallic alloy, A and B atoms in a Lennard-Jones binary mixture [30] together with experimental data concerning the popular fragile glass-former ortho-terphenyl (OTP) [74,75]. Figure 8 evidences the good collapse of the SE violation in terms of the reduced DW. . Stokes-Einstein product K SE , normalised by its high temperature value K 0 (τ α 1 ps), as a function of the reduced DW factor u 2 /u 2 g , u 2 g being the DW factor at the glass transition. In addition to unentangled polymers, the plot also considers MD data concerning atomic binary mixtures (atoms labelled as A and B) and metallic alloys made by Cu and Zr atoms, as well as experimental data for ortho-terphenyl (OTP) [74,75]. Two predictions of the master curve are presented in terms of the quantityK SE andK SE , Equations (A8) and (A10), respectively. Both quantities have no adjustable parameters.K 0 andK 0 are suitable constants to ensure the unit limit value at large u 2 /u 2 g . A third master curve, drawn from the fractional SE law τ 1−κ α with κ = 0.85 (orange curve), is superimposed to the other curves. For numerical data, u 2 g is obtained according to the procedure outlined in [47]. Data from [30]. Figure 8 offers the opportunity to test the master curve predicted by the model of Section 2 with no adjustable parametersK SE , Equation (A8), and its approximationK SE , Equation (A10). It is seen thatK SE predicts a stronger SE breakdown than actually observed. Larger deviations are exhibited by the approximantK SE . How to improve the agreement: The expression of the diffusion coefficient in Section 2.2, assuming that displacements as large as a are statistically independent, aims at a SE productK SE , Equation (A8), with no additional adjustable parameters with respect to the ones of τ α , i.e., the ones of Equation (A2). This puts severe constraints on the shape of the distribution of the square displacements needed to overcome the relevant energy barriers p(a 2 ), Equation (7). The form of the distribution is adequate for large displacements to reach the transition state governing τ α , as proven by the effective fit of the MD data by the predicted master curve shown in Figure 4. However, the findings of Figure 8 suggest that it must be improved for small displacements affecting D. Alternatively, we may also state that the distribution p(ln τ), Equation (11), should be refined as far as the short relaxation times are concerned. Figure 8 shows that better agreement occurs by assuming the fractional SE form Dτ α τ 1−κ α [15,31] with τ α as given from Equation (A2). The best fit is found with κ = 0.85, which equals the universal value found by Mallamace et al. [76], deviating from the prediction of the "obstruction model" κ = 2/3 [15].

Displacement Correlation Length
Several results of the present paper suggest that the vibrational dynamics, as sensed by DW factor, provides insight into DH. A line of attack to understand how vibrational dynamics is related to slow relaxation is provided by the model of Section 2. The model is based on the distribution of the (squared) displacements needed by a particle to rearrange in the different local environments, Equation (7), leading in turn to the distribution of relaxation times, Equation (11). As noted in Section 6.2, the model needs further development. A further aspect to be improved is the absence of any detail on the localization of the particles with a given dynamics. This prevents any prediction concerning a peculiar aspect of DH, i.e., the existence of spatial domains with characteristic length scales where the particles undergo correlated motion, e.g., see Figure 1.
To make progress, it is worthwhile to preliminarily judge whether DW exhibits some correlation with possible dynamic length scales. To pursue this task, we studied the following monomer displacement-displacement correlation (DDC) functions [53,54]: An average over all the i-th and j-th monomers spaced by r is understood.û k (t 0 , t) is the versor of the displacement vector of k-th monomer in a time interval from t 0 to t 0 + t, u k (t 0 , t) = r k (t 0 + t) − r k (t 0 ) and δu k (t 0 , t) = |u k (t 0 , t)| − |u(t 0 , t)| , where |u k (t 0 , t)| is the modulus of the displacement. Henceforth, C δu (r, τ α ) and C u (r, τ α ) will be referred to as modulus (or mobility) and direction DDC functions, respectively. Local anisotropies and collective elastic solid-like response to the rattling of the monomers in the cage of their neighbours play a central role in the DDC build-up [64].
We consider DDCs of the states presented in part of the states in Figure 2. We remind that the states (i) exhibit different DH degree, e.g., see Figures 5 and 6,and (ii) are grouped in sets labelled B through E, each set being characterized by a single value of the DW factor. Figure 9a,b shows the spatial dependence of the direction and the modulus DDC functions, respectively, for the sets of states labelled B through E in Figure 2. Both correlation functions manifest damped oscillations in-phase with the pair correlation function g(r), thus evidencing that the correlated motion of a tagged monomer and its surroundings is influenced by the structure of the latter. This agrees with previous work on DDCs in Lennard-Jones systems [72,77], hard-sphere [78] and experiments on colloids [79]. The highest correlations are reached at a distance corresponding to the bond length b = 0.97 which demonstrates the highly concerted dynamics of bonded monomers. The correlation peaks, located at the first-, second-,... neighbours shells, vanish approximately in an exponential way on increasing the distance from the tagged particle (see insets of Figure 9). In more detail, Figure 9a shows that the direction correlations do not show significant increase in their spatial extension on increasing the structural relaxation time. Figure 9b shows the modulus (mobility) correlations. Differently from the direction correlations, their spatial extension increase meaningfully with the structural relaxation time (see also the inset of Figure 9b). Figure 9b clearly shows that physical states with equal DW factor, i.e., belonging to the same set of states (B, . . . , E), exhibit the same spatial correlations. This provides further support that Equation (3) also holds if the spatial dependence of the correlation function is considered for a given time up to τ α at least. To provide additional insight, we evaluated the length scales of the exponential decays of the DDC maxima with the distance ∼ exp[−r/ξ X (τ α )] with X = u, δu, thus defining two distinct dynamic correlation lengths pertaining to direction and modulus DDCs, ξ u (τ α ) and ξ δu (τ α ), respectively. Figure 10 shows these quantities. It is seen that the spatial extension of the modulus DDC increases quite a lot with τ α and reaches distances beyond the next-nearest neighbours for the states with the slowest relaxation. Instead, the direction correlations are virtually independent of the structural relaxation. Irrespective of the correlation length under consideration, Figure 10 also shows that they are equal within the errors for states with equal DW factor, i.e., belonging to the same set of states (B, . . . , E). Note that g(r) is virtually state-independent. The insets are semi-log plots of the corresponding main panels. Note the approximate exponential decay of the peak amplitudes with slopes ξ u (τ α ) and ξ δu (τ α ), respectively. Data from [53].
We are now in a position to compare our results with previous work on DDCs. Simulations of Lennard-Jones binary mixture (BM) observed that at time t α , corresponding to maximum dynamic heterogeneity, ξ BM δu (t α ) increases as the temperature decreases, whereas ξ BM u (t α ) is almost constant [80]. This agrees with our findings in Figure 10 concerning unentangled polymers. As to the modulus DDC correlation length, one finds [53] that after suitable algebraic manipulation to allow comparison [79], our changes of ξ δu (τ α ) with τ α are in quantitative agreement with the results of Bennemann et al. reported in a study of the same polymer system investigated here [72]. States with equal DW factor, i.e., belonging to the same set B, · · · , E exhibit equal directional and mobility correlation lengths. Data from [53].

Discussion
There is wide experimental, numerical and theoretical evidence that the fast vibrational dynamics, as sensed by the Debye-Waller factor u 2 , and the time scale τ α of the slow microscopic reorganisation of a liquid close to the transition to the glassy state are correlated in an universal way. Potential applicative implications concerning the quick characterisation of the stability of disordered structures with ultraslow relaxation are apparent. Less attention has been paid to a series of numerical MD simulation studies concluding in favour of strong correlations also between the vibrational dynamics and the dynamical heterogeneity, the spatial distribution of long-time mobility developing when approaching the disordered solid state. We reviewed these studies, mainly concerning melts of unentangled linear polymers, unifying all the results for the first time in terms of Equation (3). The latter has been tested both in space and time. In particular, we considered time-dependent quantities accounting for transport and relaxation like MSD, ISF and NGP and showed that they are related to the self-part of the van Hove function, which reduces to suitable correlation functions. In this respect, the correlation between the Debye-Waller factor and the breakdown of the SE law, a hallmark of DH presence, is seen as an ancillary consequence of the extension of Equation (3), which at times is much longer than τ α , where the motion is diffusive. We also inspected Equation (3) in space by considering both the van Hove function and DDC functions. Notably, DDC functions are collective in nature, differently from the self-part of the van Hove function, which is a single-particle quantity. This suggests that Equation (3) also holds for collective correlation functions. A further validation of this conclusion is offered by the collective stress-stress correlation function, which has been presented elsewhere [32] and not discussed in this review.
The understanding of the microscopic origin of the correlation between the vibrational dynamics and the heterogeneous dynamics close to the glass transition is still unsatisfactory in several respects. In particular, both the model discussed here, as well as other ones reported in the literature, even if successful in relating the Debye-Waller factor u 2 with the time scale τ α , are currently unable to account for the fact evidenced by the numerical simulations that the vibrational dynamics conveys also information on the spatial correlations between the mobility of different particles.

Methods
Most results discussed in this review concern a coarse-grained model of a linear polymer chain with M monomers is adopted. Bending and torsional potentials are neglected, i.e., the chain is fully flexible. While addressing the interested reader to the referenced papers for further details, we provide here some general aspect of the numerical model. We considered systems with total number of monomer N ≥ 2000. Non-bonded monomers at a distance r interact via the truncated parametric potential: where σ * = 2 1/6 σ and the value of the constant U cut are chosen to ensure U p,q (r) = 0 at r ≥ r c = 2.5σ. The minimum of the potential U p,q (r) is at r = σ * , with a constant depth U(r = σ * ) = .
Note that U q,p (r) = U p,q (r). Bonded monomers interact with a potential which is the sum of the Finitely Extendible Nonlinear Elastic (FENE) potential and the Lennard-Jones (LJ) potential [71]. The resulting bond length is b = 0.97σ, within a few percent. We set σ = 1 and = 1. The time unit is τ MD = (mσ 2 / ) 1/2 , with m being the mass of the monomer. Temperature is in units of /k B , where k B is the Boltzmann constant. We set m = k B = 1. All the data presented in this work are expressed in reduced MD units. It is interesting to map the reduced MD units to real physical units. The procedure involves the comparison of the experiment with simulations and provide the basic length (σ), temperature ( /k B ) and time (τ MD ) units [69,71,[81][82][83]. For polyethylene and polystyrene, it was found σ = 5.3 Å, /k B = 443 K, τ MD = 1.8 ps and σ = 9.7 Å, /k B = 490 K, τ MD = 9 ps, respectively [69]. For poly(vinyl alcohol) σ = 5.2 Å, /k B = 550 K and τ MD = 1.63 ps [83].
For polyisoprene σ = 6.7 Å, /k B = 307 K and τ MD = 10 ps [81]. The densities used in this and other studies are lower than the densities at atmospheric pressure, e.g., when mapping our model to polyethylene and polystyrene, one finds ∼0.5 and ∼0.7 g/cm 3 , to be compared to the actual values 0.78 and 0.92 g/cm 3 , respectively [69].
Author Contributions: F.P., A.T. and D.L. wrote the manuscript together. All authors have read and approved the final version.
Funding: This research was funded by the project PRA-2018-34 ("ANISE") from the University of Pisa.
Acknowledgments: A generous grant of computing time from IT Center, University of Pisa, and Dell EMC R Italia is gratefully acknowledged.

Conflicts of Interest:
The authors declare no conflicts of interest.

Abbreviations
The following abbreviations are used in this manuscript.
As τ (HW) α (a 2 ) depends on the parameter a in a much more marked way than a 2 , see Equation (4), an effective approximation of the diffusivity, as expressed by Equation (12) where a 0 is a constant.

Appendix B
The van Hove function for a uniform fluid is defined as Physically, Gdr is the probability of finding a particle j in a region dr around a point r at time t if the particle i was at the origin at time 0. We may recast the van Hove function by resorting to the time-dependent, microscopic particle density [26] ρ(r, t) = We rewrite Equation (A11) as where ρ is the average number density. Equation (A15) shows that the van Hove function is proportional to the density correlation function. It is easily shown that the van Hove function may be written as [26] G(r, t) = G s (r, t) + G d (r, t) where G s (r, t) and G d (r, t) are usually called the "self" and "distinct" parts. Equation (18) provides the explicit expression of G s (r, t). The distinct part is written as Finally, we show that both G s (r, t) and G d (r, t) may be expressed in terms of suitable correlation functions. To this aim, we define the auxiliary function B i (r, t) ≡ δ[r − x i (t)]. By repeating the same passages leading from Equation (A11) to Equation (A13), one has The last passage follows from the uniformity of the fluid. Analogously, one finds