An alternative, dynamic density functional-like theory for time-dependent density fluctuations in glass-forming fluids

Slow dynamics in glass-forming systems is usually monitored through time-dependent density correlation functions.¹ The commonly used theoretical description of the density correlation function starts with the memory function representation, which absorbs the complicated dynamics of density fluctuations into a time-delayed friction kernel called the memory function.² The exact expression for the memory function involves a pair-density correlation function evolving with the so-called projected dynamics³ that is difficult to analyze or even simulate.⁴ To proceed, a factorization approximation is used, which expresses the pair-density correlation function evolving with projected dynamics as a product of two density correlation functions evolving with standard dynamics.³ The resulting closed set of equations constitutes the mode-coupling theory of glass-forming liquids,^3,5,6 which describes well several features of glassy dynamics. Its most often quoted success is the description of the cage effect, which manifests itself through a characteristic two-step decay of the density correlation function, with an intermediate-time plateau whose extension grows upon cooling. The theory reproduces well the wavevector dependence of the plateau and non-trivial power laws describing the approach to and the departure from the plateau.³ The factorization approximation can be pushed to higher-level correlation functions, leading to one of generalized mode-coupling theories.^7–9 Qualitatively, predictions of these theories are similar to those of the standard mode-coupling theory. However, these more sophisticated theories significantly improve upon the standard theory: they describe well the time-dependence of the density correlation function for a larger range of control parameters.⁹ 

One common feature of all mode-coupling-like approaches is the prediction of an ergodicity breaking transition, referred to as the mode-coupling transition, at which the plateau of the density correlation function extends infinitely. In two- or three-dimensional systems, the mode-coupling transition is replaced by a smooth crossover. The dynamics beyond the crossover is called “activated,”¹ but there is no commonly accepted quantitative description of the postulated activated processes.¹⁰ 

In addition to the fundamental problems of the avoided transition and the activated dynamics, the standard mode-coupling theory and its generalizations suffer from another important conceptual problem. Sometime after the mode-coupling theory was formulated, a number of static descriptions of the glass formation were proposed.^11–14 Most of these descriptions rely upon a free-energy-like functional derived using the replica method. They predict a transition that, on the basis of the analogy with the results obtained for exactly solvable spin-glass models, is referred to as the dynamic transition. At this transition, non-trivial replica off-diagonal density correlations appear, which are conceptually identified with the plateau of the time-dependent correlation function at the mode-coupling transition.^14,15 However, the location of the dynamic transition and the wavevector dependence of the replica off-diagonal density correlations are different from the location of the mode-coupling transition and the wavevector dependence of the infinite time plateau, respectively. In particular, in static theories of the glass formation, replica off-diagonal correlations are obtained by minimizing a free-energy-like functional. In contrast, the self-consistent equation for the wavevector dependence of the plateau derived within mode-coupling approaches cannot be obtained by a functional minimization. Thus, the two approaches are fundamentally incompatible.^16,52 This is a major conceptual problem since the static and dynamic (mode-coupling) approaches are considered to be the two facets of a unified description of glassy phenomena referred to as random first-order transition theory.^17,18

Recall that, in the limit of large spatial dimensions, a consistent description of the glass formation and glassy dynamics has been formulated.¹⁹ Both exact static²⁰ and dynamic²¹ theories have been developed and it was shown that the dynamic transition predicted by the static approach coincides with the mode-coupling-like transition predicted by the dynamic theory.¹⁹ 

Importantly, while there exist static approaches that with increasing spatial dimension reproduce the exact large dimensional static theory, in the same limit, the mode-coupling theory becomes drastically different from the exact large dimensional dynamic theory.^22,23 To rationalize this fact, we recall that the mode-coupling theory and the exact large dimensional dynamic theory, while similar in spirit, are nevertheless technically quite different. In particular, within the mode-coupling theory, one derives a self-consistent equation for the time-dependent density correlation function, whereas the exact large dimensional dynamic theory imposes a self-consistency condition on a stochastic process.²⁴ This condition cannot be expressed in terms of a self-consistent equation for a correlation function.

