text:
# A consequence of local equilibration and heterogeneity in glassy materials

## I Generalized FDT: what do we know?

In this section, we introduce the basic concepts and quantities of interest for this paper. The exact results obtained in infinite-range glass models are first presented. We then review the argument given to support a possible larger validity of the mean-field results, and conclude the section with the experimental and numerical evidence available. The goal of this section is also to recall what is known, what is only supposed, and most importantly to emphasize what is not known.

### I.1 Infinite-range glass models

The solution of the non-equilibrium dynamics of infinite-range glass models consists of the asymptotic analysis of coupled dynamical equations involving two-time correlation functions, $`C(t_1,t_2)`$, and their thermodynamically conjugated response functions, $`R(t_1,t_2)`$ cuku . Alternatively, integrated responses, or susceptibilities, $`\chi (t_1,t_2)=_{t_2}^{t_1}dt^{}R(t_1,t^{})`$, can be studied, with the advantage that $`\chi (t_1,t_2)`$ is easier to measure in a simulation or an experiment. Adopting notations appropriate to magnetic systems, those functions read for instance

$`C(t_1,t_2)`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}s_i(t_1)s_i(t_2),`$ (1)
$`R(t_1,t_2)`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{s_i(t_1)}{h_i(t_2)}}|_{h=0},`$ (2)
for a system composed of $`N`$ spins $`s_i(i=1,\mathrm{},N)`$. The fields $`h_i`$ are thermodynamically conjugated to the spins. Since (1) and (2) are self-averaging, no average of any kind is required if the thermodynamic limit, $`N\mathrm{}`$, is taken.

It turns out that in infinite-range glass models, a generalized form of the FDT is satisfied in the asymptotic limit of large times in the aging case, or small driving forces for the driven case. The limits respectively read $`t_1,t_2\mathrm{}`$ with $`C(t_1,t_2)`$ fixed, or $`\epsilon 0`$ with $`C(t_1t_2)`$ fixed, where $`\epsilon `$ generically refers to the amplitude of the driving force (see below for examples). This asymptotic form of the FDT reads

$$R(t_1,t_2)=\frac{X[C(t_1,t_2)]}{T}\frac{C(t_1,t_2)}{t_2},$$
(3)
where $`T`$ is the temperature, and $`X(C)`$ is called the fluctuation-dissipation ratio (FDR). The equilibrium FDT is recovered when $`X(C)=1`$. The finding of a generalized FDT for non-equilibrated systems suggests the possibility to develop a non-equilibrium extension of standard statistical mechanics or thermodynamics. In particular, the quantity $`T/X`$ has been shown to play the role of a non-equilibrium effective temperature for the slow modes of the system cukupe , leading to an important body of related works reviewteff .

Moreover, in infinite-range glass models which statically exhibit a full replica symmetry breaking pattern, the new dynamical quantity $`X(C)`$ was found to be directly related to the static order parameter involved in the replica symmetry breaking, the distribution of overlaps $`P(q)`$, through the very simple relation cuku ,

$$X(C)=_0^Cdq^{}P(q^{}).$$
(4)
The ‘Parisi function’ is defined as

$$P(q)=\underset{N\mathrm{}}{lim}\delta \left(\frac{1}{N}\underset{i=1}{\overset{N}{}}s_i^\alpha s_i^\beta q\right),$$
(5)
where $`s_i^\alpha `$ denotes the value of spin $`i`$ in the configuration $`\alpha `$, and $`\mathrm{}`$ represents an average over configurations $`(\alpha ,\beta )`$ weighted by their Boltzmann probability.

Equation (4) has a remarkable form, since it relates two quantities of completely different origin. The FDR in the left hand side of (4) is measured in the off-equilibrium dynamics, while the Parisi function in the right hand side is computed statically using the equilibrium Gibbs measure. Note also that the computation of $`X(C)`$ does not involve replicas.

Finally, we recall the remarkable result that a generalized FDT as in Eq. (3) is asymptotically obeyed in the case of the driven dynamics cukupe ; BBKrheo , with the same value of the FDR, so that Eq. (4) remains unchanged in that case BBKultra .

### I.2 Stochastic stability

Further developments have taken a more speculative form, in the sense that no exact result has been obtained for finite dimensional models. However, it was argued that Eq. (4) is generally true for finite dimensional glassy systems silvio . We now briefly review the formal derivation of this result to emphasize the main hypothesis which are made. Following Refs. cuku , the generalized susceptibilities $`\chi ^r`$ are first introduced:

$$\chi ^r=\frac{r!}{N^{r1}}\underset{i_1<\mathrm{}<i_r}{}\frac{s_{i_1}\mathrm{}s_{i_r}}{h_{i_1\mathrm{}i_r}}|_{h=0}.$$
(6)
The ‘derivation’ of Eq. (4) simply consists of the comparison between a static and a dynamic computation of $`\chi ^r`$ silvio . A static average gives

$$\chi ^r_{\mathrm{stat}}=\underset{N\mathrm{}}{lim}\underset{t\mathrm{}}{lim}\chi ^r(t)=\frac{1}{T}_0^1dqP(q)(q^r1).$$
(7)
A dynamical average leads instead to

$$\chi ^r_{\mathrm{dyn}}=\underset{t\mathrm{}}{lim}\underset{N\mathrm{}}{lim}\chi ^r(t)=\frac{1}{T}_0^1dC\frac{\mathrm{d}X(C)}{\mathrm{d}C}(C^r1).$$
(8)
Equation (4) directly follows by requiring the equality $`\chi ^r_{\mathrm{stat}}=\chi ^r_{\mathrm{dyn}}`$ for all $`r`$.

There are, however, several assumptions in this derivation. First, one assumes that a non-trivial form of the FDT given by Eq. (3) is valid in the large time limit. Similarly, the static computation is performed in the thermodynamic limit, $`N\mathrm{}`$. A third assumption is that static and dynamic calculations give the same result. This can be justified by a standard nucleation argument silvio , although exceptions are known cuku . The last and strongest hypothesis has been called ‘stochastic stability’ silvio . The name originates from the fact that the susceptibility (6) probes the behaviour of the system under the random perturbation

