Configuration-tree Theoretical Calculation of the Mean-Squared Displacement of Particles in Glass Formers

TL;DR

This paper presents the first analytical calculation of the mean-squared displacement in glass formers using a configuration-tree theory, showing excellent agreement with simulations across various temperatures.

Contribution

It introduces a novel analytical approach based on a local random configuration-tree theory to evaluate particle MSD in glasses, validated by numerical simulations.

Findings

Good quantitative agreement with simulations across temperatures

First analytical calculation of MSD in glasses using this method

No free parameters used in the model

Abstract

We report an analytical evaluation of the mean-squared displacement (MSD) of the particles in glasses based on their coarse grained trajectories. The calculation is conducted by means of a local random configuration-tree theory that was recently proposed by one of us [C.-H. Lam, J. Stat. Mech. \textbf{2018}, 023301 (2018)]. Results are compared with the numerical simulations of a lattice glass model, and good quantitative agreement has been obtained over a wide range of temperatures in the entire region of time with virtually no free parameters. To the best of our knowledge, the calculation is the first in its kind.

Equations

$$\text{mbox}MSD(t)={\mathcal{N}}^{-1}\sum _j{\mathbf{d}}_j^2(t),$$

$${\tilde{w}}_0={\int }_{-\infty }^{\infty }d(\Delta E)w(\Delta E)P(\Delta E),$$

$${\tilde{w}}^{-1}=\frac{{\int }_{-\infty }^{\mathcal{C}{k}_{B}T}d(\Delta E)w(\Delta E{)}^{-1}P(\Delta E)}{{\int }_{-\infty }^{\mathcal{C}{k}_{B}T}d(\Delta E)P(\Delta E)},$$

$$N{p}_0(t+dt;m)-N{p}_0(t;m)=-N{p}_0(t;m)dtmz{\tilde{w}}_0+N{p}_1(t;m)dt\tilde{w},$$

$${\partial }_t{p}_0(t;m)=-mz{\tilde{w}}_0{p}_0(t;m)+\tilde{w}{p}_1(t;m).$$

$${\partial }_t{p}_1(t;m)$$

$${\partial }_t{p}_{n>1}(t;m)$$

$$\bar{n}(t;m)=\sum _{n}n{p}_n(t;m),$$

$${\partial }_t\bar{n}(t;m)=(c_m-1)\tilde{w}\left[1-p_0(t;m)\right]+mz{\tilde{w}}_0{p}_0(t;m),$$

$$\bar{n}(t;m)={\int }_0^{t}d{t}^{\prime }\{(c_m-1)\tilde{w}\left[1-p_0(t^{\prime };m)\right]+mz{\tilde{w}}_0{p}_0(t^{\prime };m)\}.$$

$$\bar{n}(t;m)\approx (c_m-1)\tilde{w}t$$

$$P_n(s;m)={\int }_0^{\infty }dt{e}^{-st}{p}_n(t;m),$$

$$s{P}_0(s;m)-1=-mz{\tilde{w}}_0{P}_0(s;m)+\tilde{w}{P}_1(s;m),$$

$$s{P}_1(s;m)=-(1+c_m)\tilde{w}{P}_1(s;m)+mz{\tilde{w}}_0{P}_0(s;m)+\tilde{w}{P}_2(s;m).$$

$$P_n(s;m)={\alpha }_n(s;m)+{\beta }_n(s;m)P_2(s;m),$$

$$\displaystyle\left[\begin{tabular}[]{ccc}$\alpha_{0}(s;m)$&$\beta_{0}(s;m)$\\ $\alpha_{1}(s;m)$&$\beta_{1}(s;m)$\end{tabular}\right]$$

$$s{P}_n(s;m)=c_m\tilde{w}\left[P_{n-1}(s;m)-P_n(s;m)\right]+\tilde{w}\left[P_{n+1}(s;m)-P_n(s;m)\right].$$

$$P_n(s;m)=P_2(s;m)c^{n-2},n>1.$$

$$c^2-(1+c_m+s{\tilde{w}}^{-1})c+c_m=0.$$

$$c(s;m)=\frac{1}{2}\left(1+c_m+s{\tilde{w}}^{-1}-\sqrt{(1+c_m+s{\tilde{w}}^{-1}{)}^2-4c_m}\right).$$

$$\sum _{n=0}^{\infty }{P}_n(s;m)={\int }_0^{\infty }dt{e}^{-st}=s^{-1}.$$

