# Off-equilibrium dynamics in a singular diffusion model

## Abstract

We introduce a schematic non-linear diffusion model where density fluctuations induce a rich out of equilibrium dynamics. The properties of the model are studied by numerical simulations and analytically in a mean field approximation. At low temperatures and high densities we find a long off-equilibrium glassy region, where the system evolves out of an initially pinned state showing aging and a slow decay of the autocorrelation as an enhanced power law, along with strong spatial heterogeneities and violation of the fluctuation dissipation theorem.

As fluids are supercooled below the melting temperature their structural relaxation becomes very slow and may result in a glass transition characterized by a dramatic increase of the viscosity. Under these conditions the dynamical evolution of glass-forming liquids shows a markedly out of equilibrium behavior . In the main approach to the glassy dynamics, the Mode Coupling theory , the dynamical equations are solved by resumming a non trivial set of diagrams. This theory predicts an equilibrium relaxation time $`\tau _0`$ which diverges as a power law at a dynamical transition. Today it’s well established  that the quoted theory applies in a region located well before the point of structural arrest, where the relaxation time is found experimentally to diverge according to the Voghel-Fulcher law 

$$\tau _0\mathrm{exp}[v(\rho _c\overline{\rho })^1]$$
(1)   
Therefore we do not have information from the theoretical point of view on the kinetics close to the dynamical transition where standard theories do not apply. In this paper we introduce a schematic diffusion equation with a phenomenologically chosen mobility which reproduces the equilibrium relaxation time (1) observed in this region. We then study the consequences of such an assumption on the out of equilibrium dynamics. Due to its relative simplicity the model is amenable of both numerical and analytical investigations allowing a complete description of its features. Despite the simplicity of the equilibrium properties, the out of equilibrium behavior is complex and in qualitative agreement with the known properties of real systems. Our analytical results could in principle be tested quantitatively in experimental data and in molecular dynamics simulations.

Glassy systems are usually schematized as composed of particles rattling inside “cages” of typical size $`a`$ formed by the neighbors. Although diffusion is unimpeded inside the cells, motion over larger distances is strongly suppressed at high molecular densities because a global rearrangement of many particles is required. This glassy behavior is known to be reproduced also in the simple case of hard spheres . For suitably coarse graining of space and time scales, the dynamics of supercooled liquids is Brownian in its microscopic origin. Therefore we consider a diffusion equation for the variable $`\rho (\stackrel{}{r},t)`$ which represents the coarse grained particle density over distances of order $`a`$

$$\frac{\rho (\stackrel{}{r},t)}{t}=\left[M(\rho )\frac{\delta F\{\rho \}}{\delta \rho }\right]+\eta (\stackrel{}{r},t)$$
(2)   
In Eq. (2) $`F[\rho ]=𝑑\stackrel{}{r}\left[\rho \mathrm{ln}\rho +(1\rho )\mathrm{ln}(1\rho )\right]`$ is the entropy of the lattice gas-like model , which we consider for sake of simplicity. More realistic forms of $`F[\rho ]`$, which take into account the attractive interactions between particles can also be considered; here we find, however, that the qualitative features of the model do not change . $`M(\rho )`$ is a mobility, specified below, which is supposed to capture the main features of the constrained cooperative dynamics of a dense fluid. $`\eta (\stackrel{}{r},t)`$ is a gaussianly distributed random field, representing the thermal noise, whose expectations are given  by $`<\eta (\stackrel{}{r},t)>=0`$, $`<\eta (\stackrel{}{r},t)\eta (\stackrel{}{r}^{},t^{})>=2T\left\{M(\rho )\left[\delta (\stackrel{}{r}\stackrel{}{r}^{})\delta (tt^{})\right]\right\}`$, where $`<\mathrm{}>`$ means ensemble averages and $`T`$ is the temperature in units of the Boltzmann constant $`k_B`$.

Now we specify $`M(\rho )`$. As shown below, from Eq. (2) the characteristic equilibrium relaxation time $`\tau _0`$ behaves as $`M^1(\overline{\rho })`$, $`\overline{\rho }`$ being the average density, for low $`T`$. Then, in order to reproduce the behavior (1) of $`\tau _0`$, we assume a local mobility of the form 

