Does the Adam-Gibbs relation hold in simulated supercooled liquids?

TL;DR

This study tests the validity of the Adam-Gibbs relation in simulated supercooled liquids and finds it generally violated, suggesting thermodynamics alone may not explain slow glassy dynamics.

Contribution

The paper provides the first comprehensive test of the Adam-Gibbs relation in simulations and experiments, challenging its universality in glass transition theories.

Findings

Adam-Gibbs relation is violated in simulated models at relevant timescales.

Experimental data show similar deviations, questioning thermodynamics as the sole driver.

Deviations are consistent with random first order transition theory or alternative mechanisms.

Abstract

We perform stringent tests of thermodynamic theories of the glass transition over the experimentally relevant temperature regime for several simulated glass-formers. The swap Monte Carlo algorithm is used to estimate the configurational entropy and static point-to-set lengthscale, and careful extrapolations are used for the relaxation times. We first quantify the relation between configurational entropy and the point-to-set lengthscale in two and three dimensions. We then show that the Adam-Gibbs relation is generally violated in simulated models for the experimentally relevant time window. Collecting experimental data for several supercooled molecular liquids, we show that the same trends are observed experimentally. Deviations from the Adam-Gibbs relation remain compatible with random first order transition theory, and may account for the reported discrepancies between Kauzmann and Vogel-Fulcher-Tammann temperatures. Alternatively, they may also indicate that even near $T_g$ thermodynamics is not the only driving force for slow dynamics.

Tables

Table 1. Table 1: Fitting parameters for the configurational entropy ( A 𝐴 A , B 𝐵 B , and T K subscript 𝑇 K T_{\rm K} ), for the relaxation time ( τ o subscript 𝜏 𝑜 \tau_{o} , T o subscript 𝑇 𝑜 T_{o} , and C 𝐶 C ) and kinetic fragility index m 𝑚 m for the simulated models. Note that Monte Carlo dynamics (HS3D, SSV3D, SSV2D) and molecular dynamics (SSP3D) have different time units.

Model	$A$	$B$	$T_K$	${\log }_{10}{\tau }_o$	$T_o$	$C$	$m$
HS3D	3.208	-37.33	0.0251	3.88	0.063	22.72	45.5
SSV3D	1.495	-1.92	0.0386	3.02	0.266	3.15	32.0
SSP3D	2.082	-1.74	0.2902	0.41	0.961	16.77	42.4
SSV2D	0.453	1.89	-	2.40	1.006	0.25	31.2

Equations

$$\log ({\tau }_{\alpha }/{\tau }_0)\propto \frac{1}{T{S}_{\mathrm{conf}}},$$

$${\xi }_{\mathrm{pts}}\propto S_{\mathrm{conf}}^{-1/(d-\theta )},$$

$$\log ({\tau }_{\alpha }/{\tau }_0)\propto {\xi }_{\mathrm{pts}}^{\psi }/T,$$

$$\log ({\tau }_{\alpha }/{\tau }_0)\propto \frac{1}{T{S}_{\mathrm{conf}}^{\alpha }},$$

$$\alpha =\frac{\psi }{d-\theta }.$$

$$v_{ij}(r)$$

$${\sigma }_{ij}$$

$$\log ({\tau }_{\alpha }/{\tau }_0)\propto (T-T_{\mathrm{VFT}}{)}^{-1}.$$

$${\tau }_{\alpha }^{\mathrm{para}}={\tau }_o\exp [C(T_o/T-1{)}^2]$$

$${\tau }_{\alpha }={\tau }_0\exp \left[\frac{K}{T}\exp [C/T]\right],$$

$${\tau }_{\alpha }={\tau }_0\exp [A/T^n],$$

$$Y(T,P)={\int }_0^{\infty }\mathrm{d}V{e}^{-\beta \left(PV+F(T,V)\right)}.$$

$$\rho (V|T,P)=\frac{e^{-\beta \left(PV+F(T,V)\right)}}{Y(T,P)}.$$

$$\rho (V|T,P)=\frac{1}{\sqrt{2\pi {\sigma }_V^2}}\exp \left[-\frac{(V-V_{\ast }{)}^2}{2{\sigma }_V^2}\right],$$

$$S_{\mathrm{tot}}(T,P)$$

$$S(T,P)=⟨S(T,V){⟩}_{T,P}-⟨\ln \rho (V|T,P){⟩}_{T,P}.$$

$$-\frac{1}{N}⟨\ln \rho (V|T,P){⟩}_{T,P}=\frac{1}{N}\ln \sqrt{2\pi e{\sigma }_V^2}.$$

Full text

I Introduction

Since its first derivation in 1965 Adam and Gibbs (1965), the Adam-Gibbs relation has played a central role in glass transition studies Berthier and Biroli (2011), since it is at the core of thermodynamic approaches to the glass problem Adam and Gibbs (1965); Berthier and Biroli (2011); Kirkpatrick et al. (1989); Bouchaud and Biroli (2004); Lubchenko and Wolynes (2007); Biroli and Bouchaud (2012); Lubchenko (2015); Charbonneau et al. (2017); Dudowicz et al. (2008). The Adam-Gibbs relation captures in a simple mathematical form the physical idea that the decrease of the configurational entropy $S_{\rm conf}$ controls the growth of the relaxation time $\tau_{\alpha}$ as the experimental glass transition temperature $T_{g}$ is approached:

[TABLE]

where $\tau_{0}$ is a microscopic timescale. Testing the Adam-Gibbs relation has almost become synonymous to testing the thermodynamic nature of glass formation Angell (1997); Richert and Angell (1998); Sastry (2001); Dyre et al. (2009).