$$P_2(s;m)=\frac{s^{-1}-{\alpha }_0(s;m)-{\alpha }_1(s;m)}{{\left[1-c(s;m)\right]}^{-1}+{\beta }_0(s;m)+{\beta }_1(s;m)},$$

$$p_n(t;m)=\frac{1}{2\pi i}{\int }_{\gamma -i\infty }^{\gamma +i\infty }ds{e}^{st}{P}_n(s;m),$$

$$\text{mbox}MSD(t;m)=\frac{\tilde{l}{a}^2}{\mathcal{V}\varphi }\bar{n}(t;m).$$

$$f(m;\bar{m})=\frac{{\bar{m}}^m{e}^{-\bar{m}}}{m!}.$$

$$\text{mbox}MSD(t)=\sum _{m}f(m;\bar{m})\text{mbox}MSD(t;m)=\frac{\tilde{l}{a}^2}{\mathcal{V}({\Omega }^{-1}-{\varphi }_v)}\sum _{m}f(m;\bar{m})\bar{n}(t;m).$$

$$D=\underset{t\to \infty }{\lim }\frac{\text{mbox}MSD(t)}{2dt},$$

$$D=\frac{\tilde{w}\tilde{l}{a}^2\Omega }{2d\mathcal{V}(1-{\varphi }_v\Omega )}{\sum _m}^{\prime }f(m;\bar{m})(c_m-1),$$

$$w(\Delta E)=w_0{e}^{-(E_0+\Delta E/2)/k_{B}T},$$

Full text

1 Introduction

Upon rapid cooling, many liquids undergo a dramatic – often more than ten decades – increase in their viscosity and eventually become glasses at a temperature $T_{g}$. The question of how the constituent particles of a liquid lose their mobility in this process remains a source of dispute, having generated a library of theories of varying sophistication and complexity [1, 2, 3].

Glasses are notoriously complex systems [4]. Analytical yet accurate evaluations of any of their properties are very rare and valuable. The mean-squared displacement (MSD) of particles, defined as

[TABLE]

where $\mathbf{d}_{j}(t)$ is the displacement of the $j$-th particle during a time $t$ and $\mathcal{N}$ is the total number of particles in a glass, is often computed to study glassy transition. As a ubiquitous signature of glass formation, upon approaching $T_{g}$, MSD$(t)$ develops a plateau, which is a precursor of the so-called $\alpha$ relaxation characterizing any glasses [5].

In this contribution, we present an analytical calculation of MSD$(t)$ on the basis of a theory that was recently proposed by one of us [6], belonging to the class of dynamic facilitation descriptions [7, 8, 9, 2, 10, 11, 12]. Going beyond related calculations in Refs. [13, 14], our calculation covers the full range of time scales, spanning both the time window before and that after the plateau and hence providing a complete picture of the glassy dynamics. The theoretical result is then benchmarked against the ‘distinguishable particle lattice model’ (DPLM) [15], which is a lattice glass model realizing a dynamic facilitation mechanism motivated by glassy film simulation observations [16]. Satisfactory agreement has been achieved – including the position of the MSD plateau – for various particle densities and over a wide range of temperatures by adjusting only two parameters of the model, which can further be fixed by the diffusion coefficient $D$ of the particles resulting in no additional free parameter.

2 Theory

2.1 The concept of a random local configuration tree

Let us start with a brief recapitulation of the basic concepts of the theory that has already been described in details in Ref. [6]. For the completeness of this paper, however, we have also provided in the Appendix a more comprehensive discourse of the main points and a comparison with the common approaches in the literature.

At temperatures close to $T_{g}$, a particle in a liquid mostly vibrates about a temporary local equilibrium center, and occasionally makes a transient hop from one center to another. In the coarse-grained picture with the fast vibrations averaged out, only such hops are interesting and the phase space of the system reduces to the space of particle configurations. The system would then transit from one configuration to another via one or a sequence of particle hops – which are assisted by the vibrations acting as an effective heat bath – and relax its structure. The hopping distance is typically comparable to the particle size as measured in many works (see Ref. [17] and references therein) from the van Hove correlation function exhibiting a peak. Each hop is primarily assumed to be initiated by a void. The void propagates along a string of particles that simultaneously displace one another during the propagation. These motions, called micro-strings, generally come with different lengths obeying a distribution narrower than an exponential distribution [18]. As an approximation, we then assume that the micro-strings are of the same length and each possesses $\tilde{l}$ particles. In general, a void can propagate along one of $z$ strings lying in different directions. However, these strings each sense their own local energy landscape and they are not equally favorable.