$$M(\rho )=e^{v[\rho (\stackrel{}{r},t)1]^1}$$
(3)   
where a rescaled density, so that $`\rho _c=1`$, has been considered. Due to its generality, Eq. (2) is also suited for the description of different physical systems where a constrained cooperative dynamics is believed to play a fundamental role, such as granular materials . We have studied Eq. (2) for a temperature quench, by simulations and in a mean field approach. We sketch our numerical results before entering a detailed mean field analysis.

Eq. (2) has been simulated by a standard first order Euler discretization scheme on a $`128`$x$`128`$ two-dimensional square lattice with periodic boundary conditions, starting from an uncorrelated high temperature initial state. The system is quenched to a very low temperature (results will be presented for $`T=10^4`$, but similar behaviors are found for different temperatures).

The decay of the average density fluctuations $`S^2(t)=\sqrt{<(\rho \overline{\rho })^2>}`$, is plotted in Fig. $`1^b`$. For low densities $`\overline{\rho }`$ a normal liquid region is observed; $`S^2(t)`$ very quickly decays to the constant value characteristic of the equilibrium state. By raising the density, one enters a different region. Here the vanishing of $`M`$ slows the dynamics: this produces a long-lasting off-equilibrium glassy behavior before the equilibrium state is reached. The behavior of $`S^2(t)`$ shows that for large densities the dynamics can be divided in three regimes. Initially, for $`t`$ smaller than a characteristic time $`\tau _p`$, $`S^2(t)`$ remains constant. This is the first regime. An analysis of the system configuration $`\rho (\stackrel{}{r},t)`$ in this time domain does not show any appreciable evolution: the system is pinned in the initial configuration. Then, for $`t>\tau _p`$ a second regime is entered and less dense regions equilibrates whereas high density zones are still practically frozen. This is characterized by the decrease of $`S^2(t)`$. In this regime, that will be referred to as slow evolution, one observes pronounced correlated spatial heterogeneities in the system. This spatial pattern is outlined by a slow decay, as a function of $`k`$, of the structure factor $`C(\stackrel{}{k},t)=<\rho (\stackrel{}{k},t)\rho (\stackrel{}{k},t)>`$, that is consistent with a stretched exponential fit (see Fig. 2) $`C(\stackrel{}{k},t)𝒞\mathrm{exp}[l(t)k]^{2\mu }`$ with $`\mu 1/6`$ (at variance with the Gaussian decay of standard diffusion), similarly to some experimental observations . Finally the system enters the equilibrium state characterized by a constant value of $`S^2(t)`$. This whole pattern is reflected by the behavior of the particle mean square displacement $`R^2(t)`$, shown in Fig.$`1^a`$, calculated through

$$R^2(t)=_0^tD(t^{})𝑑t^{}$$
(4)   
where $`D(t)=<M(\rho )>`$ is the average mobility . For low densities $`R^2(t)t`$ in the whole time domain, as expected for simple diffusion. As the density is increased toward the limiting value $`\overline{\rho }=1`$ three regimes can again be distinguished. After an initial linear increase (regime 1) a progressively more pronounced inflection is observed in an intermediate time domain (regime 2) whose duration is enhanced as $`\overline{\rho }`$ is increased. The same pattern is observed in both spin-glass like lattice gas and Lennard-Jones molecular dynamics simulation . Asymptotically, in equilibrium, $`R^2(t)t`$, as for simple diffusion.

In order to get analytical results we now introduce an approximation on Eq. (2) by first expanding the logarithm on the r.h.s. of Eq. (2) to lowest order  and then by replacing the mobility $`M(\rho )`$ with the effective diffusivity $`D(\rho )=<M(\rho )>`$. Since average quantities do not depend on the position, due to space homogeneity, one has $`D(\rho )=D(t)`$. Eq. (2) then becomes

