Is the glassy dynamics same in 2D as in 3D? The Adam Gibbs relation test

The binary mixture ($A_{80}B_{20}$) of Lennard Jones particles introduced by Kob and Andersen (KA) Kob1995 is a model glass forming liquid in three dimensions. In the case of the modified KA (MKA) model in two dimensions (2D), the composition is set to $A_{65}B_{35}$ Bruning2008 . The interaction potential, including correction terms that make both potential energy and force to go to zero smoothly at the cutoff, is given by:

where $\alpha,\beta\in\{A,B\}$. Units of length, energy are $\sigma_{AA}$, $\epsilon_{AA}$ respectively. $\epsilon_{AB}=1.5$, $\epsilon_{BB}=0.5$, $\sigma_{AB}=0.80$, $\sigma_{BB}=0.88$. The cutoff is at $2.5\,\sigma_{\alpha\beta}$. Simulations were done at a fixed number density $\rho=1.20\,\sigma_{AA}^{-3}$ with system size $N=2000$ for a temperature range ${k_{B}T\over\epsilon_{AA}}=0.36-2.00$. The MD integration time step was in the range $dt/(\sqrt{\sigma_{AA}^{2}/\epsilon_{AA}})=0.003-0.005$, with smaller $dt$ values used at higher temperatures. For this model the onset temperature is $T_{onset}=0.80$. It serves as the reference temperature to differentiate between the low and the high $T$ regimes. For $T>0.90$ (reduced units), only single trajectory has been chosen and for $T<0.90$, 4-5 independent trajectories were analyzed. The runlength for all trajectories analyzed is more than $100\,\tau_{\alpha}$.

It is a 50:50 binary mixture of purely repulsive soft disks in two dimensions Karmakar2010 , with a inverse power law exponent 10. The interaction potential is given by,

where $\alpha,\beta\in\{A,B\}$. $\epsilon=1.0$, $\sigma_{AA}=1.0$ and $\sqrt{\frac{\sigma_{AA}^{2}}{\epsilon_{AA}}}$ set the units of energy, length and time respectively. $\sigma_{BB}$ = 1.40, $\sigma_{AB}$ = 1.18 $\neq$ $\frac{\sigma_{AA}+\sigma_{BB}}{2}$, $c_{0}=-0.806140903539923$, $c_{2}=0.7$, $c_{4}=-0.156300219287607$ represent the correction to make the potential vanish at cutoff$=1.385418025\,\sigma_{\alpha\beta}$ continuously upto the second derivative. Simulations were done at a number density of $\rho\sigma_{AA}^{3}=0.85$ over a range of temperatures ${k_{B}T\over\epsilon_{AA}}=0.46-2.00$, for a system size of $N=2048$. The onset temperature is at $T_{onset}=0.80$ (reduced units) which demarcates the high and low temperature regimes. For ${k_{B}T\over\epsilon_{AA}}\geq 1.00$, a single trajectory has been chosen with integration time step $dt/(\sqrt{\sigma_{AA}^{2}/\epsilon_{AA}})=0.003$ and for $T\leq T_{onset}$ 3-5 independent trajectories were used with $dt/(\sqrt{\sigma_{AA}^{2}/\epsilon_{AA}})=0.005$. The runlength for all trajectories analyzed is more than $100\,\tau_{\alpha}$.

To compute the anharmonic correction by method 2 described in Sec. II, MD runs were performed for short time duration $\sim$ MSD plateau time. These runs were started from inherent structures obtained from various parent temperatures $T_{p}$. The range of $T_{p}$ values were chosen to be $T_{P}=0.36-0.90$ for the 2DMKA model and $T_{p}=0.46-0.80$ for the 2DR10 model respectively. Number of runs were 26-926 for the 2DMKA model and 1000-3985 for the 2DR10 model depending on the temperature. For the same initial configuration, runs were performed at various target temperatures $T$ from $10^{-6}$ upto $T_{p}$. At all target temperatures $T$, runlengths were chosen to be for 1000 $MD$ steps for the 2DMKA model and $3000$ MD steps for the 2DR10 model. Integration time step was set to $dt/(\sqrt{\sigma_{AA}^{2}/\epsilon_{AA}})=0.003$ for both models.