To simplify the analysis, let us divide the liquid into a number of local regions, each of volume $\mathcal{V}$. We assume that $\mathcal{V}$ is big enough so that a void propagation in one region does not affect the structural relaxation of neighboring regions, yet small enough so that one void propagation energetically affects the propagation of other voids in the same region (i.e. the fully interacting approximation [6]). Under this approximation, it suffices to consider the dynamics of one local region. Let us represent the configurations of this region by the nodes of a graph and possible transitions by the edges. Any pair of nodes connected by one edge signify the presence of a micro-string. In a region containing $m$ voids, a node has $mz$ edges linked to it. The fully interacting approximation allows us to ignore the correlations between the edges. The region can revisit a previous node by retracing of a path of adjoining edges in reversed direction. However, we neglect the possibility of a loop of edges, so that the graph becomes a tree with its root corresponding to the initial configuration. The validity of this approximation stems from the very high dimensions of the particle configuration space and this has been discussed in details in Ref. [6]. It should be further noted that, even with the inclusion of some loops, one may still represent the graph as a tree. This amounts to some additional correlations between edges, which should not alter our results significantly.

It should be seen that the configuration tree is random in the sense that, in general its energetically favorable edges are randomly distributed. In this paper, we assume that the transition probabilities for the energetically favorable edges are finite and fixed but negligible for the unfavorable edges. In Figure 1, an illustration of the energy landscape is shown for the local configuration tree with one void and $z=4$.

We should also mention that, the tree theory sketched above does not consider vibrations explicitly and more precisely captures the dynamics of the so-called inherent structures [19]. An inherent structure is nothing but a microscopic configuration of meta-stable vibration centers. A void propagation then takes the system from one inherent structure to another. As to be discussed later, our theory is consistent with the numerical simulations carried out in Ref. [20] of the evolution of the inherent structures of a binary Lennard-Jones liquid.

In the rest of this section, we derive and solve a set of rate equations governing the motions of a local region as described by the tree.

2.2 Level dynamics on the tree: rate equations

Let us consider a local region with $m$ voids and suppose the region is at equilibrium initially sitting at the root of its configuration tree defined to be level $n=0$. In time $t$, the region will descend down the tree to some level $n$ with a probability denoted by $p_{n}(t;m)$. Our purpose here is to derive the rate equations satisfied by this probability. Our formulation is analogous to previous calculations for random walks on regular trees [21, 22]. It is emphasized again that, in our theory the random walk takes place primarily in the configuration space rather than in the real space.

To ease the discussion, let us look at an ensemble of $N$ identical regions all initially at the root. At instant $t$, the number of members in this ensemble at the root ($n=0$) is $Np_{0}(t;m)$. After a short interval $dt$, this number becomes $Np_{0}(t+dt;m)$, the difference due to a loss of those members transiting to level $n=1$ and a gain of members at level $n=1$ transiting to the root. The transition from the root to level $n=1$ can take place through $mz$ edges [6], which we label by $\nu=1,2,...,mz$. Let us denote by $dt\times\tilde{w}^{\nu}_{0\rightarrow 1}$ the probability of such a transition via edge $\nu$. The total transition rate is then $\sum_{\nu}\tilde{w}^{\nu}_{0\rightarrow 1}$ which averages to $mz\tilde{w}_{0}$, with $\tilde{w}_{0}$ being the rate per edge. As the transitions are assisted by an effective heat bath of vibrations, $\tilde{w}^{\nu}_{0\rightarrow 1}$ in general depends on the energy difference $\Delta E$ between the configurations connected by the $\nu$-th edge. Denoting by $P(\Delta E)$ the probability for the occurrence of $\Delta E$ assuming an initial equilibrium configuration, a formally exact expression for the average transition rate $\tilde{w}_{0}$ along any of the edges from the root is

[TABLE]

where the rate $w(\Delta E)$ depends on the detailed microscopic dynamics. Now the number of members lost to level $n=1$ is $Np_{0}(t;m)dt~{}mz\tilde{w}_{0}$.

