Identifying structural signature of dynamical heterogeneity via the local softness parameter

When a liquid is cooled fast, it evades crystallization, and below the onset temperature enters a supercooled liquid domain where the properties of the liquid appear quite different from that at high temperatures. Although the average structure across the onset temperature appears to change continuously, the dynamics shows a marked difference and below the onset temperature, it slows down substantially and also becomes heterogeneous ediger_spatial_heterogenity ; kob_plimpton_dynamic_heterogenity . The origin of this dynamic heterogeneity in supercooled liquids is a topic of intense research ediger_spatial_heterogenity ; kob_plimpton_dynamic_heterogenity ; harowell_nature_phy ; harrowell_prl_2006 ; cooper_free_volume ; liu_nature ; rottler ; cubuk . In analogy with crystals where under external perturbation regions with structural defects show a higher probability of rearrangements taylor_plastic_deformation ; collings_colloid_crystal , it is often suggested that in supercooled liquids such structural defects may also be present, giving rise to the dynamic heterogeneity. However, while identifying structural defects in crystals with otherwise well-defined structure is trivial, doing the same in a supercooled liquid where particles are arranged in a disordered manner is a non-trivial task jp_dyre_structure .

Harowell and coworkers investigated different properties of the initial configuration of a supercooled liquid and analyzed their correlation with irreversible rearrangements in the system harowell_nature_phy ; harrowell_prl_2006 ; cooper_free_volume . They found that the local inherent structure potential energy harrowell_prl_2006 and free volumecooper_free_volume does not have any correlation with these rearrangements, but the Debye-Waller Factorharrowell_prl_2006 and the low-frequency normal modes of the systems are spatially correlated with the irreversible structural rearrangementsharowell_nature_phy . Liu and Manning have shown that under shear, the rearrangements of particles take place in the regions which contribute to these low-frequency normal modesmanning_vib_modes . Smessaert and Rottler did a quantitative analysis of these low-frequency soft modes and showed that they are long-lived rottler . There are also studies where the dynamics is connected to the elastic properties of the system shoving ; local_elasticity_jean_louis ; lerne ; zaccone_pnas . The heterogeneity of the local elastic modulus was found to be correlated with the dynamic heterogeneity showing regions with low shear modulus have higher plastic activity local_elasticity_jean_louis . A recent study has proposed that the slowing down of the dynamics is controlled by the mesoscopic elastic stiffness parameter, which is more sensitive than the shear modulus lerne .

Apart from the above mentioned correlation between vibrational and elastic modes with the dynamics, there are also other studies that calculate properties that are purely structural in nature and correlate them with the dynamics. Using mutual information technique, it was shown that coarse-grained energy and density correlates with the dynamics paddy . Another technique that has been quite successful in the recent time is Machine learning (ML). ML techniques allow for the identification of a softness parameter that encodes in a single number a large number of structural descriptors cubuk ; liu_nature ; olivier_PRE . It was shown that below the onset temperature, the dynamics is controlled by this softness parameterliu_nature and this softness parameter for attractive and repulsive systems are different, which eventually leads to the difference in their dynamics olivier_PRE . It was also shown that this softness parameter could identify structural defects which are similar to dislocation in the crystalline solids, and ductile materials have a large number of soft spots compared to brittle materials richard_prm . A recent experimental study of colloidal glasses showed that this ML softness parameter can predict the devitrification process ganpati_NP_2021 . Unfortunately, this softness consists of a linear reweighting of the local pair correlation, the weights being blindly found by the ML algorithm and its physical connection to structure remains unclear.