Since computational methods have become available in the early 2000’s to measure the configurational entropy in numerical simulations Sciortino et al. (1999); Sastry (2000); Berthier et al. (2019a), the Adam-Gibbs relation has been tested in a large number of studies using many different models of glass-forming materials Sastry (2001); Mossa et al. (2002); Sciortino (2005); Saika-Voivod et al. (2001); Angelani et al. (2005); Sengupta et al. (2012); Starr et al. (2013); Parmar et al. (2017); Handle and Sciortino (2018). Importantly, these simulations are all restricted to a high temperature regime (typically above the mode-coupling crossover temperature $T_{\rm mct}$ Götze (2008)) that barely overlaps with the corresponding experimental studies. In addition simulations typically cover a dynamic window of at most 3-4 decades, much narrower than in experimental studies. Despite these caveats, the general consensus is that the Adam-Gibbs relation is generally valid in the regime accessed by the simulations. In experiments, which typically analyse temperatures close to $T_{g}$, the Adam-Gibbs relation seems again to be well obeyed for a range of materials Magill (1967); Takahara et al. (1995); Angell (1997); Richert and Angell (1998); Ngai (1999); Alba-Simionesco (2001); Roland et al. (2004); Cangialosi et al. (2005); Masiewicz et al. (2015). Yet, experiments indicate as well that the Adam-Gibbs relation does not hold anymore above a temperature scale close to $T_{\rm mct}$ Richert and Angell (1998); Ngai (1999), in stark contrast with the numerical results. Systematic deviations from the Adam-Gibbs relation were also reported below $T_{\rm mct}$ for some systems Ngai (1999); Roland et al. (2004), but imprecise entropy measurements or inappropriate timescale determinations have been invoked to rationalise them.

In the last three decades, the random first order transition (RFOT) theory of the glass transition Kirkpatrick et al. (1989); Lubchenko and Wolynes (2007) has revisited the Adam-Gibbs relation in greater depth Lubchenko and Wolynes (2007); Bouchaud and Biroli (2004); Biroli and Bouchaud (2012); Lubchenko (2015) to provide an increasingly precise description of the connection between thermodynamics and dynamics in supercooled liquids. This connection can be decomposed in two steps. First, the decrease of the configurational entropy is shown, by a purely thermodynamic reasoning Bouchaud and Biroli (2004), to give rise to a growing ‘point-to-set’ static correlation lengthscale:

[TABLE]

where an interface exponent $\theta$ is introduced. In the simplest approximation, one has $\theta=d-1$ which corresponds to a (hyper-)surface in a space of dimension $d$. The value $\theta=d/2$ was also proposed Kirkpatrick et al. (1989); Lubchenko (2015), to take into account finite dimensional surface fluctuations due to the disordered nature of the amorphous phase. More generally, the inequality $\theta\leq d-1$ is expected to hold. Second, the connection to dynamics is made via the assumption that relaxation in the liquid for $T<T_{\rm mct}$ proceeds via thermally activated events correlated over a lengthscale $\xi_{\rm pts}$, resulting in the general relation Kirkpatrick et al. (1989); Bouchaud and Biroli (2004),

[TABLE]

where $\psi$ is a dynamical exponent. Various theoretical and numerical estimates of $\psi$ have been proposed Bouchaud and Biroli (2004); Cammarota et al. (2009a); Karmakar et al. (2009); Hocky et al. (2012); Gutiérrez et al. (2015). In the original paper by Kirkpatrick et al. Kirkpatrick et al. (1989), $\psi=\theta=d/2$ was assumed and so only one exponent had been introduced.

Using Eqs. (2, 3), one finds a generalised version of the Adam-Gibbs relation,

[TABLE]

with a non-trivial exponent

[TABLE]

This shows that $\alpha$ may or may not be equal to unity, depending on the relative values of the two independent exponents $\psi$ and $\theta$. As a consequence, Eq. (4) may or may not be equivalent to Eq. (1).

To our knowledge, a direct test of Eqs. (3, 4, 5) in the theoretically-motivated temperature regime, employing appropriate observables, has never been performed. Most previous simulations have considered a temperature regime $T\gtrsim T_{\rm mct}$ Sastry (2001); Mossa et al. (2002); Cavagna et al. (2012); Sengupta et al. (2012) where the physics is expected to be non-activated and the configurational entropy and point-to-set lengthscales are not well-defined. This is of course valuable work, but theory itself suggests that the tested scaling relations have no reason to hold in this temperature regime. Experiments instead access the correct temperature regime, but cannot easily measure the point-to-set correlation lengthscale. As a proxy, Refs. Capaccioli et al. (2008); Brun et al. (2012) replaced $\xi_{\rm pts}$ by the lengthscale of dynamic heterogeneities that can be more easily estimated experimentally Berthier et al. (2005). Many other experimental studies study Eq. (1) directly near $T_{g}$ Richert and Angell (1998); Roland et al. (2004).

In this work, we take advantage of the progress allowed by the swap Monte Carlo algorithm Berthier et al. (2016a); Ninarello et al. (2017) to measure directly in several numerical models the temperature dependence of the configurational entropy and point-to-set lengthscale down to $T_{g}$. For the dynamics, we build on previous work Ninarello et al. (2017) and provide additional experimental support showing that one can safely estimate the temperature dependence of the relaxation time also down to $T_{g}$, using a careful fitting procedure. We collect data from earlier works Berthier et al. (2019b, 2017); Ozawa et al. (2018) that we extend where needed, and perform new simulations for one additional model.

As a result, we are in a position to provide for the first time stringent tests of the Adam-Gibbs relation and of RFOT theory for computer models simulated in the same regime as in experiments. Our results suggest that the Adam-Gibbs relation is generally not valid in computer models in the experimental regime $T_{g}<T<T_{\rm mct}$. To test our findings against experiments, we collect high-quality thermodynamic and dynamic data for several supercooled liquids (most of which are obtained by state-of-the-art thermodynamic measurements Tatsumi et al. (2012)), and reach similar conclusions. Overall, we find that Eq. (1) is not obeyed for most systems, while Eq. (4) is obeyed with an exponent $\alpha$ that fluctuates weakly from system to system, with typically $\alpha<1$. Our findings can be taken either as a confirmation that RFOT theory works well, with a non-trivial set of critical exponents, or that a small $\alpha<1$ exponent indicates that thermodynamics is not the only driving force for the dynamic slowdown near $T_{g}$.

This paper is organised as follows. In Sec. II we present the numerical methods used to obtain the configurational entropy, the point-to-set lengthscale, and the relaxation time. We also describe our choice of experimental data to reliably test the Adam-Gibbs relation over a broad range of temperatures. In Sec. III we present the results of our analysis of the exponents $\theta$ and $\alpha$ in simulations, then in experiments. We discuss the physical meaning of our results in Sec. IV.

II Description of the data

In order to analyse quantitatively the connection between dynamic and thermodynamic properties, we collect and extend data from previous numerical works. We also collect data from selected published experimental works, and motivate our selection.