On the other hand, the transition for any member from level $n=1$ to the root can take place only through the edge linking the root and the node where the member is currently located. The probability of this transition for a member sitting on the node of the $\nu$-th edge is then $dt\times\tilde{w}^{\nu}_{1\rightarrow 0}$. Denote by $N^{\nu}_{1}$ the number of such members, which satisfies $\sum_{\nu}N^{\nu}_{1}=Np_{1}(t;m)$. The number of members gained to the root is then $\sum_{\nu}N^{\nu}_{1}dt~{}\tilde{w}^{\nu}_{1\rightarrow 0}$. At low temperatures, very few edges are energetically favorable and $N^{\nu}_{1}$ are significant only for the nodes linked to these edges. The rest can then be ignored. As $\Delta E$ for energetically favorable edges linking only low energy configurations are expected to be small, we can neglect the difference in their transition rates. As such, we shall replace $\tilde{w}^{\nu}_{1\rightarrow 0}$ with an average rate $\tilde{w}$ for them. To actually compute $\tilde{w}$, let us define the energetically favorable edges as those for which $\Delta E\leq\mathcal{C}k_{B}T$, where $\mathcal{C}$ is a constant of order unity [6]. Now, the rate $\tilde{w}$ for a favorable edge is approximated as follows

[TABLE]

which will be further explained below. As such, the number of members gained amounts to $Np_{1}(t;m)dt~{}\tilde{w}$.

Combining these considerations, we obtain according to the detailed balance principle that

[TABLE]

or equivalently

[TABLE]

By applying the same reasoning to other levels, we find

[TABLE]

Here we have introduced the quantity $c_{m}=(mz-1)q$, with $q=\int^{\mathcal{C}k_{B}T}_{-\infty}d(\Delta E)P(\Delta E)$ giving the likelihood of an edge being energetically favorable. $c_{m}$ gives the number of children of a node at level 1 or beyond. One can readily show that $\partial_{t}\sum_{n}p_{n}(t;m)=0$, i.e. the total probability ($=1$) is conserved. The equations imply slightly different properties of the root node, which is however an artifact that arises because of all the transition rates only that (i.e. $\tilde{w}_{0}$) from level 0 to 1 can be computed exactly while the rest ($\tilde{w}$) are approximations.

Equation (3) provides $\tilde{w}$ by averaging the transition time $w(\Delta E)^{-1}$. At low temperatures, the long-time dynamics is dominated by a series of transitions along the levels on a thin tree branch with a degree $c_{m}$ close to unity for some dominant value of $m$. The average of the transition time, as opposed to that of the transition rate, thus gives a more accurate description for long-time dynamics. Furthermore, Eq. (3) is accurate even at higher temperatures. This is because $w(\Delta E)$ then admits small fluctuations and averaging over time or rate indeed give similar values.

From Eqs. (4) - (6), it follows that the average level

[TABLE]

the local region could have reached must satisfy the following equation,

[TABLE]

subjected to initial conditions $\bar{n}(0;m)=0$ and $p_{0}(0;m)=1$. An integration gives

[TABLE]

It is interesting to see that $\bar{n}(t;m)$ is completely determined if $p_{0}(t;m)$ is known. For small $t$, we may take $p_{0}(t;m)\rightarrow 1$, and then $\bar{n}(t;m)\rightarrow mz\tilde{w}_{0}t$, i.e. the initial growth of $\bar{n}(t;m)$ which is exact.

Inspecting Eq. (8) or (9), we note that there are two types of solutions depending on whether $c_{m}$ is greater or less than $1$. For $c_{m}>1$, $\bar{n}$ increases without bound, implying that the local region will be descending down the tree and access all of its configurations in time. In the meantime, $p_{0}(t;m)$ must decay to zero, yielding

[TABLE]

for large $t$. We may call such regions mobile regions.

For $c_{m}<1$, however, $\bar{n}$ cannot increase arbitrarily. Instead, it will initially increase as the second term on the right-hand side of Eq. (8) is larger than the first term, but finally saturate when the two terms balance each other. In this regime, the local region is trapped in a few configurations, and cannot escape. These regions may be called immobile regions.

2.3 Solutions to the rate equations

While Eqs. (4) - (6) may be numerically integrated in a straightforward manner, here we solve them by the method of Laplace transform. So we define

[TABLE]

where $s$ resides in the upper half complex plane. In terms of $P_{n}(s;m)$, Eqs. (4) and (5) can be written as

[TABLE]