$$\frac{\rho (\stackrel{}{r},t)}{t}=D(t)^2\rho (\stackrel{}{r},t)+\eta (\stackrel{}{r},t)$$
(5)   
where the rescaling $`tt/[\overline{\rho }(1\overline{\rho })]`$ and $`T\stackrel{~}{T}=T\overline{\rho }(1\overline{\rho })`$, has been performed, and $`<\eta (\stackrel{}{r},t)\eta (\stackrel{}{r}^{},t^{})>=2\stackrel{~}{T}D(t)^2\left[\delta (\stackrel{}{r}\stackrel{}{r}^{})\delta (tt^{})\right]`$. Transforming Eq. (5) into reciprocal space, one obtains the following formal solution for the two time correlator $`C(\stackrel{}{k},t^{},t)=<\rho (\stackrel{}{k},t^{})\rho (\stackrel{}{k},t)>`$, ($`tt^{}`$)

$`C(\stackrel{}{k},t^{},t)`$ $`=`$ $`e^{[R^2(t^{})+R^2(t)]k^2}`$ (7)   
$`\left\{C(\stackrel{}{k},0,0)+\stackrel{~}{T}\left[e^{2R^2(t^{})k^2}1\right]\right\}`$   
The whole problem is now reduced to the knowledge of $`R(t)`$ which must be calculated self-consistently enforcing Eq. (4), where $`D(t)`$ is given by $`D(t)=<M(\rho )>=_0^1M(\rho )P(\rho )𝑑\rho `$. Here $`P(\rho )`$ is the probability distribution of the density field that, for Eq. (5) can be shown  to be Gaussian at all times. Then we have

$$D(t)=[2\pi S^2(t)]^{1/2}_0^1M(\rho )e^{(\rho \overline{\rho })^2/[2S^2(t)]}𝑑\rho $$
(8)   
The quantity $`S(t)`$ can be computed as $`S^2(t)=(2\pi )^d_{k<\mathrm{\Lambda }}C(\stackrel{}{k},t,t)𝑑\stackrel{}{k}`$ where $`\mathrm{\Lambda }`$ is a momentum cutoff of order $`a^1`$. From Eq. (7), $`S(t)`$ is a function of $`R(t)`$:

$$S^2(t)=h[S^2(0)q\stackrel{~}{T}]R^d(t)\mathrm{\Phi }_d[\sqrt{2}\mathrm{\Lambda }R(t)]+q\stackrel{~}{T}$$
(9)   
where $`\mathrm{\Phi }_d[x]=_0^xy^{d1}\mathrm{exp}(y^2)𝑑y`$, $`q=(\mathrm{\Sigma }_d/d)[\mathrm{\Lambda }/(2\pi )]^d`$, $`h=[d/(\mathrm{\Lambda }\sqrt{2})^d]`$ and $`\mathrm{\Sigma }_d`$ is the surface of the $`d`$-dimensional unitary hypersphere. Notice that the asymptotic value of the density fluctuations $`S^2(\mathrm{})=qT\overline{\rho }(1\overline{\rho })`$ vanishes at the point of structural arrest.

Eqs. (4, 8, 9) are a closed set of equations that can be studied analytically. In the present approximation the non linearity of $`M(\rho )`$ is accounted for by a self-consistency prescription for the calculation of $`D(t)`$. Similar approximation techniques are well developed and widely used in several field of statistical physics  producing reliable results . From Eq. (7) the normalized correlator $`\stackrel{~}{C}(\stackrel{}{k},t^{},t)=C(\stackrel{}{k},t^{},t)/C(\stackrel{}{k},t^{},t^{})`$ is given by

$$\stackrel{~}{C}(\stackrel{}{k},t^{},t)=e^{[R^2(t)R^2(t^{})]k^2}$$
(10)   
showing that a scaling form $`\stackrel{~}{C}(\stackrel{}{k},t^{},t)=𝒮[\varphi (t)/\varphi (t^{})]`$, with $`\varphi (t)=\mathrm{exp}\{R^2(t)k^2\}`$, is obeyed as suggested by a scaling approach to dynamical processes .