II.1 Numerical models

The recent development of the swap Monte Carlo algorithm allows us to access very low-temperature equilibrium configurations in computer simulations. In particular, the temperature regime $T_{g}<T<T_{\rm mct}$ can be comfortably accessed. This temperature regime is the correct one to test thermodynamic theories, as it is precisely where they should apply, and it corresponds to the regime explored experimentally.

We gather simulation data for polydisperse systems using a continuous size distribution Ninarello et al. (2017). The particle diameters $\sigma$ are distributed between $\sigma_{\rm min}$ and $\sigma_{\rm max}$ from $f(\sigma)=c\sigma^{-3}$, where $c$ is a normalization constant and $\sigma_{\rm min}/\sigma_{\rm max}=0.45$. We use the average diameter $\overline{\sigma}$ as the unit length.

We study four numerical models: three-dimensional additive hard spheres (HS3D) Berthier et al. (2016a), two and three dimensional non-additive soft disks (SSV2D) Berthier et al. (2019b) and spheres (SSV3D) Ninarello et al. (2017) under an isochoric path. We also perform new simulations of three-dimensional non-additive soft spheres (SSP3D), under an isobaric path. To thermalize the last model, we use an hybrid molecular dynamics/swap Monte Carlo scheme Berthier et al. (2018).

We use the following pairwise potential for the polydisperse soft sphere/disk models Ninarello et al. (2017),

[TABLE]

where $v_{0}$ is the energy unit, and $\epsilon$ quantifies the degree of non-additivity of the system. We set $\epsilon=0.2$ for SSV3D and SSV2D, and $\epsilon=0.1$ for SSP3D. The constants, $c_{0}$, $c_{1}$ and $c_{2}$, are chosen to smooth $v_{ij}(r)$ up to its second derivative at the cut-off distance $r_{\rm cut}=1.25\sigma_{ij}$. We set the number density $\rho=N/L^{3}=1.02$ with $N=1500$ for SSV3D, and $\rho=N/L^{2}=1.01$ with $N=1000$ for SSV2D. For SSP3D, the pressure on the isobaric path is $P=30.0$. For HS3D Berthier et al. (2016a), the pair interaction is zero for non-overlapping particles and infinite otherwise. The relevant control parameter for hard spheres is the reduced pressure $p=P/(\rho T)$. For hard spheres, $1/p$ plays precisely the same role as temperature $T$ for a dense liquid Berthier and Witten (2009), and there is no distinction between isochoric and isobaric paths.

Relaxation times for HS3D, SSV3D, and SSV2D are measured in units of MC sweeps, which comprise $N$ Monte Carlo trial moves. For SSP3D, the relaxation time is expressed in units of $\sqrt{v_{0}/m\overline{\sigma}^{2}}$, where $m$ is the mass of the particles.

II.2 Configurational entropy and point-to-set length

The configurational entropy $S_{\rm conf}$ is measured from configurations generated with swap Monte Carlo simulations. It is defined as $S_{\rm conf}=S_{\rm tot}-S_{\rm glass}$, where $S_{\rm tot}$ and $S_{\rm glass}$ are the total and glass entropies, respectively Berthier et al. (2019a). $S_{\rm tot}$ and $S_{\rm glass}$ are computed using thermodynamic integration schemes, as explained in Ref. Ozawa et al. (2018). In Appendix A we describe how to measure $S_{\rm conf}$ along an isobaric path using constant pressure simulations for SSP3D, as this was not documented before.

Figure 1 shows the configurational entropy that we use for latter analysis. The data for $S_{\rm conf}(T)$ are normalized by the values at the mode coupling crossover $T_{\rm mct}$, whose value is determined by a power law fit to the dynamic relaxation time data Götze (2008). The actual values are $T_{\rm mct}=0.0426$, $0.104$, $0.556$, and $0.123$ for HS3D, SSV3D, SSP3D, and SSV2D, respectively.

In order to increase the accuracy of the analysis, we employ empirical fitting functions. For the three-dimensional models, we use a conventional fitting function plus a quadratic correction, $TS_{\rm conf}=A(T-T_{\rm K})+B(T-T_{\rm K})^{2}$ Richert and Angell (1998); Banerjee et al. (2014). For the two-dimensional model, we use $1/S_{\rm conf}=A/T+B$ Berthier et al. (2019b). These fitting functions conveniently enable us to incorporate $S_{\rm conf}$ values in between actual data points and help us determining the exponents. The fitting parameters are presented in Table 1.

We also collect the point-to-set lengthscale $\xi_{\rm pts}$ data for SSV2D Berthier et al. (2019b) and HS3D Berthier et al. (2017), obtained from recently developed computational methods Biroli et al. (2008); Berthier et al. (2016b). Together with $S_{\rm conf}$, the data for $\xi_{\rm pts}$ will allow us to estimate the exponent $\theta$ using Eq. (2).

II.3 Relaxation times

Dynamical information is obtained using either standard Monte Carlo (for HS3D, SSV3D, SSV2D) or molecular dynamics (for SSP3D). The equivalence between the two types of dynamics is well documented Berthier and Kob (2007). Both Monte Carlo and molecular dynamics simulations are run starting from initial configurations that are obtained using the swap Monte Carlo algorithm. This procedure allows us to cover about 5 orders of magnitude of relevant slow dynamics.

The relaxation time $\tau_{\alpha}$ is measured by the self-intermediate scattering function in three dimensional models. For the two-dimensional model, we use the autocorrelation function of the bond-orientational order parameter, which is insensitive to the long-range Mermin-Wagner fluctuations that are specific to $d=2$ Flenner and Szamel (2015).

The relaxation time $\tau_{\alpha}$ for HS3D Berthier et al. (2017), SSV3D Ninarello et al. (2017), SSP3D (new to this work), and SSV2D Berthier et al. (2019b) is shown in Figure 2. The data are normalized using an onset temperature $T_{o}$ for the emergence of slow dynamics, determined from the fitting procedure described below, and define $\tau_{o}=\tau_{\alpha}(T=T_{o})$. Clearly, all simulation data show a non-Arrhenius temperature dependence of the relaxation time, which demonstrates that our models describe fragile glass-formers.