Interestingly in the liquid state theories hansen like mode coupling theory gotze_condmat and density functional theory wolynes_kirkpatrick the structure plays an important role in determining the dynamics. However, in the supercooled liquid regime, since the change in the structure is small and gradual, whereas the dynamics slows down dramatically, the validity of these theories has been questioned. Also, studies showing that in this regime, systems with similar structures have orders of magnitude differences in their dynamics further supports the idea that in this regime structure does not play a role in the dynamics berthier_tarjus_prl_2009 ; bertheir_tarjus_EPJE_2011 . However certain extensions of these theories have been found to work reasonably well in the supercooled regime sarika_PNAS ; chong_pre ; role_pair_configuration ; manoj_PRL ; indranil ; Schweizer_2003_jcp ; schweizer_2005_jcp ; mei_pnas_2021 ; ghosh_schweizer_jcp2020 ; mirigian_schweizer_jcp2014 ; mei_schweizer . In a study involving some of us the dynamic density functional theory (DDFT) wolynes_kirkpatrick ; schweizer_2005_jcp was used to develop a microscopic mean field theory role_pair_correlation . It was shown that the softness of the mean field caging potential described by the structure of the liquid is connected to the dynamics for a wide variety of systems, even for those which are out of equilibrium manoj_PRL . However all these theories have their root in the liquid state theory and thus only deal with the average properties of the system. Also note that both softness and dynamics can have a local microscopic variation and it is well known that the correlation between average quantities does not guarantee a correlation at a microscopic level. Thus using liquid state theories to study local properties and a causal relationship between structure and dynamics at the local level is a nontrivial and challenging task.

In this paper, we study the correlation between softness and dynamics at a microscopic level. Using DDFT formalism we first describe the softness parameter at a local level in terms of the structure of the liquid. We then show that this microscopic softness can capture the structural heterogeneity, and the lifetime of the softness parameter is similar to the lifetime of the cage. This, we believe, is a nontrivial. The most important result is the observed causal relationship between local softness and local rearrangements.

The rest of the paper is organized as follows. In the next section, we provide the simulation details. In Sec. III, we present the formalism we use to identify local rearrangements. Sec. IV presents the details of the calculation of the softness parameter at the local level and looks at its distribution and time evolution. In Sec. V we study the correlation between the local rearrangement events and the local softness. Finally, in Sec. VI, we present the conclusion.

The system we studied is the Kob-Andersen model for glass-forming liquid, which is a binary mixture (80:20) of Lennard-Jones (LJ) particles PRL_73_1376_1994 . The interaction between the particles i and j, where i,j=A,B (the type of the particles) is given by

We use $\sigma_{AA}=1.0$, $\sigma_{AB}=0.8$, $\sigma_{BB}=0.88$, $\epsilon_{AA}=1.0$, $\epsilon_{AB}=1.5$, $\epsilon_{BB}=0.5$, $m_{A}=m_{B}=1$ and Boltzmann constant $k_{B}=1$. In our simulations we have used periodic boundary conditions and Nosé-Hoover thermostat with integration timestep 0.0025$\tau$. The time constants for Nosé-Hoover thermostat are taken to be 100 timesteps. The total number density $\rho=N/V=1.2$ is fixed, where V is the system volume, and N=4000 is the total number of particles.

In this work, we aim to correlate structure parameter with the mobility at a local level; for identifying fast particle or an event, there are many ways which can be used, such as doing iso-configurational runs and identify irreversible reorganisationsharowell_nature_phy or tracking the mean square displacement over a period of time. Here we have used a method which was first proposed by Candelier et al.candelier ; PRE_88_022314_2013 where they calculate a quantity $p_{hop}(i,t)$ which captures for every particle i in a certain time window W = $[t_{1},t_{2}]$, the cage jumps when the average position of the particle changes rapidly. The expression for $p_{hop}(i,t)$ is

for all t $\in$ W, where averages are performed over the time intervals surrounding time t, i.e., U = [$t-\Delta t/2$, t] and V = [t, $t+\Delta t/2$] where $\Delta t$ should be a timescale over which the particles can rearrange. For a time window W the small value of $p_{hop}$ means the particle is contained within same cage and conversely if $p_{hop}$ is large this means particle is within two distinct cage (see APPENDIX A for details). To identify a rearrangement event we set a threshold value $p_{c}$, as done in Ref.olivier_PRE where $p_{c}$ is chosen as the root mean squared displacement $<\Delta r(t)^{2}>$ value, computed at a time where the non-Gaussian parameter $\alpha_{2}$ = $\frac{3<\Delta r^{4}(t)>}{5<\Delta r^{2}(t)>^{2}}-1$ has a maximum. When $p_{hop}>p_{c}$ we consider that a rearrangement event has taken place. The $p_{c}$ values are 0.115 for T=0.47, 0.130 for T=0.53, 0.141 for T=0.58, 0.159 for T=0.7, and 0.178 for T=0.8. We also vary $\Delta t$ from $15$ to $75$ LJ units and find that qualitatively the results remain quite similar. For the rest of the work we consider $\Delta t=15$ LJ unit. We also find that the neighbourhood of the fast particles are structurally less ordered (see APPENDIX B for details).