In Figs. 1 and 2 we compare the dynamics for the 2DMKA and the 2DR10 models respectively using two estimators: the intermediate scattering functions and the MSD. The estimators computed from standard displacement field are shown in the top rows and those from the cage-relative coordinates are displayed in bottom rows. Both standard and cage-relative estimators show qualitative signatures of glassy dynamics KobBinderBook . However, we note the quantitative differences: (a) the plateau regime is more prominent for cage-relative $F_{s}(K,t)$ and MSD, (b) plateau height $f_{c}$ of the cage-relative intermediate scattering function is higher and (c) the transient wiggle in the plateau regime of the cage-relative $F_{s}(K,t)$ due to finite size effect is absent in both the models. These observations are in sync with previous studies Shiba2016 .

Next we study the effect of long wavelength fluctuations on the slowdown of dynamics with reducing temperature. First we compute the usual $\tau_{\alpha}$ from the dynamics of the standard displacement field in Fig. 3 left columns. Note that Ref. Sengupta2012 used the overlap function while the present study analyzes the $F_{s}(K,t)$. We show for both the models that $\tau_{\alpha}$ extracted from either correlation function are mutually proportional. More interestingly, we extract $\tau_{\alpha CR}$ from the cage-relative intermediate scattering function. We see that the main effect of the long wavelength fluctuation is quantitative: the cage-relative displacement field has longer $\alpha$-relaxation timescale. This finding is consistent with previous observations Illing2017 . However, the $T$ dependence of the timescales are similar, as the VFT fit $\tau_{\alpha}(T)=\tau_{0}\,\exp({D\over T-T_{0}})$ to $\tau_{\alpha}$ from both the $F_{s}(K,t)$ and the CR-$F_{s}(K,t)$ yield similar divergence temperature $T_{0}$ and kinetic fragility ${T_{0}\over D}$, see Table 1. For the sake of completeness, we also show the high T Arrhenius fit to CR-$F_{s}(K,t)$ timescale to mark the onset temperature of non-Arrhenius dynamics.

In the low $T$ non-Arrhenius regime, the various characteristic timescales such as the $\alpha$-relaxation time and the diffusion coefficient becomes decoupled - a phenomenon usually labelled as the breakdown of the Stokes-Einstein relation (SEB). Consequently, in the context of testing the AG and the RFOT relations an interesting open question is the choice of timescale to characterize dynamics, which the theories does not explicitly specify. Some of us in a previous study 2017Parmar showed that for the AG relation, the diffusion coefficient is the more fundamental choice of timescale. Thus we also characterize dynamics by extracting the diffusion coefficient $D_{CR}$ from the cage-relative MSD, see middle panels of Fig. 3. In the right columns of Fig. 3, we compare the time scales $\tau_{\alpha CR}$ and the $D_{CR}$ and find a decoupling (SEB) at low temperatures. Thus taken together, these observations indicate that the long wavelength fluctuation affects the dynamics in 2D quantitatively, but qualitative features of the dynamics remains the same.