The swap numerical schemes allow us to prepare equilibrated configurations at very low temperatures. Because they involve non-physical particle dynamics, one cannot use them to measure the relaxation time of the physical dynamics in this low-temperature regime. Therefore, we need to extrapolate the relaxation time from the regime where $\tau_{\alpha}$ can be measured to the experimental regime, where this is unachievable.

We start by employing the Vogel-Fulcher-Tammann (VFT) law:

[TABLE]

where $\tau_{0}$ and $T_{\rm VFT}$ are fitting parameters. We fitted this function on our numerical data over the accessible time window and we concluded that it performs very badly when extrapolated at lower temperatures. We found for instance that the swap Monte Carlo algorithm easily thermalises at temperatures below the extrapolated VFT critical temperature $T_{\rm VFT}$, which invalidates directly its use to describe numerical data Ninarello et al. (2017). The inability of the VFT law to describe experimental data over a wide range of temperature was discussed in detail in Refs. Stickel et al. (1995); Blodgett et al. (2015).

It has been found in previous experimental studies that the parabolic law

[TABLE]

fits accurately the data over a very large temperature range Elmatad et al. (2010); Mauro et al. (2009). Its fitting parameters are $\tau_{o}$, $C$, and $T_{o}$.

In addition to the VFT and parabolic laws, we consider two other functional forms, shown in Fig. 2. One is a double exponential equation (MEYGEA) discussed in Refs. Mauro et al. (2009); Elmatad et al. (2010):

[TABLE]

where $\tau_{0}$, $K$, and $C$ are the fitting parameters. The other one is the Avramov and Milchev (AM) equation Avramov and Milchev (1988) given by

[TABLE]

where $\tau_{0}$, $A$, and $n$ (real exponent) are the fitting parameters. All the fitting functions considered in this paper have three free-fitting parameters which is the minimal number to mathematically characterize non-Arrhenius behavior. Given the small variation of the apparent activation energy over the dynamic range studied experimentally, it is not surprising that several smooth functions of temperature can describe the evolution of $\log(\tau_{\alpha})$. Figure 2 shows that different fitting functions produce slight variations in the extrapolated value for $T_{g}$. The key issue is therefore to choose the best fitting function, i.e., the one from which the low temperature data can be inferred accurately from the high temperature one.

To find the best fitting procedure, we train on experimental data (see Appendix B). We fit the above four equations to the data, restricting ourselves to a modest dynamic range, comparable to numerical timescales. We then extrapolate to temperatures close to $T_{g}$, and compare the extrapolation to the actual data. We find excellent agreement when using the parabolic law for the experimental data with kinetic fragility indexes similar to our numerical models, which validates further our procedure. Thus, we empirically find that fitting the parabolic law to the numerical time window provides an excellent description of the data close to $T_{g}$, as reported previously Elmatad et al. (2010). This is a purely practical choice, and we make no assumption about the physical mechanism which could lead to such a law.

By using the fitting parameter $\tau_{o}$ obtained from the parabolic law, we define two time windows. First we define the simulation window by $\tau_{\alpha}/\tau_{o}\in[10^{0},10^{5}]$. The upper bound of this timescale corresponds to recent simulation studies with very long timescales Berthier et al. (2017); Coslovich et al. (2018). The experimental window is defined by $\tau_{\alpha}/\tau_{o}\in[10^{3},10^{12}]$. The lower bound corresponds to a timescale around the mode-coupling crossover $T_{\rm mct}$ ($\tau_{\alpha}\simeq 10^{-7}$ s Novikov and Sokolov (2003)), and the upper bound corresponds to the timescale at the experimental glass transition $T_{\rm g}$ ($\tau_{\alpha}\simeq 100$ s). The experimental window is therefore the appropriate regime to test the predictions made by the RFOT theory. Notice that in this paper, we neither try to go below $T_{g}$, nor to examine the fate of supercooled liquids at even lower temperature Royall et al. (2018).

For numerical models, we determine the experimental glass transition temperature $T_{\rm g}$ as $\tau_{\alpha}^{\rm para}(T_{\rm g})/\tau_{o}=10^{12}$. The kinetic fragility index $m$ is determined by $m=\partial\log_{10}\tau_{\alpha}^{\rm para}/\partial(T_{\rm g}/T)|_{T=T_{\rm g}}$. The fitting parameters and fragility indexes are given in Table 1.

II.4 Experimental data

We select materials for which high-quality data for the configurational entropy and relaxation time over a broad temperature range is available in the literature. This allows for a comparison with computer simulations and an accurate determination of the exponent $\alpha$ in Eq. (4).

We select 2-methyl tetrahydrofuran (2MTHF), ethylbenzene (ETB), ethanol, glycerol, o-terphenyl (OTP), 1-propanol, propylene carbonate (PC), salol, toluene and 3-bromopentane. The configurational entropy data for 2MTHF, ETB, OTP, PC, salol, and toluene were recently obtained from accurate experiments by Tatsumi, Aso and Yamamuro. Some of the data is presented in Ref. Tatsumi et al. (2012). The data for 1-propanol is taken from Ref. Takahara et al. (1994). In these data for all the above materials, $S_{\rm conf}$ is measured by thermodynamic integration of the heat capacity difference between supercooled liquids and non-equilibrium glasses. This treatment should be conceptually better than using the crystal entropy Richert and Angell (1998), but this is still a rather crude approximation Yoshimori and Odagaki (2011), whose accuracy is expected to be material-dependent Smith et al. (2017). For ethanol Takeda et al. (1999); Haida et al. (1977), glycerol Takeda et al. (1999); Gibson and Giauque (1923), and 3-bromopentane Takahara et al. (1995), $S_{\rm conf}$ is obtained using the crystal entropy $S_{\rm cry}$, i.e., $S_{\rm conf}=S_{\rm liq}-S_{\rm cry}$. Notice that we do not seek to present thermodynamic data for ultrastable glasses prepared below $T_{g}$, even though we believe that these materials can be instrumental to test more precisely glass transition theories Beasley et al. .