An important issue to understand the off-equilibrium dynamics is the relation between the response function to a small perturbing field $`h_\stackrel{}{k}`$, $`\chi _\stackrel{}{k}(t^{},t)_t^{}^t𝑑\tau \frac{\delta <\rho _\stackrel{}{k}(t)>}{\delta h_\stackrel{}{k}(\tau )}`$, and the correlation function in the unperturbed situation, $`C(\stackrel{}{k},t^{},t)`$. In equilibrium systems, where the fluctuation-dissipation theorem holds, the quantity $`X\stackrel{~}{T}\chi _\stackrel{}{k}(t^{},t)/C(\stackrel{}{k},t^{},t)`$ is equal to one. Out of equilibrium this relation is violated, and, generally, $`X`$ is a function of $`t^{}`$ and $`t`$: $`X=X(t^{},t)`$ . In the present mean-field approximation, one may interestingly show that the generalized “fluctuation-dissipation” ratio (FDR) $`X`$ is a function of the sole $`t^{}`$: $`X(t^{})=\{[C(\stackrel{}{k},0,0)\stackrel{~}{T}]\mathrm{exp}(2k^2R^2(t^{}))+\stackrel{~}{T}\}^1`$. Only if $`t^{}\mathrm{}`$, the usual version of the FDR with $`X=1`$ is recovered (notice that $`X1`$).

In the following we will report the main results of the mean field analysis referring to a longer publication  for all the details. From the solution of the model one sees that if the density $`\overline{\rho }`$ is small or the temperature $`\stackrel{~}{T}`$ is high one immediately enters the asymptotic stationary state that will be described later on. For high densities and low temperatures, on the other hand, the evolution remains markedly far from equilibrium for a long period and three dynamical regimes are found corresponding to different behaviors of $`R(t)`$ (see Figs. $`1`$), as discussed below.

Regime 1 - Pinning:  For short times, such that $`R(t)\mathrm{\Lambda }^1`$, as shown by Eq. (10), $`\stackrel{~}{C}(\stackrel{}{k},t^{},t)`$ is essentially constant since $`|k|<\mathrm{\Lambda }`$ and the system looks pinned. In this regime we also have $`\mathrm{\Phi }_d[\sqrt{2}\mathrm{\Lambda }R(t)]R(t)^d`$, consequently $`S(t)S(0)`$ and $`D(t)D(0)`$. Therefore $`R^2(t)D(0)t`$ (see inset Fig. $`1^a`$). The duration of this regime is $`\tau _p\mathrm{\Lambda }^2D^1(0)`$. Physically $`\tau _p`$ corresponds to the time the particle spends inside its cell (cage).

Regime 2 - Slow evolution: Pinning lasts up to $`\tau _p`$. For $`t>\tau _p`$, we have $`R(t)>\mathrm{\Lambda }`$, thus particles diffuse out of the cages and the evolution starts. For sufficiently long times, computing $`R(t)`$ through Eqs. (4, 8, 9) one finds $`R^2(t)b(\mathrm{ln}t)^\delta `$ (see inset Fig. $`1^a`$), where $`b=const.`$ and $`\delta =6/d`$. Eq. (7) implies that, for fixed $`t^{}`$ the correlator decays as an enhanced power law (see inset Fig. 2)

$$\stackrel{~}{C}(\stackrel{}{k},t^{},t)=\mathrm{exp}\{R^2(t^{})k^2\}\mathrm{exp}\left\{b[\mathrm{ln}(t)]^\delta k^2\right\}$$
(11)   
and $`S(t)(\mathrm{ln}t)^{3/2}`$ (see inset Fig. $`1^b`$). A logarithmic relaxation of the density fluctuations is also observed in Molecular Dynamics simulations of out of equilibrium liquid glass formers . When also $`t^{}>\tau _p`$, one has $`\stackrel{~}{C}(\stackrel{}{k},t^{},t)=\mathrm{exp}\left\{b([\mathrm{ln}(t)]^\delta [\mathrm{ln}(t^{})]^\delta )k^2\right\}`$. The characteristic duration time, $`\tau _e`$, of the slow evolution regime can be estimated , at low $`\stackrel{~}{T}`$, to be $`\tau _eM^1(\overline{\rho })`$.