The above described unsatisfactory state of the finite-dimensional dynamic theory calls for a reformulation of the mode-coupling approach and/or for exploration of other approaches to the dynamics of glassy fluids. In this Communication, we propose a possible alternative approach. Our theory is consistent with static approaches that rely upon free-energy-like functionals. We identify a functional generalization of the Franz–Parisi potential,^25,26 which we call the Franz–Parisi functional, as the free-energy-like object generating the time dependence of the density correlation function. The equation of motion that we derive is similar in spirit to the equation of motion derived within dynamic density functional theory.²⁷ It is local in time and thus it does not capture the time-delayed friction that is considered the central feature of glassy dynamics. We believe that the present equation is just the first ingredient of a new dynamic theory and the time-delayed friction can be recovered if an analog of the memory function is identified.

In the following, we present a somewhat heuristic derivation of our theory. A more formal, projection operator-based, derivation and incorporation of time-delayed friction are left for future work.

We consider a system of N particles in volume V, with pair-wise additive interactions determined through spherically symmetric potential V(r). We assume that the microscopic dynamics is Brownian.^28,53 We denote the diffusion and friction coefficients of an isolated particle by D₀ and γ, respectively, with D₀ = T/γ, where T is the temperature and the Boltzmann constant is set to 1.

We aim to develop a theory for the density correlation function in the Fourier space, F(k; t),

where $n(k)=∑je−ik⋅rj$ is the Fourier transform of the microscopic density and Ω is the N-particle evolution operator, i.e., the Smoluchowski operator,

where β = 1/T and F_j is the force on particle j, F_j = ∑_l≠jF(r_jl) with $F(rjl)=−∂rjV(rjl)$⁠. Finally, $⋯$ in Eq. (1) denotes the equilibrium average; we use the convention that the probability distribution stands to the right of the quantity being averaged, and all operators to its left act on it as well as on everything else.

The standard approach to develop a theory for the density correlation function starts with the projection operator manipulations.²⁹ Below, we present an equivalent formulation of the first step of the standard approach, which avoids an explicit introduction of a projector operator. This formulation is inspired by the so-called linear kinetic theory developed to describe time-dependent equilibrium correlation functions.^30–32 We will contrast this approach with our new theory in Sec. IV.

We start by recognizing that the right-hand side of Eq. (1) can be interpreted as the density of a N-particle system that has time-dependent “probability distribution” of the following form:

where P_eq is the equilibrium probability distribution. We used quotation marks to emphasize that strictly speaking P_l(t) is not a probability distribution; in particular, it is not properly normalized. However, as explained by Résibois and Lebowitz,³⁰ “formally, this makes no difference.” Subscript “l” in P_l(t) indicates that distribution (3) has the same form as the linear change of the N-particle distribution due to an external potential that is turned off at the initial time, t = 0.

Next, we write the equation of motion for the density correlation function,

where $⋯l$ indicates averaging with time-dependent distribution P_l(t) [Eq. (3)]. We recall that following our convention, the evolution operator acts on distribution (3), which stands to the right of Ω.

To calculate the right-hand side of Eq. (4), we assume that at later times, t > 0, distribution P_l(t) has the same form as at the initial time but with a different magnitude,

where superscript “a” indicates the approximate character of distribution $Pla$ and function f(t) is chosen in such a way that averaging with $Pla(t)$ reproduces the density correlation function at time t,

where $⋯la$ indicates averaging with $Pla(t)$⁠.

A simple calculation gives f(t) = F(k; t)/S(k), where S(k) is the static structure factor, $S(k)=N−1n(k)n(−k)$⁠, and approximate distribution (5) depends linearly on the density correlation function at time t, $Pla(t)=n(−k)F(k;t)/S(k)Peq$⁠.

Using distribution $Pla$ in Eq. (4), we obtain the following approximate equation of motion for the density correlation function:

The same equation is obtained using the first step of the standard projection operator procedure, i.e., when one neglects the memory function term. We note that the equation of motion (7) is local in time and predicts that different Fourier components of F(k; t) relax independently. Only after the memory function is included, the relaxation of different Fourier components (modes) becomes coupled. This fact was the motivation for the name “mode-coupling” given to the first theory of this type, which was developed to describe critical dynamics.³³ 