The relaxation time data are mainly obtained from dielectric measurements, but some data are combined with other methods, such as viscosity measurements. The corresponding references are: 2MTHF Richert and Angell (1998), ETB Chen and Richert (2011); Barlow et al. (1966); Rossini (1953), ethanol Tatsumi et al. (2012), glycerol Schneider et al. (1998); Lunkenheimer et al. (2000, 2005), OTP Schmidtke et al. (2012), 1-propanol Hansen et al. (1997); Richert and Angell (1998); Sillrén et al. (2014), PC Schneider et al. (1999); Lunkenheimer et al. (2000, 2005), salol Stickel et al. (1996), toluene Schmidtke et al. (2012), and 3-bromopentane Berberian and Cole (1986).

For the experimental data, we set $\tau_{o}=10^{-10}$ s. Therefore the simulation and experimental time windows correspond to $\tau_{\alpha}\in[10^{-10}{~{}\rm s},10^{-5}{~{}\rm s}]$ and $\tau_{\alpha}\in[10^{-7}{~{}\rm s},10^{2}{~{}\rm s}]$, respectively. In particular, $T_{g}$ corresponds to the standard relaxation time $\tau_{\alpha}=100{~{}\rm s}$.

The configurational entropy and relaxation time data for the materials presented above are gathered in Fig. 3, together with empirical quadratic fits to the configurational entropy.

III Results

In this section, we perform a test of Eqs. (1, 2, 3, 4, 5) using the experimental and numerical data presented in Sec. II. We first study Eq. (2) using numerical data for $\xi_{\rm pts}$ and $S_{\rm conf}$ to estimate $\theta$. Then, we estimate $\alpha$ in Eq. (4) by comparing $\tau_{\alpha}$ and $S_{\rm conf}$ using both computer simulations and experiments to investigate the validity of the Adam-Gibbs relation in Eq. (1). Finally, the values of $\theta$ and $\alpha$ allow us to discuss that taken by $\psi=(d-\theta)\alpha$, deduced from Eq. (5).

III.1 The static exponent $\theta$

First we estimate the exponent $\theta$ in Eq. (2) combining independent data obtained for $S_{\rm conf}$ and $\xi_{\rm pts}$.

Figure 4 shows a log-log plot of $S_{\rm conf}$ versus $\xi_{\rm pts}$ for three dimensional polydisperse hard spheres (HS3D) (a) and two dimensional soft disks (SSV2D) (b). We use the fitted functional form for $S_{\rm conf}$ (obtained in Fig. 1), whereas the actual data points are used for $\xi_{\rm pts}$. We emphasize that while temperature is a running parameter in this plot, the data point in Fig. 4 correspond to the regime of interest $T<T_{\rm mct}$. Such results have never been achieved, as earlier numerical work were all performed for $T>T_{\rm mct}$, or only slightly below $T_{\rm mct}$ Cammarota et al. (2009a). Despite the larger temperature range explored in this work, we are fully aware the relative variation of $\xi_{\rm pts}$ and $S_{\rm conf}$ remain fairly modest, which makes the determination of a critical exponent quite difficult.

For HS3D, we report two estimates for $S_{\rm conf}$, obtained from different schemes. One is a generalized Frenkel-Ladd (GFL) method Frenkel and Smit (2002); Ozawa et al. (2018), and the other is the Franz-Parisi (FP) free energy method proposed earlier Franz and Parisi (1997); Berthier and Coslovich (2014); Berthier et al. (2017). The exponent $\theta$ is extracted by fits to straight lines, whose slope gives $\theta-d$, see Eq. (2). We obtain $\theta=1.35\pm 0.06$ for GFL and $\theta\simeq 1.84\pm 0.09$ for FP. These values are compatible with either the theoretical prediction $\theta=d/2$ by Kirkpatrick et al. Kirkpatrick et al. (1989), or with that of Franz $\theta=d-1$ Franz (2005).

We obtain $\theta=1.12\pm 0.02$ for SSV2D. This value is close to both theoretical predictions, $\theta=d/2$ and $\theta=d-1$, which coincide in $d=2$, giving $\theta=1$ . Obviously, one cannot discriminate between the two predictions.

Overall, we find that for $d=3$ the value measured for $\theta$ conforms with the two available predictions, which is an encouraging result from the viewpoint of RFOT theory. Unfortunately, the obtained values fall in-between the two predictions, which are too close to be discriminated. We suggest that performing point-to-set and configurational entropy measurements in $d=4$, combining recently developed tools Berthier et al. (2016b, 2019c); Ozawa et al. (2018), would be very useful to conclude on this point. Indeed, when $d=4$, the two predictions yield $\theta=d/2=2$ and $\theta=d-1=3$, which are further appart than in $d=3$.

III.2 Breakdown of the Adam-Gibbs relation and numerical estimation of $\alpha$

We next examine the validity of Eq. (4) by connecting $\tau_{\alpha}$ and $S_{\rm conf}$, and estimating the exponent $\alpha$. When $\alpha=1$, the Adam-Gibbs relation in Eq. (1) is recovered.

In Fig. 5(a,c,e,g) we show conventional Adam-Gibbs plots where the evolution of $\log_{10}(\tau_{\alpha}/\tau_{o})$ is represented as a function of $1/(Ts_{\rm conf})$, where $s_{\rm conf}=S_{\rm conf}/N$, for hard spheres (HS3D) (a), soft spheres along the isochoric path (SSV3D) (c), along the isobaric path (SSP3D) (e), and the soft disks (SSV2D) (g). We combine the dynamic and thermodynamic data described in Sec. II, restricted to the experimental time window ($\tau_{\alpha}/\tau_{o}\in[10^{3}-10^{12}]$). We use the fitted functional forms for both $\tau_{\alpha}$ (estimated in Fig. 2) and $S_{\rm conf}$ (obtained in Fig. 1), which produces “continuous curves” instead of a discrete data points. To our knowledge, this is the first time that the Adam-Gibbs relation is tested for computer models over the time window where it is actually supposed to apply.

For all three-dimensional models, we find that $\log_{10}(\tau_{\alpha}/\tau_{o})$ is a concave function of $1/Ts_{\rm conf}$, whereas it is convex for the two-dimensional model. If tested over a narrow time window close to $T_{\rm mct}$, an acceptable linear behaviour could possibly be observed, that would suggest the validity of the Adam-Gibbs relation, in agreement with many earlier findings Sastry (2001); Mossa et al. (2002); Sciortino (2005); Saika-Voivod et al. (2001); Angelani et al. (2005); Sengupta et al. (2012); Starr et al. (2013); Parmar et al. (2017); Handle and Sciortino (2018). The trend that we report here appears to contrast with recent results obtained in the Kob-Andersen model, where slight convexity and concavity are respectively observed in $d=3$ Parmar et al. (2017) and $d=2$ Sengupta et al. (2012). These results were however obtained in the numerical time window, above $T_{\rm mct}$. Our results demonstrate that when observed over a much broader range, and closer to $T_{g}$, the Adam-Gibbs relation is actually not obeyed for any of the numerical models studied here.