Regime 3 - Stationary state: For long times $`t>\tau _e`$ a simple diffusive behavior is obtained because $`D(t)`$ always attains asymptotically a constant value $`D(\mathrm{})`$. This imply $`R^2(t)D(\mathrm{})t`$, as can be seen in Fig. $`1^a`$, so that the normalized correlator exhibits the usual exponential decay as a function of $`t`$: $`\stackrel{~}{C}(\stackrel{}{k},t^{},t)=e^{D(\mathrm{})tk^2}e^{\{R^2(t^{})k^2\}}`$ (see inset Fig. 2). When $`t^{}>\tau _e`$, we have $`\stackrel{~}{C}(\stackrel{}{k},t^{},t)=e^{D(\mathrm{})[tt^{}]k^2}`$ and time translational invariance is obeyed. In the small temperature limit the density fluctuations $`S(\mathrm{})`$ can be approximatively neglected  and $`D(\mathrm{})M(\overline{\rho })`$. This leads to $`\tau _0=M^1(\overline{\rho })`$ which gives Eq. (1), as previously stated.

So far we have studied the out off-equilibrium evolution of a system governed by Eq. (2) in the presence of a vanishing mobility for which Eq. (1) holds in equilibrium. We also want to consider the case in which the mobility vanishes as a power law, as is found for instance in the Mode-Coupling Theory of supercooled liquids : $`M(\rho )=[1\rho (\stackrel{}{r},t)]^\gamma `$. This relation implies an algebraic divergence of $`\tau _0`$: $`\tau _0=M(\overline{\rho })^1`$, as experimentally found in supercooled liquids in the temperatures or densities regions far away the ideal glassy transition. The different form of the mobility, as stated before, does not changes the global picture described so far with three different regimes. In the first and third regimes $`R^2(t)`$ is linear in $`t`$ as before, while in the second one we find, for $`\gamma >1`$, an anomalous diffusion $`R(t)wt^\beta `$, with $`w`$ = const. and $`\beta =4/(\gamma d+4)`$. In this regime we also find  a stretched exponential decay of the normalized correlator $`\stackrel{~}{C}(k,t^{},t)\mathrm{exp}\left\{R^2(t^{})k^2\right\}\mathrm{exp}\left\{wt^\beta k^2\right\}`$.

In this paper we have introduced a phenomenological equation for off equilibrium glassy dynamics. The only ingredients of the model are the diffusive behavior and the request of a Vogel-Fulcher (or algebraic) divergence of $`\tau _0`$, obtained by assuming a mobility as in Eq.(3). With these sole ingredients the out of equilibrium evolution of the model is observed to be highly non trivial, even in the mean field approximation which we have studied in details. Consistently with mean field theory, also the numerical integration of the full model shows the existence of a gradual crossover from a normal liquid  to a glassy behavior by raising the density. Some properties that are observed in systems close to the glassy transition, such as the existence of strong spatial heterogeneities, anomalous diffusion, slow decay and aging of density autocorrelations, are exhibited by the model. These predictions, as long as the non trivial fluctuation dissipation ratio $`X(t^{})`$, are all amenable of experimental check. In mean field this whole richness is observed in the preasymptotic off equilibrium dynamics (which, however, may be exponentially long), whereas the asymptotic equilibrium evolution is trivial. This is an important difference with real glassy systems, where a non trivial decay of $`\stackrel{~}{C}(\stackrel{}{k},t^{},t)`$ is also observed in equilibrium. However in mean field a non exponential decay of $`\stackrel{~}{C}(\stackrel{}{k},t^{},t)`$ can be ruled out on general grounds due to the Doob’s theorem . Further studies are in progress in order to characterize the complicate fluctuations occurring in the equilibrium state.

Acknowledgments We are grateful to M.Zannetti for interesting discussions and to S.Roux for valuable comments on the manuscript. F. C. thanks M. Cirillo and R. Del Sole for their hospitality in Rome university. This work was supported with the TMR network contract ERBFMRXCT980183 and by MURST(PRIN 97).