Here we have used the initial condition that $p_{n}(0;m)=\delta_{n,0}$. Expressing $P_{n=1,0}$ in terms of $P_{2}$, we find

[TABLE]

where $\alpha_{n=0,1}$ and $\beta_{n=0,1}$ are some coefficients given as follows

[TABLE]

Analogously, Eq. (6) is transformed as

[TABLE]

This equation is invariant under a translation of $n$ and thus the general solution is of the form

[TABLE]

The quantity $c$ is obtained by substituting this expression into (21). We find

[TABLE]

This has two solutions. One solution leads to $P_{n}$ that diverges for large $n$ and must be excluded. The other solution is given by

[TABLE]

Finally, the normalization condition of the probabilities implies that

[TABLE]

Combining this with Eqs. (14) and (22), we obtain

[TABLE]

which can be used to generate all $P_{n}(s;m)$. Going back to $p_{n}(t;m)$, we need to invert the Laplace transform

[TABLE]

which can be numerically evaluated using standard algorithms [23, 24].

3 Calculating MSD$(t)$

The migration of a region on the configuration tree can be mapped onto the real-space motions of the particles in this region, thereby allowing us to calculate the MSD, MSD$(t)$ of these particles.

To this end, let us observe that, as the region moves down the tree by one level, what happens in real space is that a void propagates along a string and each particle in this string hops by an average distance $a$ comparable to the atomic separation. Suppose the region has descended to level $n$ in a period $t$. During this period, in total $n\tilde{l}$ particle hops have occurred in $n$ strings each of length $\tilde{l}$, excluding the hops due to back-and-forth repetitions between any pair of nodes, which are highly correlated and do not contribute to the MSD. Some particles may have hopped more than once. As long as the hops are uncorrelated, as is statistically reasonable for non-repetitive hops, the squared displacement of a particle over ensemble average is simply $a^{2}$ times the number $n_{hop,j}$ of non-repetitive hops particle $j$ has executed, i.e. $\langle\mathbf{d}^{2}_{j}\rangle=n_{hop,j}a^{2}$ with $\sum_{j}n_{hop,j}=n\tilde{l}$. The contribution of each non-repetitive hop to the MSD is just $a^{2}/\mathcal{V}\phi$, where $\phi$ is the particle density and $\mathcal{V}\phi$ is the number of particles in the region. As such, for a region reaching level $n$, the MSD is $n\tilde{l}a^{2}/\mathcal{V}\phi$. Ensemble averaging this by $p_{n}(t;m)$, we obtain the MSD for this region as

[TABLE]

Here $\bar{n}(t;m)$ can be directly obtained from Eq. (9) with $p_{0}(t;m)$ provided in Sec. 2.3.

To obtain the MSD for the entire liquid, we have to take an average of MSD$(t;m)$ over all local regions. The number of voids in one region may differ from that in another. Suppose the average number of voids per region is $\bar{m}=\mathcal{V}\phi_{v}$, where $\phi_{v}$ is the density of voids. With small correlations between voids, their steady-state distribution over the regions then approximately follows Poisson’s law, i.e. the probability of a region with $m$ voids is

[TABLE]

Now averaging MSD$(t;m)$ by this distribution, we finally arrive at

[TABLE]

Here $\Omega$ is the volume of a particle or of a void, i.e. $\Omega^{-1}=\phi+\phi_{v}$. Note that Eq. (30) is closed with Eq. (9), which gives $\bar{n}(t;m)$.

From Eq. (30), we can obtain the diffusion coefficient $D$ for the liquid. By definition,

[TABLE]

where $d$ is the dimension of the system. As discussed in Sec. 2.2, for mobile regions with $c_{m}=(mz-1)q>1$, $\bar{n}(t;m)=(c_{m}-1)\tilde{w}t$ for $t\rightarrow\infty$. For immobile regions, however, $\bar{n}(t;m)$ approaches a constant for $t\rightarrow\infty$. Thus, only mobile regions contribute to $D$. We find

[TABLE]

where the prime restricts the summation to mobile regions. This result on $D$ was first obtained in Ref. [6], where the behavior of $D$ has been thoroughly investigated.

4 Comparison with simulations