The clear violations of the standard Adam-Gibbs relation that we find over the experimental time window imply that the exponent $\alpha$ must deviate from the value $\alpha=1$. We varied its value around unity and used it as a free parameter to obtain generalised Adam-Gibbs plots, which are shown in Fig. 5(b,d,f,h) for the same numerical models. All plots now show a perfect straight line, suggesting that the introduction of the parameter $\alpha$ is sufficient to describe the data. We obtain $\alpha=0.24$, $0.49$, $0.72$, and $1.89$, for HS3D, SSP3D, SSV3D, and SSV2D, respectively, so that $\alpha<1$ for the three dimensional models, whereas $\alpha>1$ for the two dimensional model.

Since the four models we have simulated all display violations of the Adam-Gibbs relation, we conclude that Eq. (1) does not describe well the physics of simulated supercooled liquids when analysed over the experimental time window. Additional models should be studied and analysed before concluding about the possible universality of the exponent $\alpha$, but our initial results do not point towards a constant value. Once more, it would be very valuable to obtain data in $d=4$ to see if a different value for $\alpha$ is found in larger spatial dimensions.

III.3 Breakdown of the Adam-Gibbs relation and experimental estimation of $\alpha$

Before starting this study, we felt that there was a general consensus in the community that the Adam-Gibbs relation is well-obeyed in real materials analysed near the experimental glass transition $T_{g}$. Thus, the outcome of the computer simulations showing deviations from Eq. (1) appeared as a worrying disagreement between simulations and experiments.

Therefore, we decided to collect data sets for several molecular liquids, where high-precision dynamic and thermodynamic data would be available over both simulation and experimental time windows, in order to perform a direct comparison with computer models.

We present the results of our data collection in Fig. 6(a) using again the representation where the standard Adam-Gibbs relation would yield a straight line. We use the fitted functional form for $S_{\rm conf}$ (obtained in Fig. 3(a)), whereas the actual data points are used for $\tau_{\alpha}$. When analysed over the entire experimental time window, defined above, we again observe a clear concavity for most materials. The Adam-Gibbs relation in Eq. (1) is violated over this regime, although of course it holds if observed over a restricted time window close to $T_{g}$ Richert and Angell (1998) (almost by definition–the data is continuous!).

As for the simulations, we fit the experimental data using the exponent $\alpha$ as an additional free parameter. From the experimental data, we determine two distinct values for $\alpha$, obtained by fitting either over the simulation or the experimental time window. The typical trend that we observe is that $\alpha>1$ over the simulation time window, but $\alpha<1$ over the experimental time window. The latter fits are included in Fig. 6(a), and they describe well the data over the entire experimental time window. As with the case for the simulation models, in Fig. 6(b), we also present the generalised Adam-Gibbs plot with the fitted $\alpha$ value for each material in the experimental time window. We confirm that the linear behavior is recovered in this plot.

We notice that the concavity in the Adam-Gibbs plot in the experimental time window was already reported Ngai (1999); Roland et al. (2004). However, the concavity would be overlooked as it is less pronounced than the convexity found at much higher temperature, close to $T_{\rm mct}$ and above Ngai (1999). Moreover, Ref. Roland et al. (2004) concluded that the observed concavity was attributed to an imprecise estimate of the configurational entropy. Our results obtained from simulation data with accurate configurational entropy measurements and recent high-quality experimental data suggest instead that the observed concavity is a generic physical phenomenon reflecting the nature of glassy dynamics over the experimental time window.

IV Discussion

Our central conclusion from both simulations and experiments considered over a broad time regime $\tau_{\alpha}/\tau_{o}\in[10^{3},10^{12}]$ (defined to be both experimentally accessible and theoretically relevant) is that the conventional Adam-Gibbs relation in Eq. (1) is not obeyed. Instead, the general form predicted by RFOT theory in Eq. (4) describes numerical and experimental data well. This is maybe not so surprising, from an empirical viewpoint, given that the generalised relation has one more free fitting parameter.

We compile all our results for the values of $\alpha$ from simulations (empty points) and experiments (filled points) in Fig. 7. To organise the data, we use the kinetic fragility index $m$ as the horizontal axis. This is simply a matter of convenience (as a matter of fact, no strong trend is observed). Note that, somewhat paradoxically, we do not have values for $\alpha$ in the computer models over the simulation time window because our computational schemes to measure $S_{\rm conf}$ only become applicable for low enough temperatures, typically $T\lesssim T_{\rm mct}$ Berthier and Coslovich (2014); Ozawa et al. (2018).

The experimental data in Fig. 7 obtained by considering the simulation time window are dispersed, $\alpha=0.73-2.34$, and tend to be characterised by rather large values $\alpha>1$. By contrast, considering a broader and physically better justified experimental time window, data for both simulations and experiments are much less scattered, $\alpha\simeq 0.24-1.28$, with a preferred average value $\alpha\simeq 0.5-0.6$, except for ethanol.

Before concluding, we make a further caveat regarding the above analysis of the RFOT theory predictions. In principle, we could have introduced additional subdominant physical prefactors into the scaling relations in Eqs. (2, 3) that could also be temperature dependent quantities. In particular, a surface tension could enter the relation between $S_{\rm conf}$ and $\xi_{\rm pts}$ Biroli and Bouchaud (2012); Lubchenko and Rabochiy (2014), and an energy scale could enter the activated scaling relation in Eq. (3). These prefactors would become irrelevant if some asymptotic regime could be reached with extremely long relaxation times and very small configurational entropy values, but it is understood that experimental glasses are not in this regime Tarjus (2011). In the absence of strong theoretical insights into these quantities, we decided to ignore them. They could of course very well affect the measured values of the reported exponents. Thus, a better determination of these quantities is an important research goal Cammarota et al. (2009a, b); Ganapathi et al. (2018), in particular in the experimental time window.