After accounting for the long wavelength fluctuation in dynamics, we pose the next natural question: how does the relation between dynamics and thermodynamics gets affected by such long wavelength fluctuations? To answer it, one is required to test the AG relation without (standard) and with (cage-relative) correction due to the long-wavelength fluctuations. The AG relation predicts a linear dependence of $\ln\tau_{\alpha}$ on $(TS_{c})^{-1}$ at all spatial dimension $D$, Eqn.1. According to the RFOT theory Kirkpatrick2015 ; Lubchenko2007 ; Bouchaud2004 ; Biroli2012 ; Starr2013 ; Karmakar2015 , the growing many body correlation at lowering temperature is characterized by a growing static correlation length $\xi$. The free energy barrier for structural relaxation $\Delta G$ scales with $\xi$ as $\Delta G\sim\xi^{\psi}$. Thus the characteristic time scale for structural relaxation scales with $\xi$ as: $\tau_{\alpha}\sim\exp({\Delta G\over k_{B}T})\sim\exp({\xi^{\psi}\over k_{B}T})$. Further, in the RFOT scenario, the correlation length is related to the configuration entropy as: $\xi^{D-\theta}\sim{1\over S_{c}}$. Hence the relation between time and entropy becomes explicitly dimension dependent:

The RFOT exponents $\psi$ and $\theta$ are in general unknown. However, they should satisfy the inequality Starr2013 : $\theta\leq\psi\leq D-1$. Note that testing the RFOT relations require computation of the lengthscale $\xi$ which is beyond the scope of the present study. However, guided by these scaling relations we fit the data to a non-linear, generalized AG form:

The AG relation is recovered Starr2013 ; Kirkpatrick1989 for $\psi=D/2=\theta$.

Fig. 4 tests the relationship between the $\alpha$-relaxation timescale and the configurational entropy for the 2DMKA and the 2DR10 models. In all cases the temperature range is chosen upto the onset temperature. Here the configurational entropy is computed using the harmonic approximation, to isolate the effect of long wavelength fluctuation on the AG relation. The inset panels show the timescales $\tau_{\alpha}$ without correcting for the long wavelength fluctuation, obtained from both the overlap function (Ref Sengupta2012 ) as well as the $F_{s}(K,t)$. Comparing the timescales from the intermediate scattering function without (standard) and with (cage-relative) correction due to the long wavelength fluctuation, we see its main effect is quantitative: in both cases the AG relation still breaks down, and instead follows the non-linear form in Eqn. 18. However, removing the effects of long wavelength fluctuation changes the extent of the deviation from linearity, i.e. values of the exponent $\alpha\equiv{\psi\over D-\theta}$, see Table 3. Interestingly, for the 2DMKA model, $\alpha$ increases towards 1 wheres for the 2DR10 model, $\alpha$ increases away from 1. Thus the nature of deviation depends on the nature of the interatomic potential.

Details of configurational entropy calculation in harmonic approximation (Eqn. 7) for the 2DMKA and the 2DR10 models have been reported elsewhere Sengupta2012 . Here we describe the effect of adding anharmonic correction to $S_{c}$ (Eqn. 14). First, in Fig. 5 we show the temperature dependence of the anharmonic contribution to the potential energy characterized by polynomial fit functions, see Sec II and Eqns. 9, 10 and 12.

Analysis within the PEL framework has shown that for harmonic approximation to the vibrational entropy, $TS_{c}$ should be linear in $T$: $TS_{c}=[K_{T}(\frac{T}{T_{K}}-1)]$, where $S_{c}(T_{K})=0$ determines the Kauzmann temperature $T_{K}$ and $K_{T}$ is the thermodynamic fragility Sastry2001 ; Novikov2022 . This expectation is well tested for 3D glass-former models in experiments Richert1998 , and also in simulation studies without Banerjee2016 and with anharmonic corrections Nandi2022 although deviations have also been noted Ozawa2019 . Recently Berthier et al. has raised the intriguing possibility that the Kauzmann temperature in 2D is, in fact, zero Berthier2019 . In other words, there is no ideal glass transition at non-zero temperature in 2D. Given these background, the effect of the correction to $S_{c}$ due to vibration entropy - see Eqns. 11, 13 - on its $T$ dependence becomes especially interesting in 2D to which we turn our attention next.