While the theory is expected to be valid for a variety of glassy systems exhibiting string-like motions, here we check it with the DPLM, for which most of the required parameters are known precisely and adjustable parameters are few. As it is a lattice model without vibrations, our theory can be directly applied. For comparison with MD simulations or experiments in the presence of vibrations, as aforementioned, the MSD must be based on highly coarse-grained trajectories with vibrations completely eliminated, such as those of inherent structures [20, 25, 26, 27]. While our focus is on the DPLM, we wish to point out that simulations have been done in Ref. [20] by Liao et al on the inherent structures of a binary Lennard-Jones liquid and their results are in good qualitative agreement with our theory. A very important result of their work is that, the MSD of inherent structures initially increases linearly (rather than parabollically) in time and then goes through a shallow plateau before finally entering the diffusion regime. These features are also borne out in our theory.

As the DPLM has already been studied thoroughly in Ref. [15], here a short description is provided for the convenience of readers. The DPLM assumes a number of distinguishable hard-core particles living on a square lattice with lattice constant $a=1$. Interactions occur between neighboring particles in the form $V_{ijs_{i}s_{j}}$, each of which is sampled from a distribution $g(V)$ before a simulation begins. Here $i$ and $j$ denote a pair of neighboring sites while $s_{i}$ and $s_{j}$ denote the types of the particles sitting at the sites, respectively. Each particle is of its own type. A void is simply the absence of a particle. We have chosen $g(V)$ to be a uniform distribution in the range of $V\in[-0.5,0.5]$. The string length is chosen to be $\tilde{l}=1$. The kinetic Monte Carlo rules of activated hopping is employed for the dynamics of the particles. Each particle can hop to an empty neighboring site at a rate

[TABLE]

where $w_{0}$ and $E_{0}$ are constants and $\Delta E$ is the change in the system energy due to the hop. Though the equilibrium properties of this model are exactly solvable, it displays many characteristic behaviors, including the two-step relaxation of a glass former [15].

Note that $\Delta E$ is also the energy change due to a void propagation. In the configuration-tree theory, we call an edge in the local tree energetically favorable if $\Delta E$ for the corresponding void propagation is less than an upper cutoff $\mathcal{C}k_{B}T$, where $\mathcal{C}$ is a parameter of the order of unity. Thus, the probability $q$ that an edge is energetically favorable can be directly evaluated from $g(V)$ and the equilibrium properties of the DPLM, see Ref. [15]. Since $\Omega=1$ for the DPLM, this parameter $\mathcal{C}$ and the volume of a local region $\mathcal{V}$ are the only parameters needed to evaluate MSD$(t)$ by Eq. (30). These parameters can actually be further fixed by fitting the analytical expression of $D$, Eq. (32) to the numerical simulations of the DPLM. This has been carried out in Ref. [6] using $\tilde{w}=w_{0}$, where it was found that $\mathcal{C}=1.7$ and $\mathcal{V}=12$ give the best fit over a wide range of temperatures and void densities. This choice is appropriate for both the activated hopping dynamics in (33) and the Metropolis dynamics, with a wide choice of $g(V)$, where $\Delta E$ averages to 0. Details on the Metropolis dynamics will be given elsewhere. As such, virtually no free parameters are involved in the analytical calculation of MSD(t).

Using the Laplace-transformation approach, by applying Eq. (26) together with the coefficients from Eqs. (20) and (24) as well as the rates defined by Eq. (2) and Eq. (3), one can calculate $P_{2}(s;m)$. Then, $P_{0}(s;m)$ follows from Eqs. (12) and (13) and can be inserted into Eq. (27) to obtain $p_{0}(t;m)$. Substituting $p_{0}(t;m)$ into Eq. (9), we find $\bar{n}(t;m)$, which is further inserted in Eq. (30) to obtain the MSD. Eq. (30) and these associated equations thus constitute our closed form solution for the MSD. Alternatively, identical results can be obtained using an ordinary differential equation (ODE) integration approach. In this case, the rate equations in Eqs. (4) - (6) are numerically integrated to obtain, in particular, $p_{0}(t;m)$. The MSD can then be obtained similarly using Eqs. (9) and (30).