The main objective of this work is to identify a parameter that can describe structural heterogeneity and understand its role in dynamical heterogeneity. We have already shown that the softness of the mean-field potential can explain structural heterogeneity. We also find that the lifetime of this softness parameter is comparable with the timescale of the system dynamics. We now put this parameter through a more stringent test and study the causal relationship between the softness and dynamics. The dynamics expressed via local rearrangement events is computed using a well-established method olivier_PRE ; candelier ; candelier_prl ; PRE_88_022314_2013 . We find that particles that undergo local rearrangements have a less structured neighbourhood (see APPENDIX B for details). We also find that the percentage of particles with a softness value greater than the peak, which undergoes rearrangement increases with a decrease in temperature. At T=0.7 it is 65%, at T=0.58 it is 69% and at T=0.47 it is 75%. (see Fig.4).

To quantify the correlation between softness value and a rearrangement event, we calculate the fraction of particles having a specific value of softness that undergoes rearrangement, $P_{R}(S)$ as a function of $S$ at different temperatures (Fig.5(a)). We find that at all temperatures, $P_{R}(S)$ has a softness dependence. Particles with higher values of softness have a higher rearranging probability.

This dependence becomes orders of magnitude stronger at lower temperatures. This shows that at lower temperatures, the dynamics get more coupled to the local structure. We have already demonstrated that with a decrease in temperature, the structural heterogeneity increases. We believe that these two factors, i) increase in structural heterogeneity ii) larger coupling between structure and dynamics, together lead to the well-known increase in the dynamic heterogeneity at lower temperatures. Particles with small softness values (confined in a deeper mean-field potential) are less mobile compared to particles that have larger values of softness (confined in a shallower mean-field potential). We now check if the probability of a local rearrangement can be expressed in an Arrhenius form $\it i.e.$ $P_{R}(S)=P_{0}(S)exp(-\Delta E(S)/T)$ where the activation energy is a function of softness. $log_{10}P_{R}(S)$ as a function of inverse temperature for different softness values show a linear behaviour (Fig.5(b)), thus confirming the Arrhenius form. The energy barrier, $\Delta E(S)$, the slopes of the lines are a function of $S$. $P_{R}/P_{0}$ against $\Delta E(S)/T$ (inset of Fig.5(b)), shows a data collapse, thus further validating the Arrhenius form. At smaller values of softness, $P_{R}(S)$ strongly depends on $T$. In Fig.5(b), we find that the extrapolated lines cross each other at a certain temperature, and above this temperature, it appears that $P_{R}(S)$ increases with a decrease in $S$, which is an unphysical result. Thus above this temperature, where the extrapolated lines cross the correlation between $P_{R}(S)$ and softness does not remain valid. Interestingly the extrapolated lines cross at $T=0.801$ which is similar to the onset temperature of the system atreyee_onset . Below we argue why it should be related to the onset of the glassy dynamics.

The soft spots are equivalent to the defects in the crystals. Thus they should come from a solid-like description of the system. When describing the softness, we assume a well-defined cage that remains static over the time period that we are calculating the dynamics. We have also shown that the lifetime of the softness has a correlation with the lifetime of the cages. It is well-known that in supercooled Kob-Andersen model liquid, these cages appear below the onset temperature, where the decoupling between the short and long time dynamics starts sastry_debenedetti ; kob_mctheory1995 . Thus our theoretical formalism rightly predicts that the correlation between softness and dynamics starts around the onset temperature, and at lower temperatures where there is a larger decoupling between the short and long time dynamics leading to a longer lifetime of the cages, the local softness becomes a better predictor of the local rearrangement events.

