On mean-field theories of dynamics in supercooled liquids

TL;DR

This paper introduces a hybrid numerical method to precisely determine the memory function in supercooled liquids and compares it with existing mean-field theories, revealing both similarities and differences.

Contribution

It develops a hybrid numerical approach to extract the exact memory function and compares it with two mean-field models, providing insights into their accuracy.

Findings

The exact memory function shares traits with mean-field models.

Quantitative differences highlight limitations of mean-field approximations.

The approach aids in understanding dynamics of supercooled liquids.

Abstract

We develop a hybrid numerical approach to extract the exact memory function K(t) of a tagged particle in three-dimensional glass-forming liquids. We compare the behavior of the exact memory kernel to two mean-field approaches, namely the standard mode-coupling theory and a recently proposed ansatz for the memory function that forms the basis of a new derivation of the exact form of K(t) for a fluid with short-ranged interactions in infinite dimensions. Each of the mean-field functions qualitatively and quantitatively share traits with the exact K(t), although several important quantitative differences are manifest.

Figures

Tables

Table 1. Table 1: Temperatures used in our simulations. For each temperature, data was averaged over n t subscript 𝑛 t n_{\mathrm{t}} trajectories.

$T$	5.0	2.0	1.0	0.8	0.7	0.6	0.55	0.52	0.49	0.47	0.46	0.45
$n_t$	2000	1850	1750	1000	980	683	256	339	500	500	497	454

Equations

$$U_{ij}(r)=4{ϵ}_{ij}\left[{\left(\frac{{\sigma }_{ij}}{r}\right)}^{12}-{\left(\frac{{\sigma }_{ij}}{r}\right)}^6\right],$$