$$H_{\mathrm{perturbation}}=\delta \underset{i_1<\mathrm{}<i_r}{}h_{i_1\mathrm{}i_r}s_{i_1}\mathrm{}s_{i_r},$$
(9)
in the limit $`\delta 0`$. Here, the field $`h_{i_1\mathrm{}i_r}`$ is a random Gaussian variable of mean 0 and variance $`r!/(2N^{r1})`$. Stochastic stability implies thus that the limits $`\delta 0`$ \[or equivalently $`h0`$ in Eq. (6)\], $`N\mathrm{}`$, and $`t\mathrm{}`$ can all be inverted. This last assumption is completely uncontrolled, so that Eq. (4) is made physically plausible by this calculation, but has no ‘theorem’ character. A measurable test of this last hypothesis in a realistic system is to check the relation (4) it implies, thus making the argument a bit circular. Interesting examples and counterexamples were discussed in detail in Ref. silvio .

### I.3 Numerical and experimental results

Since the generalization of the FDT through Eqs. (3) and (4) cannot be rigorously proved, numerical simulations and experiments are necessary in order to make some progress. The advantage of simulations is that both equations can be separately checked, while experiments can obviously not compare dynamical to static data, since there are no experimental way to directly measure the Parisi function. Its indirect determination through Eq. (4) would thus be a major result, as pointed out in Ref. silvio . It is the principal aim of this paper to discuss the issue of what is precisely determined from a dynamical measurement performed in a physically accessible time window.

Dynamic results have been numerically obtained in many different glassy models, as reviewed in Ref. reviewteff . More recently, experiments have also been able to independently probe two-time susceptibilities and correlations in a number of glassy materials manip1 ; manip2 ; manip3 ; manip4 . The standard way to present data and test Eq. (3) is to build an ‘FD plot’ of the susceptibility $`\chi (t_1,t_2)`$ as a function of the correlation function $`C(t_1,t_2)`$, parameterized by the time difference $`t_1t_2`$, conventionally taken as positive cuku . If Eq. (3) holds, one gets indeed

$$\chi (t_1,t_2)=\frac{1}{T}_{C(t_1,t_2)}^1dqX(q),$$
(10)
implying that the $`\chi (C)`$ relation is independent of time, and has a slope related to the FDR,

$$\frac{\chi }{C}=\frac{X(C)}{T}.$$
(11)
We have implicitly assumed that $`C(t_1,t_1)=1`$ to simplify the notations. It is important to notice that neither simulations nor experiments have ever reported an asymptotic FD plot, i.e., a time independent relation between $`\chi `$ and $`C`$. Quite strikingly then, physically accessible time scales are such that Eq. (3) is not valid. As a consequence, it is necessary to theoretically investigate preasymptotic behaviours where times are large, but finite. Equivalently for the driven dynamics, the case of small but finite driving forces has to be considered.

Far less works have investigated the validity of Eq. (4). Numerically, it is possible to compute the Parisi function in an equilibrium simulation, obtaining thus the function $`P(q,L)`$, where $`L`$ is the—necessarily finite—linear size of the system. The $`L`$-dependence is explicitly kept since it is again important to note that no numerical work has yet been able to report the asymptotic form of this function, $`P(q,L\mathrm{})`$, in a glassy system. Moreover, the system size will play an important role below.

Early studies have then empirically compared static and dynamic data through the comparison of the curves $`\chi (C,t_2)`$ obtained in the dynamics for times $`t_2`$ ‘as large as possible’, to the static curve $`S(C,L)`$ measured statically in a system of size $`L`$ ‘as large as possible’ fede2 ; fede3 . We have defined

$$S(C,L)=\frac{1}{T}_C^1dq_0^qdq^{}P(q^{},L).$$
(12)
With these definitions, Eq. (4) can be synthetically rewritten as

$$\chi (C)=S(C),$$
(13)
where $`S(C)=lim_L\mathrm{}S(C,L)`$.

The remarkable coincidence between static and dynamic data initially found in numerical simulations of finite dimensional spin glasses in dimensions $`d=3`$ and $`d=4`$ was first taken as the sign that both quantities had converged to their asymptotic limits fede2 . Conclusions on the correct description of the low-temperature phase in spin glasses phase were then drawn fede2 .

These conclusions were however premature, as strikingly demonstrated in recent experiments which have shown that even in the experimental time window the FD plots still have a clear time dependence manip3 . This again emphasizes the necessity of a careful theoretical study of preasymptotic behaviours in order to get a correct interpretation of experimental data. This issue is discussed throughout the rest of the paper.

## II Generalized FDT with space

### II.1 What about space? A scaling argument

To understand preasymptotic behaviours, a physical picture of the slow dynamical processes involved in the system is needed, since infinite-range glass models provide us with no relevant prediction in that regime. Unfortunately, due to their mean-field nature, they do not provide us with a clear description of the physics either. One has thus to resort to more phenomenological descriptions.

It is now well established, at least in the case of spin glasses, that slow dynamics is accompanied by the existence of a ‘coherence length’ scale, $`\mathrm{}`$: the larger the coherence length the slower the dynamics. Here, an analogy with coarsening phenomena in non-disordered ferromagnetic systems is quite illuminating. When the coarsening dynamics proceeds, larger and larger magnetized domains are present, and the dynamics indeed becomes slower and slower. The extent to which slow dynamics is always accompanied by a large coherence length in any glassy material is an important open problem in this field. In this paper, we adopt the point of view that this is indeed the case. See Refs. tarjus ; zouches  for related discussions.

Using this spatial description of slow dynamics, preasymptotic behaviours were conjectured in Refs. behose ; babe  to obey a generalization of Eq. (4). The physical picture is the following. At a given wait time, $`t_2`$, of the aging dynamics, the system has become locally equilibrated up to a coherence length scale, $`\mathrm{}(t_2)`$. It can thus be seen as a heterogeneous mosaic of independent and quasi-equilibrated sub-systems of size $`\mathrm{}(t_2)`$. Now, a dynamic measurement performed on this mosaic is in fact probing an ensemble average over equilibrated systems of finite size $`\mathrm{}(t_2)`$. Translating this idea into an equation, one arrives at the conjecture that the following generalization of Eq. (4) can be valid behose ; babe :