We also discuss potential sources of uncertainty in terms of experimental measurements of $\tau_{\alpha}$ and $S_{\rm conf}$, whose accuracy would affect the determination of the scaling exponent $\alpha$. Regarding $\tau_{\alpha}$, we note that the dielectric relaxation measurement for ethanol involves a Debye relaxation process which is distinct from the structural $\alpha$ relaxation process, as recently clarified in an experiment Chua et al. (2017). Indeed, the relaxation time extracted from the main peak that we used in this paper Brand et al. (2000) corresponds to the former process in ethanol whereas instead we should use the $\alpha$ relaxation process, but it is also found that the overall temperature dependence of two relaxation processes are very similar Chua et al. (2017). It could be that the unusual behavior in ethanol, showing $\alpha>1$, is related to this issue.

The experimental measurement of $S_{\rm conf}$ also involves approximations. First, using the excess entropy, $S_{\rm exe}=S_{\rm liq}-S_{\rm cry}$, instead of $S_{\rm conf}=S_{\rm liq}-S_{\rm glass}$ is an approximation, in general. The validity of $S_{\rm glass}\approx S_{\rm cry}$ has been widely studied Goldstein (1976); Yamamuro et al. (1998); Johari (2000); Martinez and Angell (2001); Angell and Borick (2002); Smith et al. (2017); Alvarez-Donado and Antonelli (2019); Han et al. (2019) and its validity seems to be non-universal Smith et al. (2017). Typically, $S_{\rm glass}$ is determined from the heat capacity of the non-equilibrium glass state, and so it still involves some approximations compared to its theoretical definition Yoshimori and Odagaki (2011).

Second, the measurements of $S_{\rm conf}$ in Ref. Tatsumi et al., 2012 were performed by doping the different materials to avoid crystallization. In particular, measurements for toluene and ETB involve $10$ wt % doping with benzene. Therefore, mixing effects possibly contribute to the absolute value of $S_{\rm conf}$, and to its temperature dependence Ozawa and Berthier (2017).

To summarize our results in terms of numerical values for the critical exponents introduced within RFOT theory, we observe in $d=3$ that the combination $\theta\simeq 3/2$ and $\alpha\simeq 0.5-0.6$ works well, which would then result in $\psi$ falling in the range $\psi\simeq 0.75-0.90$. If we use instead the value $\theta=2$, we would obtain a somewhat larger value for the dynamic exponent $\psi\simeq 1.0-1.2$, which agrees well with earlier indirect analysis Capaccioli et al. (2008); Brun et al. (2012). Both values violate the general bound $\psi\geq\theta$ discussed in the context of spin glasses Fisher and Huse (1988), the equality $\psi=\theta$ found for the random field Ising model Balog and Tarjus (2015), and the prediction $\psi=\theta=d/2$ in Ref. Kirkpatrick et al. (1989). In the absence of stronger theoretical constraints, we tentatively conclude that the measured $\psi$ value that we observe appears somewhat small, i.e., smaller than all known theoretical predictions. In $d=2$, we get $\theta\simeq 1.1$ and $\alpha\simeq 1.9$, which in turns implies that $\psi\simeq 1.7$, which appears somewhat large, by contrast with $d=3$.

Our conclusion that $\alpha<1$ is favored by the data over the experimental time window sheds some new light on an old debate in the glass literature Tanaka (2003); Lubchenko and Wolynes (2007); Hecksher et al. (2008); Elmatad et al. (2010). Assuming the existence of an ideal glass transition at equilibrium where $S_{\rm conf}\to 0$ and $\tau_{\alpha}\to\infty$, one is naturally led to the determination of two critical temperatures: the Kauzmann temperature $T_{K}$ where $S_{\rm conf}$ vanishes, and the critical temperature $T_{0}$ where the relaxation time diverges (not to be confused with onset temperature $T_{o}$ used above). Typically, the latter is obtained from a Vogel-Fulcher-Tammann fit ($T_{0}=T_{\rm VFT}$ in Eq. (8)) to the relaxation time. The possible equality $T_{0}=T_{K}$ would provide a strong empirical sign for the existence of an ideal glass transition underlying glass formation Lubchenko and Wolynes (2007). A large data set collected by Tanaka suggests the existence of systematic differences between the two temperatures Tanaka (2003), with the tendency that $T_{K}>T_{0}$, and an apparent correlation with kinetic fragility. In our analysis using Eq. (4) to describe the data, the connection between thermodynamics and dynamics becomes automatically satisfied, and thus by construction thermodynamic and dynamic singularities necessarily coincide. Assuming that the determination of $T_{K}$ is the most robust one, we conclude that it is the experimental determination of $T_{0}$ which should be questioned. In particular, using $\alpha<1$ in Eqs. (4) and assuming an asymptotically linear vanishing of $S_{\rm conf}$, one would predict that $\log(\tau_{\alpha}/\tau_{0})\propto(T-T_{0})^{-\alpha}$, which is distinct from the standard Vogel-Fulcher-Tamman fit and would automatically produce the equality $T_{K}=T_{0}$.

From a broader perspective, we conclude that the Adam-Gibbs relation, which is an important milestone in the field of glass transition studies, is generally violated in both computer models and real materials when tested over a broad, experimentally-relevant temperature range. We nevertheless argued that the failure of Eq. (1) cannot be taken as evidence that thermodynamic theories of the glass transition are incorrect. The RFOT theory prediction of a connexion between statics and dynamics in Eq. (4) is obeyed by all materials, with exponent values that are reasonable, but remain to be predicted from first principles. A larger concern, perhaps, is the apparent lack of universality in the data shown in Fig. 7 which clearly display variations from one system to another. This may still be rationalised by invoking the fact that $\alpha$ is obtained from the analysis of a finite time window where additional preasymptotic effects and temperature dependent prefactors may influence the reported results.