The agreement of the theory with the simulations for a wide range of temperatures is shown in Figure 2. In the plot, we have chosen a void density $\phi_{v}=0.01$ for the left panel while $0.05$ for the right panel. The temperature values range from the non-glassy liquid phase at $T=4$ to the deeply supercooled phase at $T=0.15$ [15]. Here, glassy behaviors at low $T$ manifest as a plateau, which is relatively shallow only because of the lack of vibrations in the model and by no means implying a slight or moderate supercooling. Indeed, this shallow plateau is also a property of the MSD observed in lattice models [28] and of inherent structures of a binary Lennard-Jones liquid [20]. As shown in the figure, the agreement is very good at high temperatures, while at low temperatures the theory deviates from the simulations but the agreement is still impressive. The nearly exact agreement for MSD$(t)\leq 2\times 10^{-3}$ for all curves also shows the validity of the rates $\tilde{w}_{0}$ used in our theory. More importantly, the theory captures the location of the plateau very well for the full wide range of temperatures studied. Only at very low temperatures $T\leq 0.2$, the predicted plateaus become too flat compared with the more shallow ones from simulations. The discrepancy may be due to neglecting random fluctuations in the theory as it adopts averaged transition rates and assumes local regions completely isolated from other parts of the system. At any finite $T$, our theory predicts no ideal phase transition [6]. Therefore, even at arbitrarily low but finite $T$, an increasingly broad plateau in the MSD is predicted, which eventually must cross over to the long-time diffusive regime.

5 Conclusion and outlook

Based on the local random configuration-tree theory, we have analytically calculated the MSD based on the inherent structures of glass formers. The results have been compared with numerical simulations of a lattice glass model – the DPLM, and an overall satisfactory agreement has been achieved at the quantitative level with few free parameters.

In this paper we have only applied the theory to the DPLM, for the required parameters are precisely known here. In the near future, we shall apply it to other glass formers, including computational (e.g. polydispersed Lennard-Jones liquids) and experimental systems. For such systems, the MSD should be based on coarse-grained trajectories [17] for inherent structures with fast vibrations averaged out. There will be far more parameters involved. However, these parameters fall in groups and each group can be tightly constrained by computational or experimental observations. We shall also use the theory to shed light on the effects of particle swap on particle dynamics in these systems, which has recently attracted lots of attention [29, 30, 31, 32, 33, 34, 35].

Another interesting line of investigation would be to look at other physical quantities from the configuration-tree theory point of view, such as the self-intermediate scattering function and the particle displacement distribution function. The former is especially important for its widespread use in the experimental literature. One should note that the MSD is closely related to this quantity in the long wavelength limit. We are currently working on generalizing our theory to calculate these functions.

Finally, we make a few remarks on the role of voids in our theory. There are two basic tenets in the theory. The first is the existence of a local random configuration tree, which in our opinion is quite general and robust. The second explicitly assumes the presence of voids and string-like particle motions. While string-like motions have ample evidences, the concept of voids is sometimes contested. Regarding this the following point is worthy of notice: A void capable of initiating string-like motions is not necessarily as well-defined spatially as in lattice gas models such as the DPLM. It could be a quasi-void consisting of a relatively small void coupled to free volumes fragmentary and distributed [17].

We thank Giorgio Parisi for helpful comments. We are grateful for the support of Hong Kong GRF (Grant 15330516) and Hong Kong PolyU (Grant 1-ZVGH).  

Appendix: Summary of main assumptions and comparison with the literature

Here we summarize the main concepts and implications of the theory and discuss their similarities and differences to other approaches in the literature.

The starting assumption of the tree theory is that elementary particle motions are activated particle hops initiated by voids, an idea that can be traced back to Glarum’s theory of diffusion of defects identifiable as units of free volume or voids [36] and activated motions introduced by Goldstein [37]. A hop takes the system from one inherent structure [19] to another. The hopping distance is assumed fixed and comparable to the particle diameters, as measured from the short-time van Hove correlation function [38, 39, 17]. Our theory is thus microscopic, based on elementary atomistic processes, unlike for example continuous-time random walk (CTRW) models [40, 41, 13], which assumes a broad distribution of hopping distances appropriate at a more coarse-grained level. Voids of volumes comparable to particle sizes may be too few to account for observed motions [42]. Noting the non-rigidity of fluids, in our theory we conceive that the a void may also be a quasi-void which consists of a relatively small interstitial void and fragmented free volumes close by, so that the squeezing-in of another particle can happen [17].

In addition, particle hops can be correlated in the form of microstrings [18, 10] each describing simultaneous coherent hops of a line of particles. The length of microstrings is short and narrowly distributed [18] and thus approximately fixed. This is in contrast with the length of incoherent strings, which has been phenomenologically assumed to diverge upon cooling[43]. 2. 2.