$$X(C,t_2)=_0^Cdq^{}P(q^{},L),L=\mathrm{}(t_2),$$
(14)
which implies that non-equilibrium properties of the system at time $`t_2`$, as encoded in the FDR, are related to the equilibrium properties of a system of finite size $`L=\mathrm{}(t_2)`$, as encoded in the Parisi function. Note also that static and dynamic preasymptotic effects are explicitly present in Eq. (14), since both finite times, $`t_2`$, and sizes, $`L`$, are taken into account.

The conjecture (14) was first formulated and tested in Ref. behose  in the context of the non-equilibrium critical dynamics of the 2$`d`$ XY model, and further studied in Refs. babe ; silvio2  in the numerical simulation of the Ising spin glass in dimension $`d=2`$. In both cases, the observed behaviours were by construction preasymptotic, since for both models the true asymptotic behaviour was simple equilibrium with $`X(C)=1`$ and $`P(q)=\delta (q)`$. However, on numerically accessible time scales, a non-trivial behaviour reminiscent of the one observed, say, in realistic spin glasses was found, and the conjectured relation (14) was indeed obeyed to a very good precision.

It is straightforward to obtain the relation corresponding to Eq. (14) in the case of a driven glassy material:

$$X(C,\epsilon )=_0^Cdq^{}P(q^{},L),L=\mathrm{}(\epsilon ),$$
(15)
where the relevant control parameter in the dynamics is now the amplitude of the driving force, $`\epsilon `$. The advantage of this formulation, as exemplified in the numerical simulations below, is that the dynamics is stationary, so that the coherence length itself does not change during the dynamical measurement, as is the case in the aging regime.

### II.2 Single-site quantities

A different generalization of Eq. (4) was recently discussed in the literature fede ; parisi . Inspired by numerical results obtained in a disordered spin system fede , the argument reviewed in section I.2 was reformulated using single-site quantities defined for a given sample parisi . Therefore, disordered systems only are concerned by this approach.

Consider the alternative single-site definition of a generalized susceptibility parisi ,

$$\chi _i^r=\frac{r!}{M^{r1}}\underset{a_1<\mathrm{}<a_r}{\overset{M}{}}\frac{s_i^{a_1}\mathrm{}s_i^{a_r}}{h_{a_1\mathrm{}a_r}}|_{h=0},$$
(16)
where $`M`$ thermodynamically coupled copies of the same system have been introduced, so that $`s_i^a`$ denotes the value of spin $`i`$ in the copy $`a`$. The role of a small coupling of amplitude $`k`$ between the $`M`$ copies is discussed in Ref. parisi . Its presence is necessary to properly define single-site overlap distributions $`P_i(q)`$, although its precise form is inessential.

If an inversion of the limits $`M\mathrm{}`$, $`t\mathrm{}`$, $`k0`$ and $`h0`$ is again allowed, then a reasoning analogous to the one developed in section I.2 leads to the equality

$$X_i(C)=_0^Cdq^{}P_i(q^{}),$$
(17)
between the single-site FDR and local overlap distributions. The single-site FDR is obtained by averaging local correlation, $`C_i(t_1,t_2)=s_i(t_1)s_i(t_2)`$, and response functions, $`R_i(t_1,t_2)=s_i(t_1)/h_i(t_2)`$, over various realizations of the thermal history for a single realization of the disorder. We note finally that the derivation of this local FDT makes use of the same type of (uncontrolled) hypothesis of the global one.

### II.3 Single-box quantities

Although Eqs. (17) and (14) are different, they are also consistent, because Eq. (17) physically follows from the concept of a local equilibration in space, so that dynamical properties at site $`i`$ can be linked to static properties at the same site $`i`$.

The disadvantage of this single-site formulation is however evident, since it entirely relies on the presence of quenched disorder in the system. A slight generalization of Eq. (17) would be to consider single-box quantities instead of single site to get similarly:

$$X_i^v(C)=_0^Cdq^{}P_i^v(q^{}),$$
(18)
where the box FDR $`X_i^v(C)`$ is defined from box dynamical functions, $`C_i^v(t_1,t_2)=v^1_js_j(t_1)s_j(t_2)`$ and $`\chi _i^v(t_1,t_2)=v^1_js_j(t_1)/h_j(t_2)`$, and the sums run over a box of finite volume $`v`$ centered around site $`i`$ malcom1 ; malcom2 ; malcom3 . The corresponding static box overlap distribution is defined as boxdata :

$$P_i^v(q)=\underset{N\mathrm{}}{lim}\delta \left(\frac{1}{v}\underset{j=1}{\overset{v}{}}s_j^\alpha s_j^\beta q\right),$$
(19)
using the same notations as in (5). The interest of the formulation (18) is that it should become independent of the considered site $`i`$ for a sufficiently large box volume $`v`$. A second interest is that it suggests an alternative to Eqs. (14) and (15), namely

$$X(C,t_2)=_0^Cdq^{}P^v(q),v=\mathrm{}^d(t_2),$$
(20)
for the aging case, and

$$X(C,\epsilon )=_0^Cdq^{}P^v(q^{}),v=\mathrm{}^d(\epsilon ),$$
(21)
for the driven case.

## III Simulations of aging and driven dynamics

We now turn to a numerical investigation of the two relations (14) and (15) discussed in the preceding section. We first present the model studied and some technical details, before giving our results both for a driven and an aging dynamics.

### III.1 Model and details of the simulations

In this section, we study numerically the Edwards-Anderson model of a spin glass defined by the Hamiltonian EA

$$H=\underset{i,j}{\overset{N}{}}J_{ij}s_is_j,$$
(22)
where $`s_i(i=1\mathrm{}N)`$ are $`N=L^d`$ Ising spins located at the sites of a cubic lattice when $`d=3`$, or a hypercubic lattice when $`d=4`$, of linear size $`L`$. The coupling constants $`J_{ij}`$ are random Gaussian variables of mean 0 and variance 1 and the sum in (22) runs over nearest neighbours. We use a standard Monte Carlo algorithm where the spins are randomly updated. Times are given in Monte Carlo steps, where one step represents $`N`$ attempts to update a spin. Sizes are given in units of the lattice spacing. To study the dynamics of the system, a large system size, $`L=40`$ ($`d=3`$) and $`L=24`$ ($`d=4`$), is chosen, so that few realizations of the disorder are needed, typically 10. Such large sizes are necessary, as the Fig. 1 discussed below will clearly demonstrate. This implies that conclusions drawn for systems as small as $`L=10`$ have to be taken with much caution star .

