Relation between dynamic heterogeneities observed in scattering experiments and four-body correlations

Dynamic heterogeneity is expected to be a key concept for understanding the origin of slow dynamics near the glass transition. In previous studies, quantitative evaluations of dynamic heterogeneity have been attempted using two different routes, i.e., the speckle patterns in scattering experiments or the four-body correlation functions of microscopic configuration data obtained from molecular dynamics simulations or real-space observations using confocal microscopy. However, the physical relationship between these dynamic heterogeneities obtained using different methods has not been clarified. This study proposes a connection between dynamic heterogeneities characterized based on speckle patterns and those obtained from four-body correlations. The validity of the relationship is also clarified.

The glass transition is a common phenomenon observed in metals, polymers, molecular and ionic liquids, and colloidal dispersions [1,2]. However, this phenomenon has remained a mystery in many areas for decades, such as the dramatically slow dynamics near the glass transition upon cooling; consequently, it is recognized as an open problem in physics. As an example of this problem, a dramatic increase in the viscosity of a liquid near its glass transition temperature is experimentally observed. One suggestion for the cause of this dynamic behavior is the cooperative movement of the molecules in a supercooled liquid. Dynamic heterogeneity (DH) has been proposed in previous studies [3,4] to understand the mechanism underlying these molecular dynamics. This concept characterizes the heterogeneity of molecules' mobility dependent on their positions in a liquid. The idea of DH is helpful for understanding the origin of the slow dynamics near the glass transition.
Two methods for evaluating DH have been reported to date to quantify the phenomenon of the cooperative movement of molecules; one is based on a four-body correlation function S 4 (Approach 1), and the other is based on the scattering intensity I of the speckle patterns measured in a scattering experiment (Approach 2).
In Approach 1, one calculates the function S 4 from microscopic configuration data generated in molecular dynamics (MD) simulations or obtained from real-space observations using confocal microscopy. Many studies on this topic have been reported since approximately 2000 [5][6][7][8][9][10][11][12][13][14][15], where several observables were used to measure the mobility of each particle to evaluate the DH. However, this method is not applicable in experiments because the information of each individual particle cannot be measured. Therefore, researchers have not yet clearly determined whether the results obtained using Approach 1 are valid in reality due to a lack of any method for confirmation. On the other hand, in Approach 2, one evaluates DH based on the scattering intensity I measured in a scattering experiment, such as X-ray photon correla-tion spectroscopy (XPCS). According to previous studies [4,[16][17][18], DH is quantitatively evaluated using the scattering intensity of speckle patterns measured in this type of experiment. Using this method, one can determine the actual extent of the heterogeneity of the dynamics. This approach has also been applied to computer simulations, such as MD simulations. However, to our knowledge, no reports have directly compared the two approaches yet after some earlier related studies, found for example in [19]. Based on the background described above, we have studied the relation between Approaches 1 and 2. The relation between these two approaches must be determined because they may focus on different physical origins. We applied the two approaches individually to the same system via MD simulations and evaluated the DH intensity, representing the amplitude of the spatio-temporal variations of slow/fast domains, using each approach independently to investigate this relation. We examined the dependence of the extent to which a state is supercooled on the intensity by setting different temperatures, representing states ranging from a normal liquid to a highly supercooled liquid. In this paper, we will present the two types of DH intensity calculated using Approaches 1 and 2 and discuss the relation between them by comparing the corresponding results. This paper is organized as described below. After our MD simulation model is briefly explained, we quantitatively show the methods we used to evaluate DH with this model. The results of the two evaluation approaches introduced above, i.e., by calculating the four-body correlation function and using the scattering intensity, are also presented. Finally, we discuss the relation between the two approaches by comparing their results and summarize our findings.
We have conducted MD simulations in three dimensions to investigate the origin of the slow dynamics in a supercooled state. Our simulation model is composed of two types of particles, 1 and 2, whose sizes and masses differ. The mass ratio between these two types of particles is set to m 2 /m 1 = 2, and the size ratio is set to σ 2 /σ 1 = 1.2, which is effective for preventing crystallization of the system at low temperatures [20]. In this model, a repulsive soft-sphere potential exists between particles: v αβ (r) = (σ αβ /r) 12 where r is the distance between two particles and σ α and σ β are the radii of particles α, β ∈ 1, 2, respectively. represents the strength of the pair-interaction, which is truncated at r = 3σ αβ . Spatial distance, time and temperature are measured in units of σ 1 , τ 0 = (m 1 σ 2 1 / ) 1/2 and /k B , respectively. We set the number of particles in the whole system to N = 1 × 10 5 , which includes equal numbers of particles of each type, i.e., N 1 = N 2 = 5×10 4 , and we fix the density to ρ = N/V = 0.8/σ 3 1 . Then, the system length is L = V 1/3 = 50.0 σ 1 . We set different temperatures of k B T / = 0.772, 0.473, 0.352, 0.306, and 0.267 to assess the dependence of DH on the degree of supercooling of the system. Note that the freezing point of the corresponding binary mixture is approximately T = 0.772 [20], below which the system is in a supercooled state.
In previous studies [5], the intensity of DH has been quantitatively evaluated by calculating the four-body correlation function S 4,k . In this approach, one obtains S 4,k from an order parameter Q k , which indicates the mobility (or immobility) of each particle in the system. Several forms of the order parameter have been suggested [5][6][7][8][9][10][11][12][13][14][15]. In this study, we define it as where r j is the position vector of particle j in real space and the bracket · · · denotes the ensemble average over all particles. The hatˆdenotes a value in Fourier space, and q is the wave vector corresponding to the position vector r. In this order parameter, the particle mobility deviation δD j is weighted by each particle position considering each volume ratio σ 3 j / σ 3 j . We quantified the mobility deviation by defining the immobility of particle as Here, ∆r j (t, τ ) is the displacement of particle j between time = t and t + τ . k represents the wave vector cor-responding to the particle displacement ∆r. The average · · · km is calculated for all wave vectors k that are consistent with the periodic boundary condition and satisfy k m − δk m ≤ |k| ≤ k m + δk m , with k m ≡ 2π and δk m = 0.001k m for Eq. (4) or 0.01k m for Eq. (8), assuming isotropy of the system. D j is defined as any quantity that represents the mobility (or immobility) of particle j from time t to t + τ . However, in this study, we use the definition in Eq. (4) for better consistency with the latter method that uses the scattering intensity. In this definition, the value of D j changes from 1 to 0 when the displacement of the particle |∆r j (t, τ )| becomes comparable to its radius in the time interval τ . The mean value averaged over time t is presented as · · · t . Therefore, δD j is the deviation of the mobility of particle j from the mean value. We provide the reader a more intuitive understanding of the order parameter Q k (r, t, τ ) by illustrating the spatial distributions of D j at three time intervals, τ /τ α = 10 −2 , 1, and 10, from a fixed initial time t in Fig. 1 for a normal liquid (k B T / = 0.772) and a highly supercooled liquid (k B T / = 0.267). Here the α relaxation time τ α is defined using the self-intermediate Note that D j is discrete in real space, and thus the plotted values are obtained through local spatial averaging. A darker (or lighter) color corresponds to the regions with less (or more) mobile particles; hence, an increase in the color variation indicates that more significant DH appears in the system. Accordingly, from these examples, DH becomes more intense at lower temperatures, i.e., when the system is highly supercooled. Based on the order parameter Q k , we define the fourbody correlation function as where · · · t,q is the average over time t and the angular components of q. The function S 4,k represents the spatial correlation of the mobility of each particle in the time interval τ . The correlation length ξ 4,k (τ ) and the intensity χ 4,k (τ ) are obtained by fitting S 4,k (q, τ ) to the Orstein-Zernike (OZ) form at small wavenumbers q [8,21]. The values ξ 4,k and χ 4,k thus represent the characteristic size and the intensity of DH, respectively. We have quantified χ 4,k by calculating the four-body correlation function S 4,k using the method described in the previous paragraphs, and we plot the results against the time interval τ in Fig. 2 (a) at different temperatures k B T / = 0.267, 0.306, 0.352, 0.473, and 0.772. This figure shows that the DH intensity χ 4,k reaches a peak χ * 4,k , which increases in height with decreasing The spatial distributions of Dj are drawn on an xy cross-sectional plane of the 3D system as a function of the time interval τ from a fixed initial time t. The upper panels correspond to kBT / = 0.772, and the lower panels correspond to kBT / = 0.267. According to Eq. (4), the value of Dj changes from 1 to 0 as the particle displacements increase with time interval τ (increasing from left to right). Although Dj uniformly is assigned a value of 1 (black) at τ = 0 or 0 (white) for τ τα, notable heterogeneity appears at intermediate time intervals τ τα, where darker (or lighter) colors represent the heterogeneous domains in which the dynamics of the particles are faster (or slower) than the average. Clearly, the variance (intensity) of the heterogeneity is larger at kBT / = 0.267 than at 0.772 around τ τα. temperature T . Therefore, the dynamics of the system become more heterogeneous when it is more supercooled. In addition, we show the dependence of the particle immobility averaged over all particles and time, D t = 1 N N j=1 D j t , which is identical to the selfintermediate scattering function F s (k m , τ ), in Fig.2 (b). As the temperature decreases, F s (k m , τ ) exhibits drastic slowing with a more stretched form, consistent with previous studies [21,22]. In the same figure, the standard deviation σ D = D 2 t − D 2 t is shown with the light colored areas. The relaxation time τ α clearly corresponds to the time interval when the DH intensity peaks, τ * 4,k . In previous experimental studies [4,[16][17][18], the DH intensity was reportedly quantified based on speckle patterns measured in scattering experiments, such as XPCS. The method discussed in the previous paragraph cannot be applied to these experiments in practice; therefore, DH must be quantified from speckle patterns as an alternative. While the scattering intensity I is feasible to measure experimentally, we calculate it in our MD simulations using the following equation: As the temperature decreases, the peak intensity, χ * 4,k , becomes more prominent, and the time interval when the peak occurs, τ * 4,k , is prolonged. (b) The particle immobility averaged over all particles in the system and time t, D t = Fs(km, τ ). Ten times the standard deviation at each time interval ±10σD is shown with the light colored areas around the average immobility D t.
where ρ k (t) is the Fourier transform of the particle density ρ(t) = N j=1 σ 3 j δ [r − r j (t)] and k is a wave vector. Fig. 3 (a) illustrates the scattering intensity in a speckle pattern calculated using Eq. (7). The pattern is plotted on a 2D plane with k z = 0 at k B T / = 0.267. The scattering intensity I varies with the vector k, and it reaches a sharp peak with a magnitude of approximately |k| = k m = 2π. One can confirm this peak in Fig. 3 (b), where the intensity averaged over the angular components of k, I(k, t) k , is plotted.
Using the scattering intensity I(k, t), we calculate the correlation function between the speckle patterns at two different times using the following equation [23,24]: This function C I represents the relaxation process of the scattering intensity at time t. In the presence of DH, C I fluctuates with time t, and thus the DH intensity is quantified by calculating the variance of this relaxation function [17,18]. Then, the normalized variance is defined as follows: The value χ k (τ ) represents the intensity of DH at time interval τ and indicates the extent to which the particle dynamics vary locally. We have quantified the DH intensity χ k (τ ) based on the scattering intensity I(k, t) as described above, and the results are shown in Fig. 4 (a). We have investigated the values at different temperatures k B T / = 0.267, 0.306, 0.352, 0.473, and 0.772. Note that the value χ k directly calculated from the equation contains some statistical noise due to the limited number of sampling points on the speckle pattern, n k . We have applied a correlation procedure to the values by extrapolation to the case of 1/n k → 0 (n k → ∞) to overcome this problem of spatial resolution [17,18,25]. The values approximated using the method described above are plotted in the figure. χ k (τ ) reaches a peak at each temperature, and the height of this peak becomes more prominent in a more supercooled state. These results are similar to those obtained using the approach based on the fourbody correlation function and in previous experimental studies [18,25,26]. We also show the two-time correla- As the temperature decreases, the intensity peak, χ * k , becomes more prominent, and the time interval when the peak occurs, τ * , becomes later. (b) The two-time correlation function averaged over time t with the same time interval τ , CI t − 1. Ten times the standard deviation at each time interval ±10σC I is shown with the light colored areas around the average immobility CI t − 1.
tion function calculated using Eq. (8) and averaged over time t, C I t − 1, in Fig. 4 (b). This averaged correlation represents the relaxation process of the scattering intensity I. Similar features are observed when comparing these results and the intermediate scattering functions shown in Fig. 2 (b). We also show the standard deviation σ C I = C 2 I t − C I 2 t with the light colored areas for each temperature. Note that the values are multiplied to make the standard deviations easier to visualize. As shown in this figure, DH becomes most significant around the relaxation time τ * of C I t − 1, which is almost equivalent to the α relaxation time τ α in Fig.=2. Additionally, the peak intensity χ * k = χ k (τ * ) increases when the system is in a more supercooled state, as also shown in Fig. 4 (a).
Finally, we discuss the relation between the two approaches for evaluating the DH presented above. Ac-cording to previous experimental studies [17,18,25], the temporal two-time correlation function at t, C I −1, is approximated as a stretched exponential exp [− [γ(t)τ ] µ(t) ], and the variance of C I (k m , t, τ ) arises from the fluctuations of the two dynamic parameters, γ(t) and µ(t). Additionally, the temporal self-intermediate scattering function D(k m , t, τ ) has also been approximated in a similar form [21]. Therefore, we assume that both the mean value and the variance of C I − 1 are approximately equal to those of D, i.e., Fig. 2 (b) and Fig. 4 (b), the present simulation data support the validity of these assumptions. From Eqs. (3), (4), (5) and (6), the DH intensity is approximated as From Eqs. (9), (10) and (11), we can relate the two types of DH intensity obtained from the four-body correlation function and from the scattering intensity as follows: As summarized in the Supplemental Material, the same result derived from Eq. (12) is also obtained using an equation based on the Siegert relation [27] with a slightly different choice for D j from Eq. (4).
We have checked the relation presented in Eq. (12) by comparing the results obtained using both approaches, as shown in Fig. 5. We show the dependence of χ 4,k and 4χ k on the time interval τ , with different colors corresponding to different temperatures. Note that χ 4,k is illustrated using circles, and χ k is illustrated using crosses. As shown in Fig. 5, the heights of the peaks, χ * 4,k and 4χ * k , are approximately the same at each temperature. Therefore, we conclude that Eq. (12) is valid to some extent and that the two approaches discussed above focus on the same physical property of the system.
In conclusion, in previous studies, only one of the following two routes has been used to quantitatively characterize DH in any particular case: 1) using four-body correlation functions or 2) using the speckle patterns observed in scattering experiments. In the present study, we successfully computed the intensity of DH using both methods by analyzing the same simulation data and compared the results in detail as functions of the temperature and the separation time. We confirmed a high level of agreement between the two sets of results throughout the whole parameter range of the present MD simulations. The present findings provide strong evidence for physical consistency between the DHs characterized using routes 1) and 2). This work was supported by Grants-in-Aid for Scientific Research (JSPS KAKENHI) under grant nos. JP 20H00129 and 20H05619.  We also discuss the relation between the two approaches for evaluating DH using the Siegert relation [27]: where g 2 is the time correlation function of scattering intensity I and g 1 is the normalized full-intermediate scattering function. Considering and Eq.(13) leads to Eq. (16) is apparently different from the empirically supported assumption, which we introduced as Eq. (10), but still holds quite well as we will see in Fig.6. Here, we have further assumed the instantaneous relation between the two functions as follows: Therefore, On the other hand, if the particle immobility in Eq. (2) is defined slightly different from Eq. (4) as the DH intensity is obtained from the equation Note that if we accept a crude order estimation assuming no correlations exist between ∆r i and ∆r j for i = j and use Eqs. (18) and (20), the two types of DH intensity obtained from the four-body correlation function and from the scattering intensity are again related as follows: which is equivalent to Eq. (12).
To test the validity of the above discussion numerically, we have assessed whether the Siegert relation is valid for the current system. We show the comparison between the time correlation function of scattering intensity C I t − 1 and the squared intermediate scattering function |g 1 | 2 in Fig. 6. The two functions show good agreement with each other, and Eq. (16) works in the MD simulations. We also have assessed the relation presented in Eq. (22) by comparing the results obtained using both approaches, as shown in Fig. 7. We show the dependence of χ 4,k and 4χ k on the time interval τ , with different colors corresponding to different temperatures. Note that χ 4,k is illustrated using circles, and χ k is illustrated using crosses. As shown in Fig. 7, the heights of the peaks, χ * 4,k and 4χ * k , are approximately the same at each temperature. However, the gaps between the two DH intensities in Fig. 5 are smaller than those in Fig. 7, while the gaps in Fig. 7 can be reduced if an adjustable parameter is introduced in the crude estimation Eq. (21).