We propose the following alternative approach to approximately evaluate the time evolution of the density correlation function. First, we rewrite definition (1) of the density correlation function by distinguishing between averaging over the time evolution (which corresponds to averaging over noise in the Langevin formulation of Brownian dynamics) and averaging over the initial positions of the particles that will be denoted by $rj0$⁠, j = 1,…, N,

where $⋯n$ denotes averaging with the following time-dependent distribution:

where ∑_π denotes the sum over permutations of particle labels, π(j) is the label of particle j in permutation π, and subscript “n” emphasizes that distribution P_n(t) cannot be interpreted as a linear change from the equilibrium distribution. We note that P_n(t) is the distribution of particles’ positions at time t, which are denoted by r_j, j = 1,…, N, but it depends parametrically on the initial positions $rj0$⁠, j = 1,…, N. Furthermore, in Eq. (8), n⁰(−k) denotes the microscopic density calculated for the initial configuration, $n0(k)=∑je−ik⋅rj0$⁠, and $⋯0$ denotes averaging over the equilibrium distribution of the initial positions of the particles.

Equations (8) and (9) allow for the following procedure to evaluate the density correlation function. First, one calculates average density, $n(k)n$⁠, using distribution (9). This average density depends on the initial positions of all the particles. Then, one calculates the joint average of $n(k)n$ and density of the initial configuration n⁰(−k) over the equilibrium distribution of the initial positions.

The advantage of this formulation is that to describe the evolution of density $n(k)n$⁠, one can (in fact, one should) go beyond the expansion used in the linear kinetic theory. In other words, the approximate theory for $n(k)n$ should be nonlinear.

Since one has to calculate the time evolution of the average density, one could try using the dynamic density functional theory.²⁷ However, one cannot use the standard version of this approach^34–36 since the important feature of average density $n(k)n$ is that it depends on the initial positions of all the particles, i.e., on N parameters, $rj0$⁠, j = 1,…, N.

Our proposed approach is similar in spirit to that used by Dufty and Rodríguez³⁷ to elucidate an earlier result due to Hauge.³⁸ Hauge showed that one can obtain long-time mode-coupling contributions to time-dependent equilibrium correlation functions from the nonlinear Boltzmann equation. This was unexpected since mode-coupling contributions were thought to originate from correlated sequences of particles’ interactions, whereas the Boltzmann equation was known to include only uncorrelated collisions. Dufty and Rodríguez pointed out that Hauge’s result could be explained if his one-particle density were re-interpreted as a density that implicitly depends on the positions of all the particles. They showed that, for the hard-sphere system, this new density satisfies exactly a nonlinear equation that has the same form as the Boltzmann equation. Once the time dependence of the new density is evaluated, its correlation with the density of the initial configuration reproduces the exact time-dependent equilibrium correlation function. Our average density $n(k)n$ is analogous to the one-particle density introduced by Dufty and Rodríguez.

To proceed, we write down the equation of motion for the density correlation function,

To evaluate the right-hand side of Eq. (10), we need an approximate expression for P_n(t).

At this point, we recall the central assumption of the dynamic density functional theory, which is known as the adiabatic approximation.^27,35 It states that the correlations in a non-equilibrium system are the same as those in an equilibrium system with a non-uniform density equal to the instantaneous density of the non-equilibrium system. We follow this approximation in spirit and assume that non-equilibrium correlations embodied in distribution (9) are the same as in an equilibrium system in which averaged density averaged over the initial conditions is uniform and the correlation of the averaged density with the density of the initial configuration reproduces the time-dependent density correlation function.

Specifically, we propose the following approximate formula for P_n(t):

where $Z[V1ext,V2ext]$ is the partition function,

and where one- and two-body time-dependent external potentials, $V1ext$ and $V2ext$⁠, are determined by the two conditions described above. Explicitly, we require that density $n(k)na$⁠, where $⋯na$ denotes averaging with distribution (11), is on average uniform,

Second, we require that the correlation of density $n(k)na$ and the density of the initial configuration reproduce the density correlation function at time t,

We assume that conditions (13) and (14) uniquely determine potentials $V1ext$ and $V2ext$⁠.