The aging dynamics is simulated by preparing the system at initial time in a random configuration, thus mimicking an infinite temperature state. The temperature is then changed at time $`t=0`$ to its final value. We will present below data for the temperature $`T=0.7T_c`$ in both $`d=3`$ and $`d=4`$, where $`T_c=0.95`$ and $`T_c=1.8`$, respectively. Dynamical measurements are then performed during the resulting aging process.

To simulate the driven dynamics, we use the same type of perturbation as in Ref. BBKultra . On each link $`(i,j)`$ of the lattice, a coupling $`\stackrel{~}{J}_{ij}`$ is added, where $`\stackrel{~}{J}_{ij}`$ is drawn from a Gaussian distribution of mean 0 and variance $`\epsilon `$, which defines then the amplitude of the driving mechanism. The key information about the $`\stackrel{~}{J}`$’s is that they are chosen to be antisymmetric, $`\stackrel{~}{J}_{ij}=\stackrel{~}{J}_{ji}`$, which implies that their effect cannot be incorporated in a Hamiltonian force of the type (22). Hence, the effect of the $`\stackrel{~}{J}`$’s is that of a non-conservative force, and the system is thus externally driven. As a result, a driven stationary state is reached at sufficiently large times, even for temperatures below $`T_c`$: aging is stopped by the driving force BBKultra . In that case, measurements are performed in the stationary regime, meaning that a considerable amount of the data corresponding to transient behaviours has to be discarded.

In both cases, dynamical measurements consist of the measurements of the global functions (1) and (2). Since these measurements are by now classic, we do not discuss them further here reviewteff .

In order to perform a comparison with static data, we shall make use of the published equilibrium data of Ref. pqdata3  in $`d=3`$ and Ref. pqdata4  in $`d=4`$. Fortunately, box-overlap data exist in $`d=3`$ and we shall also make use of those boxdata .

### III.2 Time is length

To accurately test the relations (14) and (15), the following steps must be followed.

* Identify and measure in the dynamics the coherence length $`\mathrm{}`$.
* Measure global correlation and response functions, from which the corresponding FD plot is built.
* Measure in equilibrium simulations the overlap distribution functions.
* Compare static and equilibrium data with carefully chosen values of the parameters (size, time, driving force), as suggested by Eqs. (14) and (15).

The first point, the identification of a coherence length, has been discussed for the aging dynamics of the Ising spin glass in various dimensions in Refs. heiko ; hajime ; BB , so that the existence, temperature and wait time behaviour of this quantity is well-characterized.

To the best of our knowledge, no such results are available for the driven dynamics. We have thus measured in the driven stationary state the standard four-spin correlation function heiko :

$$C_4(r)=\frac{1}{N}\underset{i=1}{\overset{N}{}}\overline{s_i^as_{i+r}^as_i^bs_{i+r}^b},$$
(23)
where the overline means an average over the disorder, and $`(a,b)`$ are two independent copies of the system. Some corresponding curves are shown in Fig. 1, for $`d=3`$. Very similar results are also obtained for $`d=4`$. The features observed in this figure are expected. The amplitude of the driving force is known to control the relaxation time, $`t_{\mathrm{rel}}=t_{\mathrm{rel}}(\epsilon )`$, of the system, the smaller $`\epsilon `$ the larger $`t_{\mathrm{rel}}`$ BBKultra . From Fig. 1, we recognize that a slower dynamics also implies the existence of a larger coherence length, $`\mathrm{}(\epsilon )`$, as revealed by a slower spatial decay of $`C_4(r)`$.

A non-ambiguous definition of the coherence length is obtained, as in the aging regime, through the study of the scaling properties of $`C_4(r)`$. For the aging, it is known that the correlation (23) is well described by the scaling form BB ; enzo

$$C_4(r)\frac{1}{r^{\alpha (T)}}𝒞\left(\frac{r}{\mathrm{}}\right),$$
(24)
which can be seen as the definition of the coherence length $`\mathrm{}`$. We find that the same scaling behaviour (24) is obeyed in the driven dynamics as well, with the same value of the temperature-dependent exponent $`\alpha (T)`$ as in the aging regime. As an example, we report the scaled data for $`T=0.7T_c`$ and $`d=3`$ in Fig. 2. We find that the scaling function $`𝒞(x)`$ is well described by a simple exponential form, $`𝒞(x)=\mathrm{exp}(x)`$. This contrasts with the ‘compressed’ exponential form reported in the aging regime BB . We have no explanation for this difference, which is however inessential.

The results of Figs. 1 and 2 allow us to obtain the relationship between the coherence length $`\mathrm{}`$ and the amplitude of the driving force, $`\epsilon `$. From the time dependence of the correlator (1), it is also possible to extract a relaxation time, $`t_{\mathrm{rel}}(\epsilon )`$. It is thus natural to eliminate $`\epsilon `$ to obtain the relationship between time and length. Not surprisingly, we find that the relation $`\mathrm{}(t_{\mathrm{rel}})`$ in the driven regime, is well compatible with the relation $`\mathrm{}(t_2)`$ measured in the aging regime at the same temperature.

Although this coincidence might seem anecdotic in this context, we emphasize its physical importance. This suggests indeed that whatever the control parameter for the dynamics is, waiting time or driving force, the dynamics is such that the same dynamic scaling $`\mathrm{}(t)`$ is obeyed: ‘Time is length’.

### III.3 Generalized FDT for the driven dynamics

Following the program presented in section III.2, we compute in the driven dynamics correlations and susceptibilities, and build the corresponding FD plots. Collecting the static data of Refs. boxdata ; pqdata3 ; pqdata4 , we are in a position to start and investigate the validity of the relations (14), (15), (20) and (21).