In Fig. 6 we show the $S_{c}$ data over a range of low temperatures upto the onset temperature for 2DMKA (top row) and 2DR10 (bottom row) models. First note that when the configurational entropy is computed in harmonic approximation, left columns of Fig. 6, $TS_{c}^{har}(T)$ is indeed linear with $T$. However the Kauzmann temperature $T_{K}$ is siginificantly different from the VFT divergence temperature $T_{0}$ in both the models, see Tables 1, 2. We also fit the data using a non-linear function $TS_{c}=[K_{T}(\frac{T}{T_{K}}-1)]^{\frac{1}{\alpha_{1}}}$ with $\alpha_{1}$ as a free fit parameter. The choice of this fit function is motivated by Eqn. 18, although other fit functions have also been used in the literature Ozawa2019 . This analysis results in the exponent $\alpha_{1}\sim 1$ within numerical accuracy providing a robust test of the linearity of $TS_{c}^{har}(T)$.

To the contrary, we find that in 2D glass-formers, adding anharmonic correction makes the $TS_{c}^{anhar}$ distinctly non-linear with $T$. This behaviour is the same in both models and using both methods 1 and 2 for computing the anharmonic vibrational entropy (see Eqns. 10-13). However, estimates of $T_{K}$ depends on the models and the methods of computing the anharmonic correction to $S_{c}$. On one hand, in the 2DMKA model, the non-linear fit estimate of $T_{K}$ values are significantly lower when anharmonic contribution to vibrational entropy is considered, and approximately same for both methods 1 and 2. On the other hand, in the 2DR10 model, method 1 results in $T_{K}\sim T_{0}$ while method 2 produces a $T_{K}$ closer to zero, see Table 2. In order to gain further insight, on the right columns of Fig. 6, we show the temperature dependence of $S_{c}$. We analyze the $T$ dependence using several plausible fitting functions which can be organized into two families: one assuming $T_{K}=0$ Berthier2019 and the other assuming $T_{K}\neq 0$. However, both families produce fits of comparable quality. Hence in the range of temperature in which configurational entropy is analyzed in the present study, our data is consistent with both finite and zero temperature ideal glass transition scenarios.

Finally, in Fig. 7, we test the AG relation with (main panel) and without (inset) including the anharmonic correction to the configurational entropy. Since the timescales at low temperatures gets decoupled (see Fig. 3), we consider both $\alpha$-relaxation time $\tau_{\alpha CR}$ as well as diffusion coefficient $D_{CR}$. Note that the correction due to the effect of long wavelength fluctuations are already incorporated. Thus Fig. 7 constitutes the main result of the present work. The top row of Fig. 7 shows data for 2DMKA model and the bottom row depicts 2DR10 model data. Lines are fits to Eqn. 18 from which the scaling exponent $\alpha$ is extracted and tabulated in Table 3. First, for the 2DMKA model, comparing the inset and the main panel, we see that the main effect of anharmonic correction is to reduce the deviation from the AG relation, Eqn. 1, i.e. to increase the value of $\alpha$ towards 1. This behaviour remains same for both $\tau_{\alpha CR}$ and $D_{CR}$. Indeed, for the 2DMKA model, the AG relation is recovered within numerical accuracy when both the correction to $\tau_{\alpha CR}$ due to long wavelength fluctuation and the correction to $S_{c}$ due to anharmonic vibration are considered. Even with $D_{CR}$, the AG relation is only weakly non-linear for this model. To the contrary, the 2DR10 model shows a qualitatively different behaviour. Even after accounting for all the corrections in timescales and entropy, one still observes a breakdown of the AG relation in two dimensions. Thus we clearly reveal an intrinsic dimension dependence that indicates that the behavior of 2D glass-formers are different from those in 3D. Interestingly, the nature of inter-particle interaction has a more fundamental effect in determining the value of the exponent $\alpha$ than long wavelength fluctuations and anharmonic correction. Although, even in 2DR10 model, including the anharmonic correction reduces the extent of deviation from linearity, i.e. decreases $\alpha$ towards 1.