Physically, two-body potential $V2ext$ is the interaction necessary to maintain the correlation between the state of the system at time t and the initial configuration that is equal to F(k; t) while keeping the uniform average density, which is maintained by the additional one-body potential. We emphasize that approximate instantaneous distribution (11) depends on the density correlation function at time t in a nonlinear and complicated way.

We note that in the first step of the two-step averaging process, i.e., while evaluating $⋯na$⁠, the initial positions $rj0$⁠, j = 1,…, N, play the role of the quenched variables. Within our theory, they appear in a very natural way.

We propose the name quasi-equilibrium approximation for formula (11) to emphasize that the distribution of particle positions at time t is the same as in an equilibrium state that satisfies conditions (13) and (14). The additional motivation for this name is the conceptual similarity of our approximation with quasi-equilibrium construction for the long-time glassy dynamics proposed by Franz et al.^39,40 We comment on this point in Sec. V.

The final step is the evaluation of the average at the right-hand side of Eq. (10) with the approximate distribution (11) in terms of reduced distribution functions and external two-body potential $V2ext$⁠. After some work, we arrive at the following equation of motion:

The right-hand side of Eq. (15) is written in terms of two objects, a three-body average $n(k1)nan0(−k)n0(k−k1)0$ and the Fourier transform of two-body potential $V2ext$⁠. Both of these objects are functionals of the instantaneous value of the density correlation function, F(k; t). We discuss two possible technical approximations for the three-body average $n(k1)nan0(−k)n0(k−k1)0$ in  Appendix A. We show that, for non-interacting particles, Eq. (15) reproduces the known exact result in  Appendix B.

According to Eq. (15), the driving force for the time evolution of the density correlation function is the interaction necessary to maintain instantaneous correlations between the density at time t and the initial density. We recall that, at a dynamic transition of a static theory of the glass formation, the correlation between the density of the system and the so-called template (also known as the zeroth replica) appears spontaneously, without any system–template interaction. This implies that, at the dynamic transition of a static theory, the right-hand side of Eq. (15) vanishes and the density correlation function freezes. In other words, Eq. (15) predicts an ergodicity-breaking transition that coincides with the dynamic transition of the static theory.

To further develop the connection with static glass formation theories. we note that partition function (12) can be used to introduce a free-energy-like functional,

At a given time, functional $F$ depends on the initial positions of the particles, which play the role of quenched variables. It should be self-averaging with respect to the distribution of initial positions. By a Legendre transform of functional $F$ with respect to the one- and two-body potentials, one can obtain a functional that depends on the average density and average instantaneous value of the time-dependent density correlation function. The latter functional is a functional generalization of the Franz–Parisi potential.^25,26 We believe that two-body potential $V2ext$ can be obtained as a functional derivative of the Franz–Parisi functional with respect to the instantaneous density correlation function, at constant average density. Thus, the time evolution predicted by Eq. (15) stops when the Franz–Parisi functional reaches its local minimum.

Finally, to make explicit contact with the dynamic density functional theory, we recall that its equation for the time evolution of the non-uniform average density $n̄(k;t)$ can be written in a very similar way,⁴¹ 

where V^ext is the external potential needed to maintain non-uniform equilibrium density equal to the instantaneous average density $n̄(k;t)$⁠.

We proposed here an alternative theory for the relaxation of density fluctuations in glassy fluids and the glass transition. The main approximation of our theory is that the correlations between the state of the system at time t and the initial state of the system can be reproduced by a Boltzmann-like formula coupling these two systems in a quasi-equilibrium fashion. The central quantity that allows one to calculate the required coupling is a generalization of the Franz–Parisi potential, which gives more freedom to the coupling between the system and the template. Our equation of motion for the density correlation function is local in time but the relaxation of different Fourier components of the density correlation function is coupled. Our approach predicts an ergodicity-breaking transition identical to the dynamic transition predicted by a static theory of the glass formation based on the same Franz–Parisi functional. We believe that a numerical implementation of our theory should start with a specific approximate Franz–Parisi functional and use it to calculate the two-body potential needed to integrate Eq. (15). This task is left for future research.

It would be interesting to test the present approach on a spherical p-spin model for which the time-dependence of correlation functions can be analyzed exactly.