We present in Fig. 3 the results obtained in $`d=4`$. In this figure, there are three different FD plots, represented with points, and corresponding to three different values of the driving force, and thus to three different coherence lengths $`\mathrm{}(\epsilon )`$. The lines correspond to the double integral $`S(C,L)`$ defined in Eq. (12) of the Parisi function $`P(q,L)`$ for three different system sizes. The agreement between the two sets of curves is extremely good. That this is not a mere coincidence is supported by three important facts.

* The two types of measurements are performed entirely independently, in completely different physical situations;
* The agreement between the two sets of curves is very good on the whole range of the curve, not on a single point. If one remembers that each curve is built from a rather complex manipulation of various quantities without any fitting procedure or free parameter, the agreement becomes much more impressive;
* The coherence lengths measured in the dynamics are $`\mathrm{}(0.2)=1.55`$, $`\mathrm{}(0.16)=2.0`$, $`\mathrm{}(0.12)=2.5`$. They are thus in ratio 1:1.29:1.61. The corresponding static length scales are in ratio 1:1.33:1.66. The agreement with the prediction (15) that both lengths should scale similarly is thus clearly excellent.

### III.4 Another observable: The link overlap

To test further the hypothesis of a strong link between off-equilibrium properties and static ones resulting from the local equilibration of the material, we focus in this section on a second pair of dynamical observables. Physically, this is motivated by the picture given in section II.1, where an off-equilibrium system was described as a mosaic of equilibrated sub-systems of size given by the coherence length. A logical consequence is that the link between static and dynamic quantities should exist for all physical observables. This is equivalent to the prediction and test of the unicity of the effective temperature in structural glasses cukupe ; sollich ; shear ; mayer .

To check this hypothesis, we have investigated the following quantities:

$`C_l(t_1,t_2)`$ $`=`$ $`{\displaystyle \frac{1}{2Nd}}{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \underset{j=1}{\overset{2d}{}}^{}}s_i(t_1)s_j(t_1)s_i(t_2)s_j(t_2),`$ (25)
$`R_l(t_1,t_2)`$ $`=`$ $`{\displaystyle \frac{1}{2Nd}}{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \underset{j=1}{\overset{2d}{}}^{}}{\displaystyle \frac{s_i(t_1)s_j(t_1)}{h_{ij}(t_2)}}|_{h=0},`$ (26)
where the prime in the second sum indicates that it is restricted to the nearest neighbours of the site $`i`$. The field $`h_{ij}`$ is thermodynamically conjugated to the observable $`s_is_j`$, and can thus be seen as a random perturbation of the coupling constants. From these two dynamical functions, a new FD plot can be built.

Correspondingly, the static quantity to use is the link-overlap distribution function, $`P(q_l)`$, defined by

$$P(q_l)=\underset{N\mathrm{}}{lim}\delta \left(\frac{1}{2Nd}\underset{i=1}{\overset{N}{}}\underset{j=1}{\overset{2d}{}}^{}s_i^\alpha s_j^\alpha s_i^\beta s_j^\beta q_l\right).$$
(27)
We use here the static data published in Ref. pqdata4

The comparison between the new sets of static and dynamic quantities is reported in Fig. 4, for the same parameters as in Fig. 3. It should not come as a surprise that the range of correlators covered by Fig. 4 is smaller than the one covered in Fig. 3. The long-time limit of the correlator is indeed non-zero in that case.

It is again clear that the agreement between static and dynamic data is very good. The point to emphasize is that for a given amplitude of the driving force, the same size $`L`$ as in Fig. 3 is used to compute static data, in agreement with the physical picture described above.

These results are a striking confirmation of the local equilibration hypothesis formulated in this paper.

### III.5 Using the box overlap distributions

Having established numerically the validity of the relationship between static and dynamic properties in $`d=4`$, we now turn to $`d=3`$. This will allow us to make use of the box-overlap distribution functions published in Ref. boxdata . To compute this function, the linear size of the total system was fixed to $`L=12`$, and the linear size of the boxes, $`B`$, was varied boxdata .

Our numerical results are reproduced in Fig. 5, where the FD plots built from the spin-spin dynamical functions (1) and (2) for various amplitudes of the driving force are compared to the static data obtained from the box-overlap distributions.

As for $`d=4`$, we find a very good agreement between static and dynamic data. The coherence length scales are in ratio 1:1.35:1.59:1.88, while the size of the boxes are in ratio 1:1.33:1.66:2, the coincidence being again satisfying.

We conclude that numerical simulations in $`d=3`$ and $`d=4`$ of the driven dynamics of the Ising spin glass nicely support the theoretical propositions (15) and (21).

### III.6 Generalized FDT in aging dynamics

We conclude the report of our numerical results with a brief mention of the results we have obtained in the aging regime for $`d=3,4`$. Similar results have already been published, as mentioned above, but the crucial coincidence between time (in the dynamics) and size (in the statics) was not noted, although it was supported by the data fede2 ; enzo . We have also discussed above the results obtained for $`d=2`$, when no spin glass phase exists babe ; silvio2 .

We give an example of our own data in the aging regime in $`d=3`$ in Fig. 6. There are four important points we want to make on these data.

First, we make use in this figure of the box-overlap distribution function, which was not the case in previous works. The total overlap distribution functions give however results which are similarly good.

Second, the dynamic data can be described using the same equilibrium quantity $`P(q,L)`$ as for the driven dynamics. This is a clear support of the analytical result that driven and aging dynamics should give rise to similar deviations from the FDT cukupe ; BBKrheo ; BBKultra .

Third, we have intentionally included data for a very short wait time, $`t_2=27`$, for which the agreement with Eq. (20) is not quite as good, to underline the fact that during the dynamical measurements the coherence length itself evolves with time. This was obviously not the case when the driven dynamics was considered. For this short wait time, for instance, the final time is $`t_1=10^5`$, so that we do not expect the FD plot to follow the line built from static data obtained at a given box size corresponding to the earlier time $`L=\mathrm{}(t_2)`$. Indeed, a clear deviation at large times (corresponding to low values of the correlator) is observed. This was anticipated and noted in Refs. behose ; babe , and was another motivation to focus also on the driven dynamics in this paper.