$$S(r)=\{\begin{array}{l}1,\text{if }r\le r_{\mathrm{on}}\\ \frac{(r_{\mathrm{cut}}^2-r^2{)}^2(r_{\mathrm{cut}}^2+2r^2-3r_{\mathrm{on}}^2)}{(r_{\mathrm{cut}}^2-r_{\mathrm{on}}^2{)}^3},\text{if }r>r_{\mathrm{on}}>r_{\mathrm{cut}}\\ 0,\text{if }r\ge r_{\mathrm{cut}},\end{array}$$

$$C^P(t)$$

$$C^{FP}(t)$$

$$B(t)=e^{i\mathcal{L}t}B(0),\dot{B}(t)=i\mathcal{L}B(t),$$

$$i\mathcal{L}=\sum _{i=1}^N\left[\frac{\partial H}{\partial P_i}\frac{\partial }{\partial Q_i}-\frac{\partial H}{\partial Q_i}\frac{\partial }{\partial P_i}\right],$$

$$\mathcal{P}B=\frac{⟨P(0)B⟩}{⟨P(0{)}^2⟩}P(0)=\frac{⟨P(0)B⟩}{k_{\mathrm{B}}T}P(0),$$

$$\mathcal{Q}=1-\mathcal{P}.$$

$$\dot{P}(t)=-{\int }_0^{t}K(u)P(t-u)du+R(t),$$

$$K(t)=\frac{⟨F(0)e^{i\mathcal{Q}\mathcal{L}t}F(0)⟩}{T},$$

$${\dot{C}}^P(t)=-{\int }_0^{t}K(t-u)C^P(u)du.$$

$$e^{i\mathcal{Q}\mathcal{L}t}=e^{i\mathcal{L}t}-{\int }_0^t{e}^{i\mathcal{L}u}\mathcal{P}i\mathcal{L}{e}^{i\mathcal{Q}\mathcal{L}(t-u)}du,$$

$$K(t)=\frac{1}{T}\left[C^F(t)+{\int }_0^t{C}^{FP}(u)K(t-u)du\right].$$

$$C_{\mathrm{d}}(t)=\frac{1}{3NT}\sum _{j=1}^N\sum _{i>j}^N\sum _{\alpha =x,y,z}⟨F_{ji,\alpha }(0)F_{ji,\alpha }(t)⟩$$

$$K(t)\to C_{\mathrm{d}}(t),$$

$${\mathcal{P}}_2(B)=\sum _{\vec{k}}\frac{⟨\hat{b}(-\vec{k})B⟩}{NS(\vec{k})}\hat{b}(\vec{k}),$$

$$K_{\mathrm{MCT}}(t)=\frac{⟨[{\mathcal{P}}_{2}F(0)]e^{i\mathcal{L}\mathcal{Q}t}[{\mathcal{P}}_{2}F(0)]⟩}{T},$$

$$K_{\mathrm{MCT}}(t)=\frac{T}{8{\pi }^3\rho }\int d\vec{k}{k}_{\alpha }^{2}c(\vec{k}{)}^2{F}_{\mathrm{s}}(\vec{k},t)F_{\mathrm{c}}(\vec{k},t)$$

$$\hat{K}(s)=\frac{T-s{\hat{C}}^P(s)}{{\hat{C}}^P(s)},$$

$$C_{\mathrm{short}}(t)=\frac{a_1}{\cosh (a_{2}t)},$$

$$C_{\mathrm{long}}(t)=A\exp [(t/\tau {)}^{\beta }],$$

$${\zeta }^{(\mathrm{mem})}$$

$${\zeta }^{(\mathrm{diag})}$$

$$\hat{C}(s)={\int }_0^{\infty }{e}^{-st}C(t)dt,$$

$$C(t)=\frac{1}{2\pi i}\underset{\delta \to \infty }{\lim }{\int }_{\gamma -i\delta }^{\gamma +i\delta }{e}^{st}\hat{C}(s)ds,$$

$$C_n(t)=\frac{\log (2)}{t}\sum _{k=1}^{2M}{\omega }_k\hat{C}\left(\frac{k\log (2)}{t}\right),$$

$${\omega }_k$$

$$(-1{)}^{M+k}\sum _{j=\lfloor \frac{k+1}{2}\rfloor }^{\min (k,M)}\frac{j^M(2j)!}{(M-j)!j!(j-1)!(k-j)!(2j-k)!},$$

$$p_{\max }=\frac{2M\log (2)}{t_{\min }},$$

$$K_i=\frac{C_i^F}{T}+\frac{dt}{2T}\left(C_{i0}^{FP}{K}_0+2\sum _{j=1}^{i-1}{C}_{ij}^{FP}{K}_j+C_{ii}^{FP}{K}_i\right).$$

$$K_i=\frac{C_i^F}{T}+\frac{dt}{2T}\left(C_{i0}^{FP}{K}_0+2\sum _{j=1}^{i-1}{\mathcal{C}}_{ij}{K}_j\right).$$

$$K_i=\frac{C_i^F}{T}+\frac{dt}{2T}\left(C_{i0}^{FP}{K}_0+2\sum _{j=\max (0,i-i^{\ast })}^{i-1}{C}_{ij}^{FP}{K}_j\right).$$

Full text

I Introduction

Upon cooling, supercooled liquids (SCLs) display a dramatic increase in viscosity, the origin of which is still not fully understood Angell (1995); Berthier and Biroli (2011); Charbonneau et al. (2017). As the experimental glass transition is approached, dynamics become highly correlated and heterogeneous, making the construction of a first-principles microscopic theory a very challenging task. Historically, the first partially successful set of microscopically-based dynamical equations of SCLs was the mode-coupling theory (MCT) Götze (2008); Reichman and Charbonneau (2005); Janssen (2018), which takes as input only time-independent quantities and outputs dynamical predictions via a mean-field-like factorization closure of the memory kernel ${K}(t)$. The dynamical equations derived in this manner match experimental results in weakly supercooled liquids, but clearly deviate from them in the deeply supercooled regime. In particular MCT predicts a power-law divergence of relaxation times at a temperature $T_{\mathrm{MCT}}$, which is much higher than the value $T_{\mathrm{g}}$ of the empirically measured glass transition temperature Fuchs et al. (1998).

In the 1980s it was noticed that the dynamical equations for the MCT density-density correlation function have a form identical to that found for the spin correlation function in the $p=3$ $p$-spin model Kirkpatrick and Thirumalai (1987), a paradigmatic member of a family of mean-field spin-glass models. This observation formed the basis for the random first-order theory (RFOT) Kirkpatrick et al. (1989), which assumes that the equilibrium solution of the $p$-spin model, obtained within replica theory, also holds for supercooled liquids. A direct implication of this connection is that MCT is a dynamical mean-field approximation which should become exact in the limit of infinite spatial dimensions, $d\to\infty$. This assumption was challenged in 2010, when it was observed numerically that with increasing dimensionality the replica and the MCT approaches give completely divergent transition points Ikeda and Miyazaki (2010); Schmid and Schilling (2010).

More recently, light has been shed on this discrepancy through a mean-field replica-based first-principles calculation showing that the RFOT scenario is the exact theory for the thermodynamics of high-dimensional sphere packings Charbonneau et al. (2014). This was followed by a remarkable exact solution of the (mean-field) $d=\infty$ dynamics of simple SCLs Maimbourg et al. (2016), which confirmed explicitly that the original MCT behaves incorrectly in the $d\to\infty$ limit. It should be noted that the theory of Ref. Maimbourg et al. (2016), which is also fully consistent with the dynamics of $p$-spin models, takes a completely different form than that of standard MCT. In particular, the mean-field theory of Ref. Maimbourg et al. (2016) forms a closed theory on the dynamical trajectories, while MCT is closed in terms of correlation functions. Despite this fact, the exact $d=\infty$ solution does share many features with standard MCT.

In 2017 a transparent set of simple liquid state approximations was proposed which is capable of reproducing the correct $d\to\infty$ limit Szamel (2017). The main assumption, which holds in infinite dimensions, is that the memory function for tagged particle transport ${K}(t)$ can be approximated through the real-space diagonal part of the pair force correlation function, which we denote $C_{\mathrm{d}}(t)$. This approximation provides a direct link between the exact behavior of a fluid in $d=\infty$ (where this approximation becomes exact), and its three-dimensional counterpart. The validity of this approximation could be assessed by comparing ${K}(t)$ and $C_{\mathrm{d}}(t)$ in $3d$ numerical simulations. This is a challenging undertaking, however, because the memory function ${K}(t)$ cannot be easily calculated directly from simulated dynamics. In particular, the time dependence of ${K}(t)$ is generated with an unusual form of dynamics, and in the SCL regime the behavior of the memory function spans several orders of magnitude in time. As a consequence, any route to the numerical evaluation of ${K}(t)$ becomes a computationally daunting task.

In this paper, we use a combination of methods to accurately calculate the memory function ${K}(t)$ in a three-dimensional glass-forming liquid. Several features of ${K}(t)$ are compared with both $C_{\mathrm{d}}(t)$ and the more standard MCT approximation to $K(t)$ in three dimensions for the case of a canonical glass-forming liquid. We find clear quantitative differences between the exact $K(t)$ and these two functions. On the other hand, salient qualitative features are in agreement, implying that there are important features of liquids in $d=3$ that owe their properties to the mean-field, $d=\infty$ behavior. Such insights pave the way for more sophisticated approximations in three dimensions that could potentially quantitatively capture the collective behavior of liquids near the experimental glass transition temperature.

In Sec. II, we describe the system we simulate, and illustrate some important dynamical observables that result from our molecular dynamics simulations. In Sec. III we show how an exact self-consistent equation for the memory function ${K}(t)$ can be derived through the projector operator technique, and we discuss the main assumptions of Ref. Szamel (2017). We then outline the methods used to calculate ${K}(t)$ in Sec. IV, and present our results in Sec. V. Finally, we summarize and discuss our findings in Sec. VI.

II Model and simulations

We simulate a three-dimensional Kob-Andersen 80:20 binary mixture Kob and Andersen (1994) of 1080 particles with molecular dynamics 111All the code used for this article is available at github.com/mbaityje/STRUCTURAL-GLASS., for which no finite-size effects in the autocorrelation times are expected or observed Berthier et al. (2012). We work in Lennard-Jones units Allen and Tildesley (1989), with an integration time step of $dt=0.0025$, and a mass of all particles set to $m=1$. The box is a periodic cube of linear size $L=(N/\rho)^{1/3}\approx 9.65$, chosen so that the density $\rho=N/L^{3}=1.2$. The interparticle potential at distance $r$ is

[TABLE]

where $i,j\in\{\mathrm{A},\mathrm{B}\}$ indicate the particle type. The interaction energy amplitudes are $\epsilon_{\mathrm{AA}}=1,\epsilon_{\mathrm{AB}}=1.5$ and $\epsilon_{\mathrm{BB}}=0.5$, and the non-additive particle radii, chosen in order to suppress crystallization, are $\sigma_{\mathrm{AA}}=1,\sigma_{\mathrm{AB}}=0.8$ and $\sigma_{\mathrm{BB}}=0.88$. The potential is cut off and smoothed in order to have a continuous force at any distance 222 For our molecular dynamics simulations we have used the hoomd package, available at (http://glotzerlab.engin.umich.edu/hoomd-blue). The package is a Python API for GPU-based molecular dynamics of liquids, developed in Refs. Anderson et al. (2008); Glaser et al. (2015). The smoothing we applied follows the xplor mode described in the hoomd documentation. The smoothed potential is $U_{ij}^{\mathrm{smooth}}(r)=S(r)U_{ij}(r)$, where

$S(r)=\begin{cases}1,~{}~{}~{}~{}~{}\text{if }r\leq r_{\mathrm{on}}\\ \frac{(r_{\mathrm{cut}}^{2}-r^{2})^{2}(r_{\mathrm{cut}}^{2}+2r^{2}-3r_{\mathrm{on}}^{2})}{(r_{\mathrm{cut}}^{2}-r_{\mathrm{on}}^{2})^{3}},~{}\text{if }r>r_{\mathrm{on}}>r_{\mathrm{cut}}\\ 0,~{}~{}~{}~{}~{}\text{if }r\geq r_{\mathrm{cut}}\,,\end{cases}$

where we used $r_{\mathrm{on}}=1.2\sigma$ and $r_{\mathrm{cut}}=2.5\sigma$. The code used for these simulations is open-source and available at www.github.io/mbaityje/STRUCTURAL-GLASS. . For this model, there is an expected dynamical crossover at a temperature $T_{\mathrm{d}}\approx 0.435$ 333See Ref. Kob and Andersen (1994) for an estimation of the mode-coupling crossover. The potential used in Ref. Kob and Andersen (1994) is shifted but not smoothed. This difference can induce a small shift in the value of the mode-coupling temperature in comparison to that of our system..

The temperatures $T$ employed in our simulations are given in Table 1. For each $T$ we have thermalized 10 independent systems, and from each we run $n_{\mathrm{traj}}$ uncorrelated trajectories from which the autocorrelation functions detailed below are calculated. Observables are measured with a quasi-exponentially decreasing rate 444For each trajectory, we used an exponential succession of 1000 time points where measurements were taken. Since time steps are discrete, some times are duplicates and the grid results in having fewer points..

In Fig. 1a we show the self-intermediate scattering function $F_{\mathrm{s}}(\vec{k},t)=\frac{1}{N}\langle\sum_{j=1}^{N}e^{i\vec{k}\cdot(\vec{r}_{j}(t)-\vec{r}_{j}(0))}\rangle$, with $\vec{k}$ denoting the wave vector, and $\vec{r}_{j}(t)$ the position of particle $j$ at time $t$. Since here and throughout this paper all particles have mass $m=1$, particle velocities are equal to momenta, $P(t)$, and accelerations are equal to forces, $F(t)$. The autocorrelation functions of these observables, defined as

[TABLE]

are shown in Fig. 1b-d. Note that $\dot{C}^{P}=C^{PF}(t)$, and $\dot{C}^{FP}=C^{F}(t)$, where we use the over-dot to indicate time derivatives.

In Fig. 2 we show the mean-square displacement of the particles, $\varDelta^{2}(t)=\left\langle\frac{1}{N}\sum_{i=1}^{N}\left[\vec{r}(t)-\vec{r}(0)\right]^{2}\right\rangle$, which for long times is related to the autocorrelation functions via $\varDelta^{2}(t)=2t\int_{0}^{t}C^{P}(u)du\equiv 6tD$ Balucani and Zoppi (1995). From the long-time behavior of $\varDelta^{2}(t)$, one can extract the diffusion coefficient, $D=\lim_{t\to\infty}\frac{\varDelta^{2}(t)}{6t}$ 555 In the extrapolation, we took into account $\mathcal{O}(t^{-1})$ corrections to scaling. Thus, the fitting relation is $\frac{\varDelta^{2}(t)}{6t}=D+\kappa/t$, where $D$ and $\kappa$ are the fitting parameters.. If we fit the vanishing of the diffusion coefficient to a power-law form $f_{D}(T)=a_{D}(T-T_{D})^{\eta_{D}}$, with $T\leq 0.7$, we obtain $T_{D}=0.419(4)$, $\eta_{D}=2.13(7)$. These values should be taken cum grano salis, since both exponent and $T_{D}$ depend on the fitting range Flenner and Szamel (2005). However, the values we extract are grossly consistent with those previously reported in the literature Kob and Andersen (1994); Flenner and Szamel (2005).

III The tagged particle memory function

In a system defined by a Hamiltonian $H$, a generic observable $B(0)$ at time $t=0$ exhibits a formal time dependence for $t>0$

[TABLE]

where the Liouville operator $\mathcal{L}$ is defined via

[TABLE]

where $P_{i}$ is the momentum of particle $i$, and $Q_{i}$ is its position.

The Mori-Zwanzig projection operator technique is based on projecting the dynamics onto a set of relevant observables (here, the impulse of a tagged particle), leaving the remaining degrees of freedom to evolve in an orthogonal dynamical subspace Mori (1965); Zwanzig (2001); Balucani and Zoppi (1995). We define the projection operator

[TABLE]

where $P$ is one component of the impulse of a tagged particle. Throughout this work we set the Boltzmann constant to $k_{\mathrm{B}}=1$. The orthogonal projector is defined as

[TABLE]

With the use of these projection operators the evolution of the system is split into two complementary sets of orthogonal dynamics, $i\mathcal{L}=i\mathcal{PL}+i\mathcal{QL}$, one projected onto the relevant variable(s) and one orthogonal to it (them).

One can show with the aid of these projectors that the momentum of a tagged particle follows the exact generalized Langevin equation Balucani and Zoppi (1995)

[TABLE]

where $R(t)=e^{i\mathcal{QL}t}\mathcal{Q}P(0)$ is usually called the fluctuating force, i.e. the force arising from the components of the system orthogonal to the relevant variable, which in this formalism may be thought of as a source of stochastic noise. The function ${K}(t)$ is the memory function, which takes the form of a force-force correlation function evolving with projected dynamics,

[TABLE]

where $F(0)$ is a component of the force acting on the tagged particle at time $t=0$. ${K}(t)$ is thus the time autocorrelation function of the fluctuating force. Eq. (9) is deceptively simple, since robust methods to compute exact projected dynamics are scarce.

The definition of the projector in Eq. (6) renders the fluctuating force $R(t)$ orthogonal to the momentum at all times, so Eq. (8) can be used to express the evolution of the momentum-momentum correlation function in terms of the memory function as

[TABLE]

The memory function ${K}(t)$ can be written in a manner that does not involve orthogonal dynamics, by using in Eq. (9) the Dyson relation Shi and Geva (2003),

[TABLE]

and by noticing that $\mathcal{P}(-i\mathcal{L})e^{i\mathcal{QL}(t-u)}F(0)=\langle F(0)e^{i\mathcal{QL}(t-u)}F(0)\rangle P(0)/\langle P(0)^{2}\rangle$, which leads to an integral equation for ${K}(t)$,

[TABLE]

Eq. (12) has the advantage that the input to the integral equation, namely $C^{F}(t)$ and $C^{FP}(t)$, involve normal unprojected dynamics. In principle, the solution of the integral equation (12), generated with the exact input of $C^{F}(t)$ and $C^{FP}(t)$, provides an exact means to obtain ${K}(t)$.

In Ref. Szamel (2017) a physically-motivated set of approximations was proposed that provides a route to the the exact solution of the dynamics in $d=\infty$ for systems with short-ranged potentials Maimbourg et al. (2016). The principal assumption of Ref. Szamel (2017) is that a mean-field theory for the memory function can be obtained by discarding the components of the force-force autocorrelation function that involve more than two particles, and simultaneously replacing the projected dynamics with normal dynamics. Explicitly, this means the replacement of ${K}(t)$ with the diagonal force correlator,

[TABLE]

where we choose the tagged particle to be particle 1, and $F_{ji,\alpha}(t)=e^{i\mathcal{L}t}F_{ji,\alpha}(0)~{}[\alpha=x,y,z]$ is a component of the force that particle $i$ exerts on particle $j$ at time $t$. Clearly, $C_{\mathrm{d}}(t)$ is the force-force correlation $C^{F}(t)$ where all the non-diagonal terms which couple different particles at different times are discarded. This function is expected to decay slower than $C^{F}(t)$, due to cancellations between the diagonal and non-diagonal components of $C^{F}(t)$ Vergeles and Szamel (1999).

While the substitution proposed in Ref. Szamel (2017),

[TABLE]

is a key step for the recovery of the exact dynamics in $d=\infty$, it may be viewed as an interesting, albeit uncontrolled, mean-field-like approximation in $d=3$. It should be noted that this mean-field relationship leads to a picture distinct from that of the canonical mean-field theory of supercooled liquids, namely MCT, which is based on an uncontrolled factorization of density modes in $k$-space Götze (2008).

To calculate the tagged particle memory function within canonical MCT, one can project the forces in Eq. (9) on the pair density modes, which are expected to give the dominant contribution to the tagged particle dynamics Balucani and Zoppi (1995). In practice, this amounts to projecting the forces onto the wave-vector dependent solute self-density and solvent collective density, $\hat{b}(\vec{k})=\sum_{i=2}^{N}e^{i\vec{k}\cdot(\vec{r}_{i}-\vec{r}_{1})}$, through the projection operator Balucani and Zoppi (1995); Egorov et al. (2002)

[TABLE]

where $S(\vec{k})=\langle\frac{1}{N}\sum_{i,j=1}^{N}e^{i\vec{k}\cdot(\vec{r}_{i}-\vec{r}_{j})}\rangle$ is the static structure factor. The MCT memory function involving these modes is then expressed as

[TABLE]

where $F(0)$ is the $\alpha^{\mathrm{th}}$ component of the force at time $t=0$. By replacing Eq. (15) into Eq. (16) one obtains Balucani and Zoppi (1995); Egorov et al. (2002)

[TABLE]

where $F_{\mathrm{c}}(\vec{k},t)=\langle\sum_{j,\ell=1}^{N}e^{i\vec{k}\cdot(\vec{r}_{j}(0)-\vec{r}_{\ell}(t))}\rangle$ is the collective scattering function and $c(\vec{k})=\frac{1}{\rho}\left(1-\frac{1}{S(\vec{k})}\right)$. It is important to remark that the short-time behavior of ${K}(t)$ cannot be captured by Eq.(17). In particular, an integration over all $k$ will result in the near complete suppression of the slow dynamical regime. The usual procedure that is followed is that the short time kinetic behavior is subtracted and reintroduced phenomenologically via the addition of a function such as that in Eq. (19) Balucani and Zoppi (1995). However, we find that such a procedure is not consistent since the true amplitude of ${K}(0)$ is significantly smaller than that of ${{K}}_{\mathrm{MCT}}(0)$. In the following, we choose a cutoff value of $k$, $k_{\mathrm{max}}=28$, for which the long-time behavior of ${{K}}_{\mathrm{MCT}}(t)$ is converged and stable and for which a reasonable description of the short time behavior is captured.

Our goal for the remainder of this work will be to compare $C_{\mathrm{d}}(t)$ and ${{K}}_{\mathrm{MCT}}(t)$ to the exact ${K}(t)$ as extracted from molecular dynamics simulations. We note that our comparisons will differ in a subtle manner from a full comparison of either the mean-field theory of Ref. Maimbourg et al. (2016) or of MCT Götze (2008) to the exact behavior of ${K}(t)$ in $d=3$. In particular, here we use the exact molecular dynamics in $d=3$ to evaluate the respective mean-field expressions as opposed to solving the non-linear mean-field equations self-consistently. Thus, the comparisons we make will deviate somewhat from those associated with complete theories as outlined in Refs. Maimbourg et al. (2016); Szamel (2017) and Ref. Voigtmann et al. (2004), respectively. We will return to this point later in our work and in the conclusions. This caveat aside, we do expect the comparisons to be revealing of the gross successes and failures of the various mean-field approaches in low spatial dimensions.

IV Numerical evaluation of the memory function

A numerical evaluation of the tagged particle memory function ${K}(t)$ is challenging due to the fact that the dynamical evolution of the force is generated by projected dynamics. We have found that recent proposals to directly generate a molecular dynamics for projected dynamics Carof et al. (2014) are successful at high temperatures where the memory function is simple in structure and decays quickly in time, but the resulting equations become very stiff in dense fluids at low temperatures, rendering the procedure unstable.

To deal with the difficult task of evaluating the memory function over the full span of distinct dynamical regimes, we have adopted a piecewise concatenation of two separate methods, one of which is accurate at long times, the other of which is accurate at short times. For the evaluation of the long-time behavior of ${K}(t)$, we Laplace transform Eq. (10), obtaining

[TABLE]

where $\hat{A}(s)$ indicates the Laplace transform of the function $A(t)$, and we notice that $\hat{C}^{P}(0)=\int_{0}^{\infty}\langle P(0)P(t)\rangle=T$. We then invert the Laplace transform with the Gaver-Stehfest method Gaver Jr (1966); Stehfest (1970). This procedure, which we will call the Laplace method in the following, is capable of estimating ${K}(t)$ at long times, although it requires averaging over a large number of trajectories. In addition there are a series of caveats related to the precision of the approach which we discuss in Appendix A.

For the short-time behavior of ${K}(t)$, we note that Eq. (12) is a Volterra integral equation of the second kind. The integral on the right-hand side can be decomposed with a trapezoidal integration scheme, and solved in $\mathcal{O}(t^{2})$, by passing the last element of the integral to the left-hand side of Eq. (12) Press et al. (2007). This algorithm, which we will call the Volterra method, is very accurate for short times. Details are given in Appendix B.

The memory functions ${K}(t)$ that we report are a piecewise concatenation of the ${K}(t)$ obtained through the Volterra method for $t\lesssim 0.05$, and the Laplace method for $t\gtrsim 0.05$. As shown in Appendix C, the two methods give similar solutions around $t=0.05$, and the resulting ${K}(t)$ describes appropriately the dynamics, satisfying the consistency check of Sec. IV.1 below. Due to the instability of the inverse Laplace transform (see Appendix A), we had to average our data over a large number $n_{t}$ of trajectories (see Tab. 1).

IV.1 Consistency check on the obtained memory functions

To ensure that the obtained ${K}(t)$ is actually a faithful representation of the memory function ${K}(t)$, we calculate $C^{FP}(t)$ via Eq. (10) after obtaining ${K}(t)$ with the approach outlined above, and we compare it with the one we measured directly from molecular dynamics. We then integrate the obtained $C^{FP}(t)$ over time with a trapezoidal integration scheme, and compare it to the measured $C^{P}(t)$, since $C^{P}(t)$ contains more visually recognizable features than does $C^{FP}(t)$. In Appendix C we show the results of this consistency check.

V Results

In Fig. 3 we show the tagged particle memory function ${K}(t)$, the diagonal force autocorrelation function $C_{\mathrm{d}}(t)$ and the memory function calculated with a traditional MCT-type theory, ${{K}}_{\mathrm{MCT}}(t)$, at different temperatures. For better visual clarity, in Fig. 4 we show the three functions at two extreme temperatures. As discussed in Sec. III, the expression of Eq. (17) only provides a description of the long-time behavior of ${K}(t)$. In particular, neither the rapid decay of ${K}(t)$ nor the relative amplitude of the plateau height compared to $K(0)$ may be obtained from it. On the other hand, $C_{\mathrm{d}}(t)$ does provide an approximate and consistent description in both regimes. In a qualitative sense, the shape of the plateau region of ${K}(t)$ is also more accurately described by $C_{\mathrm{d}}(t)$ than by ${{K}}_{\mathrm{MCT}}(t)$, including the stability of the plateau and the dip in ${K}(t)$ close to $t=0.1$, which separates the short and long time regimes.

As is clear from Fig. 3 and as noted above, the short-time behaviors of ${K}(t)$ and $C_{\mathrm{d}}(t)$ are quantitatively similar, both in terms of the their respective rapid decreases as well as the appearance of an interesting (almost temperature-independent) oscillation that separates the short and long time regimes. To quantify the similarities between ${K}(t)$ and $C_{\mathrm{d}}(t)$, we fit the short-time behavior to Balucani and Zoppi (1995); Kob et al. (2002)

[TABLE]

where $a_{1}=1$ when the functions are normalized, so the only free parameter is $a_{2}(T)$, which we show in Fig. 5a. Clearly, even the mild temperature dependence of the rapid initial decay of ${K}(t)$ is captured bu the mean-field diagonal approximation $C_{\mathrm{d}}(t)$.

The long-time decay of the memory functions can be fit to a stretched exponential form,

[TABLE]

where the plateau height $A$, the autocorrelation time $\tau$, and the stretching exponent $\beta$ are free parameters.

In Fig. 5b we see that, although the plateau height for $C_{\mathrm{d}}(t)$ is larger than for ${K}(t)$, the former does not depend on temperature, whereas the latter grows mildly as $T$ decreases. The stretching exponent $\beta$ (Fig. 5c) is different for the two functions: ${K}(t)$ has exponential decay ($\beta=1$), whereas $C_{\mathrm{d}}(t)$ has a stretched exponential decay, with $\beta\approx 0.55$. It should be noted that the same stretched exponential behavior emerges in ${{K}}_{\mathrm{MCT}}(t)$. The possibility exists that, even though our numerical procedure appears to be converged and passes non-trivial consistency tests, that it is still not capable of describing accurately the tail of ${K}(t)$. Future work will be devoted to this important issue.

The autocorrelation times $\tau$ of the memory functions are shown in Fig. 6. Even though the autocorrelation times for $C_{\mathrm{d}}(t)$ are much larger at the same temperature than those of ${K}(t)$, both exhibit a form consistent with a power law growth $\tau\propto(T-T_{\mathrm{fit}})^{-\eta_{\mathrm{fit}}}$ which extrapolates to similar dynamical crossover temperatures $T_{\mathrm{fit}}$, but with different exponents $\eta_{\mathrm{fit}}$. We also include the autocorrelation times extracted from long-time fits of ${{K}}_{\mathrm{MCT}}(t)$. When fitting in the range $T\leq 0.7$, extrapolated dynamical temperatures derived from the exact memory function, the diagonal force autocorrelation and the MCT memory function are, respectively, $T_{\mathrm{fit}}^{\mathrm{(mem)}}=0.434(17)$, $T_{\mathrm{fit}}^{\mathrm{(diag)}}=0.412(3)$ and $T_{\mathrm{fit}}^{\mathrm{(MCT)}}=0.433(8)$. The exponents of the power-law growth are $\eta_{\mathrm{fit}}^{\mathrm{(mem)}}=1.1(3)$, $\eta_{\mathrm{fit}}^{\mathrm{(diag)}}=2.61(9)$ and $\eta_{\mathrm{fit}}^{\mathrm{(MCT)}}=1.7(3)$. Clearly the exact and approximate memory functions qualitatively exhibit similar forms of decay. However, the diagonal and MCT approximations produce a much more rapid growth of relaxation times.

A further measure of comparison between exact and approximate memory functions is obtained via the friction coefficient

[TABLE]

which we show in Fig. 7.

Power law fits, for $T\leq 0.7$, of the form $\zeta(T)=B_{\zeta}(T-T_{\zeta})^{-\eta_{\zeta}}$ show a similar scenario to the one found bt fitting the long-time behavior of ${K}(t)$ and $C_{\mathrm{d}}(t)$, with compatible $T_{\zeta}$ but different $\eta_{\zeta}$. The extrapolated dynamical temperatures are $T_{\zeta}^{\mathrm{(mem)}}=0.43(2)$ for the memory function, and $T_{\zeta}^{\mathrm{(diag)}}=0.414(6)$ for the diagonal function, and the exponents are $\eta_{\zeta}^{\mathrm{noise}}=1.4(3)$ and $\eta_{\zeta}^{\mathrm{diag}}=2.6(2)$. To rule out the possibility that our extrapolations are biased by preasymptotic effects arising from the nature of the short-time behavior of the memory functions which barely depends on temperature, we have recalculated the friction coefficients by removing from ${K}(t)$ and $C_{\mathrm{d}}(t)$ the short-time part of the decays fitted via Eq. (19). After this procedure, the qualitative situation remains unchanged, with the dynamical temperatures being $T_{\zeta}^{\mathrm{mem,long}}=0.42(2)$ and $T_{\zeta}^{\mathrm{diag,long}}=0.414(6)$, and the exponents $\eta_{\zeta}^{\mathrm{noise,long}}=1.6(4)$ and $\eta_{\zeta}^{\mathrm{diag,long}}=2.6(2)$. This can be understood by remarking that the plateaus of the normalized memory functions in Fig. 3 are at $\sim\mathcal{O}(1)$ height, whereas the span of the short-time regime is a tiny fraction of the total time.

VI Discussion

In this paper we have devised a numerically exact scheme to reconstruct the memory function for tagged particle motion in SCLs. The approach relies on the combination of the solution of a Volterra integral equation with numerically exact input and the numerical inversion of the solution for the memory function in Laplace space using the Gaver-Stehfest approach. The former method is accurate at short times while the latter is accurate at long times. The success of the concatenation of these approaches, which agree on intermediate time scales, is indicated via the reproduction of simulated transport data. Our approach thus opens the door for the detailed assessment of theories of SCLs that utilize the memory function approach.

We have compared the temperature dependence of the memory function extracted from molecular dynamics simulations to two distinct mean-field approaches. In particular, we consider the standard MCT tagged particle memory function filtered with exact structural and dynamical input as well as the unprojected diagonal force-force correlation function as approximate forms for the tagged particle memory function. The latter quantity is related to the exact memory function for a fluid interacting with short ranged forces in infinite dimensions. We find that the approximate approaches share qualitative (and in some cases quantitative) features with the exact memory kernel in three dimensions while also exhibiting important differences. Focusing on the diagonal force-force correlation function, we find that the short time behavior of the memory function is quantitatively reproduced, as is the extrapolated temperature of a putative power-law singularity of the long time relaxation. On the other hand the exponent of the growth of relaxation times as well as the form of the long-time relaxation of these two functions (stretched exponential versus exponential) differ markedly. This distinction highlights the importance of carrying out a fully self-consistent calculation of the mean-field dynamics of Ref. Maimbourg et al. (2016) as opposed to using exact non mean-field trajectories within the framework of the mean-field function $C_{\mathrm{d}}(t)$. The differing rate of growth of the relaxation could be attributed to this difference in the calculation procedure. On the other hand, some discrepancies, such as those associated with the stretching exponent $\beta$, might arise from the nature of our numerical determination of ${K}(t)$.

Our observations provide useful information concerning the ability of approximate mean-field approaches to capture the realistic behavior of particle motion in three-dimensional SCLs, as well as guidelines on the construction of more sophisticated theories that build on mean-field foundations. It should be noted that the comparisons made here do not employ fully self-consistent solutions of the mean-field dynamics, and instead use the exact three-dimensional trajectories in conjunction with mean-field expressions. It is possible that this inconsistency degrades the level of agreement with the exact memory kernel. Future work will be devoted to this issue. In addition, we will compare the exact memory function as presented here with advanced approaches that numerically bridge the behavior of the infinite dimensional and three dimensional limits such as “cluster dynamical mean-field” theories Biroli et al. (2019) based on Ref. Szamel (2017).

Acknowledgements.

We thank G. Szamel, P. Charbonneau, I. Dunn, B. Kloss, F.P. Landes, A. Manacorda and F. Zamponi for inspiring discussions. We thank Ian Dunn for providing his Prony function fitting package, and F.P. Landes for useful interactions about molecular dynamics simulations. This work was funded by the Simons Foundation for the collaboration “Cracking the Glass Problem” (No. 454951 to D.R. Reichman). This work benefited from access to the University of Oregon high performance computer, Talapas.

Appendix A Computing the memory function through the Laplace method

The Laplace transform of a function $C(t)$ is

[TABLE]

where $s\in\mathbb{C}$ is a complex number. Eq. (23) can be used to pass from Eq. (10) to Eq. (18) without any complication, because $C^{P}(t)$ is locally integrable.

Eq. (18) needs then to be inverted. The inverse Laplace transform, is defined by the Bromwich integral Josso and Larsen (2012)

[TABLE]

where $\gamma$ is a real number such that it is larger than the real part of all singularities.

To calculate the inverse transform we use the Gaver-Stehfest algorithm Gaver Jr (1966); Stehfest (1970); Spendier (2010), which has the additional benefit of giving a simple understanding of the errors and origins of possible instabilities Josso and Larsen (2012); Kuznetsov (2013).

The Gaver-Stehfest inversion algorithm approximates $C(t)$ through a finite series of functions

[TABLE]

where the coefficients are given

[TABLE]

with $\lfloor x\rfloor$ denoting the integer part of $x$. The values of $\omega_{k}$ can be calculated first, and stored in a hash table. In practice, the Laplace transform only needs to be evaluated for real values of $s$ 666To evaluate the functions at points off the grid of measurements, we fitted the data through a sum of exponentials using Beylkin and Monzón’s approximate Prony method Beylkin and Monzón (2005). The implementation we adopted was written for Ref. Dunn et al. (2019), and is publicly available at github.com/iansdunn/prony. , which we call $p\in\mathbb{R}$.

Noise and precision

The Gaver-Stehfest algorithm depends on a sum of $2M$ terms. The larger $M$, the more accurate the inversion. The drawback is that large-$k$ terms in this sum amplify small fluctuations, so $M$ cannot be too large because otherwise the noise is amplified, resulting in a huge loss in precision. Instructive benchmarks on the optimal $M$ are performed in Ref. Josso and Larsen (2012). A general rule of thumb is that $2M$ should correspond to the decimal digit where the errors appear, so if the only source of noise is the rounding error, then $M=7$ is the best choice in computations with double precision accuracy.

In our case, there are many sources of noise and systematic error. Even for exact functions, the fact that one interpolates between points on a grid is, per se, a source of error. In addition, our data is intrinsically noisy due to the fact that it arises from numerical simulations that must be averaged. Lastly, given the disparate time scales involved in the decay of ${K}(t)$, data needed to be collected on a (quasi) exponentially spaced time grid, implying that at large $t$ a single measurement is representative of a very large time span.

In this work we used $M=5$, which allowed for stable inverse transforms with our data at every temperature (although we needed a large number of trajectories, see Tab. 1). A low value of $M$ implies a loss in accuracy on fine-grained scales, which we show has no dramatic consequence on the reconstruction of the memory kernel (see Appendix D).

Integration ranges

Another issue is with the Gaver-Stehfest method is that it requires the evaluation of the Laplace transforms at very large values of $p$, with

[TABLE]

where $p_{\mathrm{max}}$ is the largest value of $p$, and $t_{\mathrm{min}}=dt=0.0025$ is the smallest time. When $p$ is too large, the integrand in (23) incurs numerical underflow. This implies again that $M$ cannot be too large, and that time cannot be too small. As a consequence, the Laplace inversion becomes extremely noisy at small $t$, which is why concatenation with a method suitable at short times is required.

Appendix B Computing the memory function through the Volterra integration

Here, we follow Ref. Press et al. (2007), with minor adaptations to our problem. Eq. (12) is a Volterra equation of the second kind. By defining $t_{j}=jdt~{}(j\in\mathbb{N}^{+})$, and $A_{j}=A(t_{j})$, we can rewrite it using a trapezoidal integration scheme on a linear grid

[TABLE]

Since $C^{FP}(t,u)=C^{FP}(t-u)$, then $C^{FP}_{ii}=0$, thus Eq. (28) becomes

[TABLE]

Since $C^{FP}(t)$ decays to zero, at large times, its signal-to-noise ratio is low. In order to improve accuracy, we can define a $t^{*}=i^{*}dt$ such that $\mathcal{C}(t)\approx 0~{}\forall t>t^{*}$, and $\mathcal{C}(t-u)\approx 0~{}\forall u<t-t^{*}$. Similar procedures are used to reduce the statistical error in spin-glass correlation functions Belletti et al. (2009), with $t^{*}$ the time when the signal-to-noise ratio becomes smaller than an order 1 number. Since the grid is linear, $\mathcal{C}_{ij}\approx 0~{}\forall j<i-i^{*}$, so Eq. (29) further simplifies to

[TABLE]

Generic grid

If the time grid is non-linear, as in our case, the discretized Volterra equation is

[TABLE]

where now $t_{i}$ is the $i^{\mathrm{th}}$ element on the ordered time grid. As in the case of the linear grid, the summation can be truncated in order to avoid time regions with low signal. In this case, the summation in Eq. (31) starts at $j=\max(1,i-i^{*})$, where $i^{*}$ is $i:t^{*}=t_{i^{*}}$.

In principle the same procedure may be used with other integration schemes, such as the Simpson scheme, though instabilities may arise depending on how the integration boundaries are treated Press et al. (2007). As we show in Appendix C, the Volterra method with a trapezoidal integration scheme is very stable and accurate for small $t$.

Appendix C Consistency checks

Our method for extracting the memory function consists of the piecewise concatenation of the Volterra method with the Laplace method. In Fig. 8 we show this concatenation at $T=0.47$, with jackknife error bars. Despite the small fluctuations in the Volterra curves at all temperatures, their behavior becomes unphysical at the end of the short-time regime, and consequently the memory function calculated with the Volterra method fails the test shown at the end of this section.

The Laplace method gives the correct behavior at long times, but suffers from dramatic fluctuations at short times for reasons described in Appendix A. It is to be noted that at lower temperatures the error bars become large with respect to the signal. However, due to systematic sources of error mentioned in Appendix A, the fluctuations are non-Gaussian, and are thus not a good indicator of the reliability of the calculations.

In order to assess the reliability of the computed memory function, we verify that it satisfies Eq. (10). In Fig. 9 we show the result of this consistency check for the integrated version of Eq. (10),

[TABLE]

since in our view $C^{P}(t)$ can be more readily interpreted by eye than $C^{FP}(t)$. We see that the entire non-trivial part of $C^{P}(t)$ is well-reproduced. At longer times, the $C^{P}(t)$ calculated from the reconstructed memory functions fluctuate wildly, though the time at which the wild fluctuations begin grows as the number of trajectories used for averaging is increased. clearly, however, for the number of trajectories used here we can achieve full consistency for time scales where $C^{P}(t)$ has nearly entirely decayed to zero.

Appendix D Increasing $\boldsymbol{M}$

As stated in Appendix A, the number of coefficients, $2M$, in the inverse Laplace transform calculated with the Gaver-Stehfest method has an influence on the fine structure of a function that can be resolved. Increasing $M$ allows for the access of finer details in the behavior of ${K}(t)$, but doing so also increases numerical instability.

As shown in Fig. 10, varying $M$ from 3 to 5 produces a variation in the plateau of ${K}(t)$, but the long-time behavior remains unchanged. One can expect, therefore, that increasing to $M=7$ (the maximum possible $M$ for double precision, see Appendix A) will not lead to significant changes. Note, moreover, that the plateau heights shown in Fig. 5b are evaluated from the long-time part of ${K}(t)$ which, as shown in Fig. 10, does not vary with $M$. Consequently, even though the qualitative shape of the plateau changes, we do not expect our measurements of the plateau height to change if $M$ is increased.

In Fig. 11 we show that the consistency check presented in Appendix C improves as $M$ grows, systematically up to the value of $M$ used in this work.