In Fig.5(c) we plot the activation energy, $\Delta E(S)$ as a function of S, and it shows that as softness decreases, the activation energy increases. The range of energy in our study (1.2-2.9) is narrower than that obtained in the ML study (5.0-11.0)liu_nature , but the energy values are similar. There are two possible reasons why our study predicts a narrower range. i) To calculate the activation energy for different softness values, we need to plot a temperature dependence of $P_{R}(S)$ which can be done in the range of softness where the $P(S)$ at different temperatures overlap. Since $P(S)$ in our study shows a larger shift with temperature, this overlapping range is narrower than the ML study. ii) Our calculation of softness is obtained from a mean-field microscopic theory, thus bound by the two-body microscopic correlation functions. It is possible that beyond what the theory predicts, other correlation functions contribute to the softness parameter and are picked up in the weight function of the ML study. However, as discussed in the ML study, our study reveals that the dominant contribution comes from the two-body terms. In Fig.5(d), we plot the softness dependent overlap function, $q(S,t)$ at different softness values (highest, lowest, and moderate) and at two temperatures. The time scale of the overlap function shows a stronger softness dependence at lower temperatures. This clearly shows a causal relationship between the softness and the dynamics that manifests itself more at lower temperatures.

In this work, we describe the softness at the local level. We find that the softness parameter predicts the presence of structural heterogeneity, which grows at lower temperatures and is also longer lived. The lifetime of the softness parameter is correlated with the well-known cage structure in the supercooled liquids. We establish a causal relationship between the local softness and the local rearrangement events. We further show that even for particles that do not undergo rearrangement, the softness parameter evolves in time due to the dynamics in the neighbourhood, giving rise to the well-known facilitation effect chandler_pnas2003 ; chandler_2010annurev . An immobile hard region can eventually become soft and have a higher probability of becoming mobile due to rearrangements in the neighbourhood.

The results presented here are strikingly similar to that obtained in the ML work, although the methodology of obtaining the softness in the two different studies are entirely different liu_nature . As suggested by the Authors themselves in the ML study, since the softness depends on a large number of local parameters, it is difficult to interpret the physical meaning of the softness liu_improved_softness . The conjecture is that through the multiple local parameters, the information of the local cage around particles is incorporated in the softness; however, there is no direct proof of it. The advantage of the present study is it calculates the softness from a microscopic theory that describes the caging potential in terms of the structure of the liquid. The depth of the caging potential is similar to the local energy that is obtained by reweighting the local pair correlation by the direct correlation function. We find that there exists a causal relationship even between the depth of this local caging potential and the dynamics (see APPENDIX B for details). Given the direct correlation of the present softness parameter and the liquid structure, it will be interesting to understand its correlation with other existing definitions of the soft fields like the softness parameter obtained in machine learning studies, the quasi localized vibrational modes and the elastic stiffness parameter. These problems will be addressed in future studies.

Acknowledgement

S. M. B thanks Chandan Dasgupta, Olivier Dauchot and Jeppe Dyre for discussion, and acknowledges funding support from Science and Engineering Research Board (SERB), Department of Science and Technology Government of India. M. S. thanks Ujjwal Kumar Nandi, and Palak Patel for the help with the initial setup of the system and discussions and thanks SERB for fellowship.

APPENDIX A: Identifying rearrangements

In this work, we aim to correlate structure parameter with the mobility at a local level; for identifying fast particle or an event, there are many ways which can be used, such as doing iso-configurational runs and identify irreversible reorganisationsharowell_nature_phy or tracking the mean square displacement over a period of time. Here we have used a method which was first proposed by Candelier et al.candelier ; PRE_88_022314_2013 where they calculate a quantity $p_{hop}(i,t)$ which captures for every particle i in a certain time window W = $[t_{1},t_{2}]$, the cage jumps when the average position of the particle changes rapidly. The expression for $p_{hop}(i,t)$ is

for all t $\in$ W, where averages are performed over the time intervals surrounding time t, i.e., U = [$t-\Delta t/2$, t] and V = [t, $t+\Delta t/2$] where $\Delta t$ should be a timescale over which the particles can rearrange.

For a time window W the small value of $p_{hop}$ means the particle is contained within same cage and conversely if $p_{hop}$ is large this means particle is within two distinct cage. To identify a rearrangement event we set a threshold value $p_{c}$, as done in Ref.olivier_PRE where $p_{c}$ is chosen as the root mean squared displacement $<\Delta r(t)^{2}>$ value, computed at a time where the non-Gaussian parameter $\alpha_{2}$ = $\frac{3<\Delta r^{4}(t)>}{5<\Delta r^{2}(t)>^{2}}-1$ has a maximum. When $p_{hop}>p_{c}$ we consider that a rearrangement event has taken place. The $p_{c}$ values are 0.115 or T=0.47, 0.130 for T=0.53, 0.141 for T=0.58, 0.159 for T=0.7, and 0.178 for T=0.8. We also vary $\Delta t$ from $15$ to $75$ LJ units and find that qualitatively the results remain quite similar. For the rest of the work we consider $\Delta t=15$ LJ unit. In Fig.6 we plot the $p_{hop}$ for a particle which shows few rearrangement events.