Fourth, we have discussed above the link between the coherence length $`\mathrm{}`$ and the system size or box size in terms of relative ratio. It is clear that an absolute comparison is not possible, since the coherence length has no ‘absolute value’, contrary to the box-size univocally defined in the static computation. However, the ‘large’ numerical factor between $`L`$ and $`\mathrm{}`$ reported in Ref. fede2  (about 4 in $`d=3`$) was an argument used against the theoretical proposition (14) formulated in this and previous works mezard . It can be seen from Fig. 6 that the factor between $`\mathrm{}`$ and $`B`$ is less than 2 when the box overlap distribution is used, thus clearly weakening this criticism.

## IV Conclusion

In this paper, we have given a scaling argument to support the proposition of Refs. behose ; babe  to generalize the relationship between the off-equilibrium FDR and the equilibrium overlap distribution to the physically relevant domain of finite times and sizes. Inspired by previous approaches to the asymptotic relation silvio ; parisi , we proposed the relations (14) and (15) between the FDR computed at any moment of the dynamics, and the equilibrium box overlap distribution function, where the size of the box to be used is imposed by the dynamical coherence length.

We have argued that these relations reflect the simpler notions of the local equilibration in space and the dynamical heterogeneity of the system, two ingredients which are inherent to systems with slow dynamics, independently of their space dimensionality, their dynamics, or the presence of disorder.

We have emphasized throughout the paper that the relations (14) - (15) were the physically relevant ones, as compared to the asymptotic form of Eq. (4). Hence, they have to be taken into account before drawing conclusions on the equilibrium glass phase one seeks to study fede2 . This remark, already made in Refs. behose ; babe , is repeated here because there recently appeared FD plots built from experimental dynamic measurements in spin glasses manip3 . As anticipated babe , this dynamical measurements are very similar to the ones observed in simulations, despite the enormous difference in time scales between experiments and simulations. This is because the growth of the coherence length is so slow in spin glasses that this difference in time scales reduces to a very modest factor in terms of the coherence length zouches ; BB ; surfing . This in fact implies that the ‘finite size effects’ reported in simulations, are present in experiments as well. This points towards the experimental irrelevance of the thermodynamic limit in spin glasses.

Finally, we made clear that the relations (14) - (15) had not the character of a theorem, meaning that more work is needed to assess their range of validity. In particular, coarsening systems are a simple counterexample of this relation. This is due to the fact that the dynamical behaviour is entirely dominated by topological defects which are absent in the equilibrium simulations, so that finite time dynamics and statics do not coincide behose ; babe ; 1D . There remains to be seen to what extent this alternative ‘defect’ description of glassy dynamics generally holds in systems like supercooled liquids or soft glassy materials. This is certainly a quite challenging issue.

###### Acknowledgements.

M. Mézard suggested to use the existing numerical results for the box-overlap distributions to weaken the criticism he had himself formulated mezard . I thank E. Marinari, J. Ruiz-Lorenzo, H. Katzgraber, P. Young who kindly provided me with their static data boxdata ; pqdata3 ; pqdata4 , and A. Barrat, P. Holdsworth, and M. Sellitto for early discussions and collaborations behose ; babe . J. Kurchan made important remarks on a previous version of the manuscript and I warmly thank him for this luxury service après-thèse. This work is supported by a European Marie Curie Fellowship No HPMF-CT-2002-01927, CNRS (France) and Worcester College Oxford. Numerical simulations were started at the PSMN in ENS Lyon (France) and finished on Oswell at the Oxford Supercomputing Center, Oxford University.


format: text

text:
# RNA denaturation: excluded volume, pseudoknots and transition scenarios

## Abstract

A lattice model of RNA denaturation which fully accounts for the excluded volume effects among nucleotides is proposed. A numerical study shows that interactions forming pseudoknots must be included in order to get a sharp continuous transition. Otherwise a smooth crossover occurs from the swollen linear polymer behavior to highly ramified, almost compact conformations with secondary structures. In the latter scenario, which is appropriate when these structures are much more stable than pseudoknot links, probability distributions for the lengths of both loops and main branches obey scaling with nonclassical exponents.

In recent years, considerable attention has been devoted to the problem of describing the formation of secondary structure (base pairing map) in single molecular strands of RNA Zuker (1989); Hofacker et al. (1994); Higgs (1996); Tinoco Jr. and Bustamante (1999); Bundschuh and Hwa (1999); Pagnani et al. (2000); Chen and Dill (2000); Higgs (2000). The solution of such a problem is regarded as an important step within the general program of understanding how structure is encoded in the primary sequence of biopolymers. By making use of some simplifications, like that of disregarding excluded volume effects or pseudoknots formation, some studies established the existence of a molten phase at relatively high temperatures for an RNA molecule in dilute solution Bundschuh and Hwa (2002a, b); Müller (2003). In this phase the inhomogeneities associated to a specific primary sequence should be irrelevant for the large scale behavior and should allow the coexistence of a very large number of different secondary structures of comparable free energy.

As the temperature $`T`$ increases, a long RNA molecule should pass from the molten phase to a regime in which secondary structures essentially disappear and the global behavior becomes that of a linear polymer chain in good solvent. Excluded volume should play a relevant role at such a denaturation transition. Indeed, there the entropic free energy gain associated to the formation of hairpins or of more complicated branched structures with loops is comparable with the corresponding base pair binding and staking energies, and depends crucially on the repulsive interactions. Recent studies have shown that the discontinuous nature and the universal features of double stranded DNA denaturation are determined by excluded volume interactions Causo et al. (2000); Kafri et al. (2000); Carlon et al. (2002).

To our knowledge, starting with the related pioneering work of de Gennes de Gennes (1968) on the statistics of branchings and hairpin helices in the periodic dAT copolymer, excluded volume effects were never fully taken into account in studies of RNA denaturation. This leaves open the problem of establishing the existence and of determining the possible character of this transition in the long chain limit. A realistic embedding of the system in space, taking into account excluded volume, is also a necessary condition for discussing pseudoknots and their consequences. Pseudoknots occur, e.g., when two loops locally bind to each other determining a deviation of the configuration from planar topology Higgs (2000). Normally they are not included in models of the secondary structure Nussinov and Jacobson (1980); Zuker (1989); Hofacker et al. (1994); Higgs (1996); Bundschuh and Hwa (1999); Pagnani et al. (2000), or are considered as a perturbation Orland and Zee (2002).