The second assumption is that the particles are assumed distinguishable with particle-dependent pair interactions. Particles with diverse properties have been modeled for example in Refs. [44, 45]. This is directly justifiable for polydispersive and polymer systems. It may also be applicable to identical particles as it effectively models different frustration states at various meta-stable particle positions [15]. In contrast, however, most studies including defect models usually assume particles with identical properties [36, 7, 28, 46, 43]. 3. 3.

As a consequence of (ii), the system energy change involved in a particle hop then follows from the system potential energy landscape (PEL), which is a function defined in the particle configuration space of distinguishable particles. Unlike spin-glass models, the PEL accounts for quenched disorder in the configuration space but not in the real space [47]. The configuration space contains all possible inherent structures of the system. This space is substantially larger than that of identical particles because the exchange of particles, for example, leads to a different configuration with in general a different energy. The PEL is a rugged energy landscape [37, 48]. Hops of particles and voids are thus history dependent, as further explained in point (vi), in contrast to history independent hops in many theories [36, 43, 14]. Motions of voids are then described by random walks in the configuration space, analogous to previous descriptions for particles in Refs. [48, 49]. 4. 4.

A further consequence of (ii) is that the system energy change depends on the paths of the voids and in particular on the time-order of the paths. We have described the phenomenon as ”path interactions of voids” or equivalently ”string interactions” [17]. This is because a void’s motions alter the particle configurations along its whole track and perturb subsequent motions of other voids. Therefore, the hop of a void changes the hopping rate of another void close by. In contrast, voids are assumed non-interacting usually [36]. 5. 5.

In accord with (iii), the particle configuration space for a local region with $m$ voids is approximated by the Bethe lattice with a coordination number proportional to $m$. Each edge represents a hop associated with an independent energy change. The configuration space is of a very high dimension and is much closer to a tree than the physical space. For a local region with a single void, the tree geometry is indeed exact up to the fifth tree level for our lattice model [15]. For a region with multiple voids, a fully interacting void approximation also leads to a tree geometry with independent edge properties [6]. Discrepancies result effectively at a reduced number of independent tree branches. Nevertheless, qualitative features of our theory are robust as long as such an effective tree coordination number increases with $m$. 6. 6.

As an approximation, the nodes in the Bethe lattice are classified as energetically favorable or unfavorable. The favorable nodes constitute a random tree that is assumed to have a flat energy landscape. Within the random tree which is a sub-tree of the Bethe lattice, our theory reduces to a purely kinetic one analogous to the kinetically constrained models (KCM) [2]. A key feature of our theory is that the energetics, i.e. the PEL with the disorder quenched in the configuration space, is accounted for by the geometry of the random tree. The quenched disorder and the resulting history dependence of voids motions are then treated as the confinement to a fixed sub-tree in the configuration space. It is closely analogous to, for example, studying diffusion in fractal channels constituting quenched disorder in the physical space in porous media [50]. 7. 7.

**Particle mobility in a local region with $m$ voids reduces to an $m$-dependent percolation problem of the random tree in the Bethe lattice. ** Since the underlining Bethe lattice has a coordination number increasing with $m$, the degree (average number of children per node) of the random tree increases with $m$. The void mobility and hence the particle mobility also increase with $m$, realizing a facilitation phenomenon. Percolation theory has been applied in the configuration space [49], but is more often applied in the physical space [51]. 8. 8.

A higher mobility at regions with multiple voids is a facilitation phenomenon. Physically, at a temperature at which a single void is typically described by a finite tree, the void is trapped among a few meta-stable sites. With another void, its passing by perturbs the particle configuration so that the first void may have access to additional configurations, accounted for by additional tree-branches. This can untrap the first void temporarily. Such mutual facilitation of motions between the voids can lead to a mobile pair of interacting voids. This is analogous to the two-spin facilitation model of Fredrickson and Andersen [7], one of the first KCM [2]. Our theory provides a physical explanation of the facilitation rules. In contrast, trap models [52] consider untraping via surmounting energy barriers by visiting also the energetically unfavorable configurations. They offer however slower untrapping dynamics negligible compared with the facilitation mechanism. 9. 9.

Super-Arrhenius dynamics and kinetic arrest follow from a sequence of mobility transition of increasing $m$ upon cooling. As $T$ decreases, regions with one, two, three voids and so on successively enter the immobile phase resulting at dramatic slowdown. Our theory is thus analogous to the $p$-spin facilitation model [7], with $p$ increasing upon cooling.

References