In Fig.7 we plot the rdf of the fast A particles, $g_{AA}(r)$ and $g_{AB}(r)$. For comparison, we also plot the rdf of all A particles. We find that the neighbourhood of the fast particles are structurally less ordered.

APPENDIX B: Calculation of local rdf and its effect on the softness parameter

In the calculation of $\beta\Phi_{r}(\Delta r=0)$ there are two terms, one describing the interaction with other ‘A’ particles and is a function of $g_{AA}$ and the other describing the interaction with the ‘B’ particles and is a function of $g_{AB}$. We find that the correlation between the dynamics and softness is best described when we set an upper limit to the r integration such that they correspond to the respective minimum of $g_{AA}$ and $g_{AB}$. This implies that it is the first nearest neighbours which contribute most to the softness. Since we have a finite value of $\delta$, the approximate rdf at low ‘r’ has a finite value where the actual rdf goes to zero. The depth of the mean-field caging potential, $\Phi_{r}(\Delta r=0)$ is calculated as a product of $C(r)$ and g(r) and the $C(r)$ has terms proportional to the interaction potential, $U(r)$ which diverges at small ‘r’. Thus in this range, small errors in the rdf value get magnified in the calculation of $\Phi_{r}(\Delta r=0)$. To minimize this error, we also put a lower cutoff in the ‘r’ integration. This lower cutoff in ‘r’ range is chosen as the value where the average rdf calculated from the histogram vanishes.

As expected, the value of the softness parameter is dependent on the $\delta$ we use to describe the local rdf. However, we show in Fig.11 that the variation of $\delta$ just shifts the value of the softness parameter. Since our analysis depends more on the relative value of the softness parameter rather than the exact value, the physics is expected to be independent of $\delta$. In Fig.12 we plot the correlation of the softness of the particles calculated using different values of $\delta$. We find that although the value of the softness depends on $\delta$, the relative degree of softness of any particle is almost independent of $\delta$, and the softness values show a strong correlation.

In Fig.13 we plot the probability distribution of the local softness, $P(S)$ at T=0.47. In the same figure, we also plot the probability distribution of the softness of the particles just before they undergo a rearrangement, $P(S|R)$. This comparison is done for different $\delta$ values and the percentage of particles that undergo rearrangement and has a value of softness greater than the peak varies a little with $\delta$. We find that for $\delta$ = 0.02 it is 66%, $\delta$ = 0.04 it is 74%, $\delta$ = 0.09 it is 75.2%, and for $\delta$ = 0.12 it is 71%. This is similar to that observed in the Machine learning study where they found that the rdf provides 77% prediction accuracy of rearrangements liu_nature .

In Fig.14 we plot the fraction of particles that undergo rearrangement as a function of softness at few temperatures. This plot is done at different $\delta$ values. We find that irrespective of the $\delta$ value, the study shows that the probability of rearrangement depends on softness, and this dependence is stronger at lower temperatures. For the rest of the study, we assume $\delta=0.09$, which is similar to the value used in earlier studies piaggi ; paddy .

In Fig.15 we study the causal relationship between the depth of the local caging potential $\beta\Phi_{r}(\Delta r=0)$ and the local rearrangements. We find that 62$\%$ of the particles that undergo rearrangement has a depth of potential lower than the peak value of the distribution (Fig.15(a)). We also find that the probability of rearrangement is dependent on $\Phi_{r}(\Delta r=0)$ and this dependence becomes stronger at lower temperatures (Fig.15(b)). We also find that the dynamics can be expressed in an Arrhenius form where the barrier is dependent on $\Phi_{r}(\Delta r=0)$ (Fig.15(c)). The correlation between the local rearrangement and the depth of the local potential is present only below the temperature where the $log_{10}P_{R}$ vs $1/T$ plots cross each other. This temperature T=0.80 is similar to that predicted from the softness parameter and is connected to the onset temperature of the system atreyee_onset . The barrier for the Arrhenius dynamics is linearly proportional to the depth of the caging potential.

References