In this Letter we propose a model of the large scale conformational behavior of RNA in the high $`T`$ and molten phases. Although schematic, our model takes fully into account excluded volume and allows control of the effects of pseudoknots. While providing useful informations on the behavior of finite RNA chains, an extensive numerical analysis allows to draw precise scenarios for denaturation and the associated scaling regimes in different conditions.

At coarse-grained level we model a conformation of the RNA strand as a two-tolerant trail of $`N`$ steps on the face centered cubic (FCC) lattice 111 The FCC lattice allows for closed loops also of odd length and thus enriches their sampling. . This is a random walk in which no more than two steps are allowed to overlap on a single lattice bond, forming what we call a contact. This restriction takes into account the excluded volume. In addition, by giving an orientation to the trail, we impose that only pairs of antiparallel steps can form contacts, and whenever this happens a gain in energy $`ϵ<0`$ is counted. The orientation of the trail reflects the backbone directionality or RNA Higgs (2000). In our model a contact corresponds to a sequence of bound base pairs over a distance of the order of the persistence length of the RNA double helix. This persistence length can vary with $`T`$ near a denaturation transition. However, this effect should not matter for the large scale properties. Since our model is coarse-grained and we are not interested in $`T`$’s below those of the molten phase, we neglect also the heterogeneity of base pair interactions. Thus, $`ϵ`$ represents an effective, average parameter.

Figure 1 reports schematically two possible configurations of our model. In both (a) and (b) a diagram on the right summarizes the corresponding contact map. A bridge in the diagrams connects each pair of steps forming a contact. Bridges are numbered in order of appearance if one follows the trail orientation. A main bridge is a bridge which is not inscribed within other, larger bridges. Unlike (b), (a) shows a pseudoknot, indicated by the crossing of two bridges in the diagram. This crossing means that a step forming a loop overlaps with one outside the loop. In order to investigate the role of pseudoknots, we consider two variants of the model, which we refer to as I and II. While in model I configurations with pseudoknots are allowed, in addition to those without pseudoknots, in model II the former are forbidden altogether. The choice of attributing the same energy to all kinds of contacts is a simplification of model I. Computationally it would be awkward to attribute selectively a weaker binding energy to those contacts which form pseudoknots, as physically appropriate in most situations Tinoco Jr. and Bustamante (1999).

Thermodynamic quantities and canonical averages are defined in terms of the partition function $`Z=_w\mathrm{exp}(H(w)/T)`$. The sum extends to all allowed configurations $`w`$ with $`|w|=N`$ steps, and $`H=ϵN_c(w)`$, $`N_c(w)`$ being the number of contacts in $`w`$. For both models, we sampled configurations by a multiple Markov chain Monte Carlo procedure Tesi et al. (1996) using several ($`20`$) temperatures satisfying $`0ϵ/T3.5`$ 222 To increase the mobility of the Markov chain at low $`T`$, where branched structures are expected, in addition to the usual set of local and pivot moves (see Tesi et al. (1996) and references therein), we implemented a variant of the pivot move, in which only one arm of the branched structure is rotated. In order to obtain smooth plots, data were processed by means of the multiple histogram method Ferrenberg and Swendsen (1988). .

We first computed as a function of $`T`$ the specific heat of model I and II for different $`N`$. At a continuous conformational transition with crossover exponent $`\varphi <1`$ one expects a singular behavior $`C_{\text{max}}N^{2\varphi 1}`$ for the maximum of this quantity as $`N\mathrm{}`$ Vanderzande (1998); Orlandini et al. (2000). In case $`\varphi <1/2`$ such singularity does not imply a divergence. For both models we find no evidence of a diverging $`C_{\text{max}}`$. Hence, at this level we can only conclude that for both models the denaturation transition must be continuous and with $`\varphi <1/2`$, if it exists.

We also determined two geometrical radii of the configurations, namely the end-to-end distance, $`R_e`$, and the radius of gyration with respect to the center of mass, $`R_g`$. Multicritical phenomena theory Privman et al. (1991); Vanderzande (1998) has taught us that the ratios of the averages of such radii in the $`N\mathrm{}`$ limit are universal numbers characteristic of the different regimes involved in the transition. Fig. 2 shows plots of $`R_e^2/R_g^2`$ for different $`N`$. For model I the trend of the curves gives indication of a sharp transition at $`ϵ/T1.9`$. Indeed, for high $`T`$ the ratio approaches from below the universal value $`6.25(1)`$ appropriate for a polymer in the swollen, self avoiding walk (SAW) regime Li et al. (1995). On the other hand, at very low $`T`$’s the trail should fold in double structures with maximal number of contacts ($`N_cN/2`$), in which $`R_e`$ necessarily approaches zero. This explains the trend towards zero (from above) of the curves at low $`T`$. Remarkable is the accumulation of intersection points for $`ϵ/T1.9`$. These intersections mark a change of the trend of the curves for increasing $`N`$ and suggest the presence of a peculiar transition regime with universal ratio $`4.8`$. Hence, for model I there is clear evidence of a second order transition at $`ϵ/T1.9`$. For model II there is no similar indication: the intersections are pushed towards lower and lower $`T`$’s as one compares curves corresponding to pairs of increasing $`N`$ values (Fig. 2, inset). This means that the larger $`N`$, the deeper the SAW regime extends in the low $`T`$ region. The whole pattern suggests for model II a smooth crossover, not a transition. Further insight is provided by the study of some scaling properties. The radius of gyration is expected to scale as $`R_gN^\nu `$ for large $`N`$ Vanderzande (1998). For both models we observe that at high $`T`$ the determinations of $`\nu `$ at finite $`N`$, for $`N\mathrm{}`$ approach a value $`0.59`$ appropriate for a SAW in $`d=3`$ Vanderzande (1998); Li et al. (1995). For model I, at $`T`$’s sufficiently below the transition the $`\nu `$ estimates can be extrapolated to $`0.35`$ for large $`N`$. This indicates that the configurations are very close to compact in the low $`T`$, molten phase ($`\nu 1/d=1/3`$). For model II at very low $`T`$’s we extrapolate $`\nu 0.4`$ which is also not far from $`\nu =1/2`$, as expected for branched polymers Vanderzande (1998)333Two tolerant trails have been used to model collapse from SAW to branched polymer behavior Orlandini et al. (1992). See also P. Leoni and C. Vanderzande, to be published..