Taking an orthogonal perspective, we finally ask: Do our results validate or invalidate some theories of the glass transition? After all, we just established that a slightly generalised version of the Adam-Gibbs relation with $\alpha\simeq 0.6$ describes simulations and experiments over 9 orders of magnitude in the experimentally relevant regime. This is not a small accomplishment. One can take the alternative view that the deviations from the canonical exponent values should be taken as an indirect sign that thermodynamics only contributes some part of the slowing down, in addition to other physical factors Tarjus et al. (2005); Rabochiy et al. (2013); Wyart and Cates (2017); Dyre (2006); Tanaka (2012); Ikeda et al. (2017). This view is sometimes also invoked to rationalise the “modest” growth of static correlation lengthscale observed numerically and experimentally Tarjus (2011); Yaida et al. (2016). Our finding that $\alpha<1$ suggests instead that it is the growth of the relaxation time that is actually too modest! It is therefore difficult to rationalise how another physical factor working in addition to the entropy could be invoked to explain our findings. The most radical view is in fact that thermodynamics is just a spectator to the glassy dynamics Chandler and Garrahan (2010), in which case our findings should be interpreted as purely coincidental since entropy plays in fact no role. We have no strong argument to oppose to this view, which remains perfectly admissible.

Acknowledgements.

We thank J.-P. Bouchaud, D. Coslovich and F. Zamponi for insightful discussions. We also thank S. Tatsumi and O. Yamamuro for sharing their high-quality experimental data for the configurational entropy. The research leading to these results has received funding from the Simons Foundation (#454933, Ludovic Berthier).

Appendix A Configurational entropy along an isobaric path

We wish to measure the configurational entropy $S_{\rm conf}(T,P)$ along an isobaric (constant pressure) path. It is computed as $S_{\rm conf}(T,P)=S_{\rm tot}(T,P)-S_{\rm glass}(T,P)$, where $S_{\rm tot}(T,P)$ and $S_{\rm glass}(T,P)$ are the total and glass entropies at the temperature $T$ and pressure $P$. We explain how to get $S_{\rm conf}(T,P)$ from $NPT$ simulation trajectories.

A.1 Notations

We consider the Helmholtz free energy $-\beta F(T,V)=\ln Z(T,V)$, where $\beta=1/T$ and $Z(T,V)$ is the partition function of the $NVT$ ensemble. We also consider the Gibbs free energy $-\beta G(T,P)=\ln Y(T,P)$, where $Y(T,P)$ is the partition function of the $NPT$ ensemble, given by

[TABLE]

We introduce the probability distribution of the volume $V$ for a given $T$ and $P$,

[TABLE]

In equilibrium, $\rho(V|T,P)$ is given by Gaussian distribution,

[TABLE]

where $V_{*}$ and $\sigma_{V}^{2}$ are the mean and variance of the volume, respectively. We define $\langle(\cdots)\rangle_{T,P}=\int_{0}^{\infty}\mathrm{d}V\rho(V|T,P)(\cdots)$. Using this average, we can write $V_{*}=\langle V\rangle_{T,P}$ and $\sigma_{V}^{2}=\langle(V-V_{*})^{2}\rangle_{T,P}$.

A.2 Total entropy

The total entropy $S_{\rm tot}(T,P)$ is obtained by a thermodynamic integration of the isobaric heat capacity from a reference temperature $T_{\rm ref}=1/\beta_{\rm ref}$, to the target temperature $T=1/\beta$,

[TABLE]

where $U_{*}(T,P)$ is the mean potential energy, and $V_{*}(T,P)$ is the mean volume; $U_{*}(T,P)$ and $V_{*}(T,P)$ are measured by constant pressure simulations. The entropy at the reference state is obtained by $S_{\rm tot}(T_{\rm ref},P)=\langle S_{\rm tot}(T_{\rm ref},V)\rangle_{T,P}$ using the $NVT$ ensemble scheme Berthier et al. (2017). This treatment for the reference state will be justified below.

A.3 Glass entropy

To get the glass entropy, we use the generalised Frenkel-Ladd method which relies on the $NVT$ ensemble Ozawa et al. (2018). In general, one can smoothly connect $NVT$ and $NPT$ ensembles in terms of mean values. For example, thermodynamics guarantees that $S(T,P)=S(T,\langle V\rangle_{T,P})$. However, special attention should be paid if one uses the $NVT$ ensemble scheme with trajectories generated by the $NPT$ ensemble for finite system size Lebowitz et al. (1967). A related issue is discussed in Ref. Cheng and Ceriotti (2018). Indeed, what we can compute is $\langle S(T,V)\rangle_{T,P}$. In general,

[TABLE]

Therefore, we need to consider the second term in Eq. (16) as a correction term. We can evalute this term with Eq. (14):

[TABLE]

Since $\sigma_{V}^{2}\sim N$, this term vanishes in the thermodynamic limit, as expected. Indeed, for $N=1500$ systems, we get negligible values, $\frac{1}{N}\ln\sqrt{2\pi e\sigma_{V}^{2}}\simeq 0.0026$ and $0.0013$ at $T_{\rm ref}=7.0$ and $T=0.37$, respectively. These values are small compared to the absolute value of $S_{\rm conf}/N\simeq 0.36-0.80$. Thus we can safely use $S(T,P)=\langle S(T,V)\rangle_{T,P}$. Especially we use the following equation, $S_{\rm glass}(T,P)=\langle S_{\rm glass}(T,V)\rangle_{T,P}$.

Appendix B Extrapolation of relaxation times towards $T_{g}$

Here we test the validity of the extrapolation of relaxation time from the numerical to the experimental timescale using various fitting functions. We employ 1-propanol, propylene glycol, glycerol, OTP, and PC. Among these, 1-propanol, propylene glycol, and glycerol have kinetic fragility indexes similar to the simulation models.

Figure 8 shows various fits of the data performed over the simulation time window, $\tau_{\alpha}\leq 10^{-5}$ s, and then extrapolated to lower temperatures down to $T_{g}$ where $\tau_{\alpha}=100$ s. In the cases for 1-propanol, propylene glycol, glycerol, and OTP, shown in Fig. 8, the parabolic law is the best functional form that predicts the actual data well over the experimental time window. All other functional forms, when fitted over the simulation time window, tend to deviate from the actual data at low temperatures. For the most fragile material, PC, underestimates the actual data when the parabolic law is applied, whereas MYGEA and AM predict the data better. Notice that the uncertainty on the determination of $T_{g}$ using the numerical time window and a parabolic fit is very small for the systems whose fragility is comparable to typical simulations models. This is the strategy we have used in previous numerical studies Ninarello et al. (2017); Berthier et al. (2017, 2019b).