We envision two directions to extend our theory. First, the present theory uses an order parameter, density correlation function, that is uniform in space. The time evolution of this order parameter stops at the ergodicity-breaking transition, which corresponds to a local minimum of the Franz–Parisi functional. However, in finite dimensions, relaxation beyond the local minimum can happen via nucleation and growth processes. To investigate such processes, one needs to allow for inhomogeneous order parameter fields, in the spirit of Franz^42,43 and Dzero et al.⁴⁴ We also note that inhomogeneous order parameter fields were observed in numerical investigation of Cammarota et al.⁴⁵ They were also introduced within dynamic field theory developed by Rizzo.^46,47

Second, the local in time equation of motion implies that our theory misses time-delayed friction. This is in contrast to the mode-coupling approach, which prominently features non-local in time relaxation processes. We believe that a more formal derivation of our theory, based on a Kawasaki–Gunton-style projection operator,⁴⁸ can result in a generalization of the equation of motion (15), which will include a memory function term describing time-delayed friction. Since the projection operator term will involve time-dependent quasi-equilibrium distribution (11), we anticipate that the memory function expression will feature time-dependent vertices. We hope that as a result, the incorporation of the memory function will modify the relaxation near the ergodicity-breaking transition but will not change its location.

In the context of the non-local in time relaxation processes, we would like to comment on the relation of our approximation and the quasi-equilibrium construction of Franz et al.^39,40 We assumed that, to calculate the time derivative of the density correlation function, we can approximate the exact probability distribution (9) by an equilibrium distribution in an external potential that assures that the correlation between the density of the system at time t and its initial density is equal to F(k; t). In contrast, the starting assumption of Franz et al. is a quasi-equilibrium condition for the transition probability. Franz et al. showed that this assumption leads to non-local in time equations describing the long-time behavior of the density correlations in the vicinity of the plateau.

Finally, we note that the present theory shares some features with the co-called “naive” mode-coupling theory, which is the starting point of the nonlinear Langevin equation theory of activated hopping proposed by Schweizer and co-workers.^49–51 This theory’s equation of motion for the averaged order parameter is local in time, features a non-equilibrium free-energy-like function, and predicts an ergodicity-breaking transition at the point at which this function develops a local minimum. The main difference is that our order parameter is a collective quantity, whereas that of Schweizer and co-workers^49–51 is a single-particle quantity. In addition, we anticipate that, in our case, barrier crossing is facilitated by non-uniform order parameter fields, whereas the theory of Schweizer and co-workers^49–51 implicitly assumes a uniform, thermally facilitated barrier hopping.

I thank Elijah Flenner and Francesco Zamponi for comments on the manuscript. I acknowledge support of the NSF (Grant No. CHE 1800282).

The author has no conflicts to disclose.

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

In a numerical implementation of our theory, one would like to start with a scheme that uses only two-particle correlations. To this end, one would express three-body average $n(k1)nan0(−k)n0(k−k1)0$ in terms of the density correlation function and static pair correlation functions. Here, we propose two possible approximations that are linear in density correlation function. We note that the right-hand side on Eq. (15) is still a nonlinear functional of F(k; t) due to the presence of the two-body potential.

First, one can argue that at short times the two-body potential $V2ext$ is extremely short range, and therefore in approximate expression (11), a given r_j is correlated only with one of $rl0$s. This suggests that one may neglect off-diagonal terms in n⁰(−k)n⁰(k − k₁), which results in

Alternatively, one may resort to a convolution-like approximation,

where we used the equilibrium distribution of the initial positions, which implies $N−1n0(k)n0(−k)0=S(k)$⁠.

As a “sanity check,” which was also performed while deriving dynamic density functional theory,³⁴ we consider here Eq. (15) for non-interacting particles. In this case, the short-time approximation (A1) is exact and Eq. (15) in the direct space reads

where N(|r − r⁰|; t) − n is the inverse Fourier transform of F(k; t). Then, one can show that, for non-interacting particles, $βV2ext(|r−r0|)=−lnN(|r−r0|;t)$⁠, and thus Eq. (B1) reproduces the exact equation of motion for density correlation function of non-interacting particles. We note that the same equation of motion is obtained from Eq. (7) since for non-interacting particles S(k) = 1.