The different behaviors of the two models are due to the presence of pseudoknots in model I. Indeed, the fraction of sampled configurations with pseudoknots in model I is already substantial and increases with $`N`$ at high $`T`$. It reaches soon values close to $`1`$ near the transition and below. Pseudoknots correspond to the formation of extra binding contacts and thus can lead to more compact configurations with respect to the case of model II. These extra contacts trigger the sharp transition observed at $`ϵ/T1.9`$. For model II, if present, a transition should be located at much lower $`T`$’s, most likely below the range of applicability of the model 444 At low enough $`T`$ the molten phase is supposed to leave place to a glassy one, in which the inhomogeneities due to primary structure become relevant Higgs (1996); Bundschuh and Hwa (1999); Bundschuh and Hwa (2002a, b); Pagnani et al. (2000)..

The possibility of forming pseudoknots is a driving factor in tertiary structure formation Tinoco Jr. and Bustamante (1999). However, model I somehow overamplifies this factor, because it gives pseudoknot forming contacts an energy equal to that of the other contacts. In fact contacts forming pseudoknots should correspond most often to the weak binding of quite short portions of single strand loops, like for kissing hairpins. Longer bindings giving rise to pseudoknots are expected to be kinetically inhibited Bundschuh and Hwa (1999). One way to make the energies of contacts forming pseudoknots closer to those of the other contacts is to introduce sufficiently high concentrations of divalent metal ions, like $`\text{Mg}^{2+}`$ in solution Misra and Draper (1998); Tinoco Jr. and Bustamante (1999). On the other hand, in model II pseudoknot forming contacts would appear if one would look at the details of the configurations at somewhat more coarse-grained level. This means that this model can be interpreted as one similar to model I, but giving essentially zero energy to such contacts. Thus, it is reasonable to expect the description of model II to be most appropriate for not too low $`T`$’s, and expecially when, e.g., a low concentration of $`\text{Mg}^{2+}`$ ions in solution enlarges the stability gap between secondary structure and pseudoknot forming links Misra and Draper (1998); Tinoco Jr. and Bustamante (1999).

RNA denaturation corresponds to a substantial suppression with increasing $`T`$ of the highly ramified structure of loops and branches characterizing the molten phase. The analysis of the loops is feasible and particularly instructive in model II. For DNA, the distribution of the lengths of denaturated loops, corresponding to openings of the double helix, follows a power law whose exponent $`c`$ determines the character of the transition Kafri et al. (2000); Carlon et al. (2002). A simple example of loop in RNA is given by the closure of an isolated hairpin. In this case the loop is connected to the rest of the structure by a single branch of double steps in model II. Of course, more complicated situations may occur (thick loop in Fig. 1(b)). Even at high $`T`$ an extensive number of minor spike-like branches is present along the RNA backbone. Thus, loops with a fixed number of branches naturally have length distributions with rather short cut-off. Therefore, we decided to sample the length of all loops identified in the various configurations, irrespective of the number of outgoing double step branches. The various lengths can be obtained from the contact map of each configuration, by using a recursive algorithm that identifies all the loops inside each main bridge in the diagram. Another interesting quantity is the return time, i.e. the total number of steps comprised within a main bridge. In the assumed planar topology of model II this time is the total arc length corresponding to each departure of the configuration from a contact-free, linear polymer behavior. Probability distributions of the loop lengths, $`\mathrm{}`$, and of the return times, $`n`$, are plotted in Fig. 3 for different $`T`$’s and for $`N=800`$. At high $`T`$, after transients both distributions behave as power laws with approximately identical exponents: $`P_{\text{loop}}(\mathrm{})\mathrm{}^c_{\mathrm{}}`$, $`P_{\text{ret}}(n)n^{c_r}`$, with $`c_{\mathrm{}}c_r=1.1(1).`$ The peaks at small $`\mathrm{}`$ and $`n`$ in the distributions indicate that loops at high $`T`$ mostly occur within isolated small hairpins in model II. The identity of $`c_{\mathrm{}}`$ and $`c_r`$ means that almost all large bridges are also main bridges. Thus, the return time essentially coincides with the loop length at high $`T`$. At lower $`T`$, while $`P_{\text{loop}}`$ becomes shorter and shorter ranged for decreasing $`T`$, the behavior of $`P_{\text{ret}}`$ remains of power law type at large arguments. The value of $`c_r`$ remains stable and close to that estimated for $`ϵ/T=0`$. This means that as the RNA molecule enters deeper and deeper into the molten phase with developed secondary structures, the loops become shorter and shorter. On the other hand, the main branches departing from the contact-free backbone encompass all accessible length scales, as appropriate for a branched polymer. This could also explain why in this range of $`T`$ the exponent $`\nu `$ discussed above is not far from $`1/2`$, as for branched polymers Vanderzande (1998). The exponent $`c_r`$ obtained here definitely deviates from the mean field value $`3/2`$ Bundschuh and Hwa (2002a).

Summarizing, for model I we could establish the existence of a sharp denaturation transition which is due to the presence of pseudoknots. Unlike the melting of DNA, this is a second order transition. Model I applies to experimental situations in which the stability gap between secondary and tertiary folding is sensibly reduced, e.g., by a high concentration of $`\text{Mg}^{2+}`$ ions Misra and Draper (1998); Tinoco Jr. and Bustamante (1999). On the other hand, one can regard model II as a more adequate description of RNA when the same stability gap is large, e.g. with low $`\text{Mg}^{2+}`$ concentrations Misra and Draper (1998); Tinoco Jr. and Bustamante (1999). In this model there is no sharp transition and the denaturation occurs as a crossover from linear to branched-compact polymer behavior. The geometry of this crossover is well described by the distributions $`P_{\text{loop}}`$ and $`P_{\text{ret}}`$ and by their exponents, whose nonclassical values are a further consequence of excluded volume.

We acknowledge discussions with E. Carlon. The work was supported by MIUR-COFIN01 and by INFM-PAIS02.