# Slow relaxation, dynamic transitions and extreme value statistics in disordered systems

Kristina van Duijvendijk, Grégory Schehr, Frédéric van Wijland

► To cite this version:

Kristina van Duijvendijk, Grégory Schehr, Frédéric van Wijland. Slow relaxation, dynamic transitions and extreme value statistics in disordered systems. Physical Review E: Statistical, Nonlinear, and Soft Matter Physics [2001-2015], 2008, 78 (1), pp.011120. 10.1103/PhysRevE.78.011120. hal-00284367v2

HAL Id: hal-00284367

https://hal.science/hal-00284367v2

Submitted on 3 Sep 2008

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Slow relaxation, dynamic transitions and extreme value statistics in disordered systems

K. van Duijvendijk,¹ G. Schehr,² and F. van Wijland ¹

¹Laboratoire Matière et Systèmes Complexes (CNRS UMR 7057), Université de Paris VII, 10 rue Alice Domon et Léonie Duquet, 75205 Paris Cedex 13, France

²Laboratoire de Physique Théorique (CNRS UMR 8627), Université de Paris-Sud, 91405 Orsay Cedex, France

We show that the dynamics of simple disordered models, like the directed Trap Model and the Random Energy Model, takes place at a coexistence point between active and inactive dynamical phases. We relate the presence of a dynamic phase transition in these models to the extreme value statistics of the associated random energy landscape.

PACS numbers: 05.40.-a, 75.10.Nr, 64.70.Pf

# I. INTRODUCTION

The dynamics of glassy systems is by definition out-of-equilibrium over experimental time-scales. Glassiness manifests itself through a great variety of dynamical features such as aging, non-exponential relaxation of correlation functions or super-Arrhenius slowing down of the dynamics [1]. However the definition of glassiness remains an open problem because, in many situations, no static parameters were found to indicate whether or not a system is in a glassy state. There is indeed a general agreement about the fact that the glassiness of a system does not necessarily arise from an underlying static transition [2].

To shed light on this question several dynamical approaches have been developed [1]. Here we explore the suggestion made in [3] that glassiness arises in a system when a coexistence between active and inactive regions of space-time takes place. The idea is that these *dynamical heterogeneities* are a defining feature of glassy systems. A method to inquire into the space-time character of a system is the application of the thermodynamic formalism of histories developed by Ruelle and coworkers [4]. While the equilibrium statistical formalism studies the fluctuations in the configuration space of the system, Ruelle’s formalism focuses on the time realizations the system can follow in configuration space. A central parameter in this formalism is the activity $K(t)$ of a history, which is the number of changes of configurations between the initial time, set to zero, and time $t$. This parameter is a physical, time-extensive observable.

For systems without quenched disorder, it was suggested that when for large times a system is found to have two well separated sets of histories – one where the activity is extensive in the system size, and the other where the activity is subex- tensive – there will be a coexistence of active and inactive regions of space-time in the system, sep- arated by sharp interfaces. Thus the system will exhibit slow dynamics. This was indeed success-

fully shown in the case of Kinetically Constrained Models of glasses [5] and for the non-equilibrium steady state of some Markov processes [6].

Here we investigate the space-time properties of a family of disordered systems. Besides the fact that these exhibit experimental features akin to those of structural glasses, there is some intrinsic interest to focus on the effect of quenched disorder. As will become clearer in the next sections, the central part of our study bears on the statistics of the activity $K(t)$, which is defined for each realization of disorder. From the probability distribution of $K$ one can build a mathematical object –its large deviation function– which in many respects can be argued to play the role of a dynamical free energy. As is known from the statics of disordered systems, the presence of quenched disorder induces additional difficulties when it comes to averaging over the disorder degrees of freedom. We have had to deal with similar ones for the dynamical free energy we are after in this work.

In this paper we focus on the directed Trap Model and the Random Energy Model. Although they are simpler to study they reproduce experimental features of more complex disordered and glassy systems like the super-activated behavior of the viscosity in glasses [7], aging [8, 9] or non-trivial violations of the fluctuation-dissipation theorem [10]. In both cases, we show analytically that the aforementioned large deviation function displays a discontinuity between the sets of active and inactive histories: this suggests that also in these disordered systems there is a signature of the glassy behavior using the thermodynamic formalism of histories. Furthermore we study the influence of the distribution of the disorder on the presence of a dynamic transition and discuss its occurrence in connection [11, 12] with the distributions of extrema of disordered energy landscapes.
## II. DIRECTED TRAP MODEL

A trap model is defined by $N$ independent traps labeled by an integer $i$, each trap being characterized by an energy $E_i$. We consider a continuous time Markov dynamics among traps: the dynamical evolution is specified by the probability $P_i(t)$ that the system stays in trap $i$ at time $t$, and by the transition rates $W_{ij}$ for jumping from trap $j$ to trap $i$. In the directed trap model the transition rates take the form

$$W_{ij} = \delta_{i,j+1} \frac{1}{B_j} \quad \text{where} \quad B_j = e^{-\beta E_j}, \quad (1)$$

$i = 1, \dots, N$. Thus the system evolves through the following Master Equation:

$$\frac{dP_i(t)}{dt} = -\frac{P_i(t)}{B_i} + \frac{P_{i-1}(t)}{B_{i-1}}. \quad (2)$$

The trapping times $B_j$ are random variables distributed according to a Lévy distribution

$$p(B_j) = \theta(B_j - 1)\mu B_j^{-1-\mu}, \quad \mu \in (0, 1), \quad (3)$$

where $\theta(x)$ is the Heaviside step function. In the following we will use the notation $\mathcal{O}(B_j) = \int p(B_j)\mathcal{O}(B_j)dB_j$ to denote an average over $B_j$. In particular, the distribution in Eq. (3) implies that the mean trapping time is infinite, which causes anomalous diffusion.

### A. Large Deviation Function

Ruelle's thermodynamic formalism enables to investigate the dynamics of a system by studying the histories the system can follow in the configuration space. In the directed Trap Model configurations are represented by traps so that a history in the configuration space between time 0 and time $t$ is specified by a sequence of traps visited during this time and by the time intervals elapsed between each jump from one trap to the next. To perform a statistics over histories one has to classify them according to a time-extensive (and history-dependent) parameter. A suitable one is the activity $K(t)$, defined as the number of configurations, here the number of traps, visited between time zero and time $t$. In the directed trap model, if the dynamics starts at time $t=0$ in the trap $i=0$, $K(t)$ is simply the trap in which the system stands at time $t$.

We now consider the probability $P_i(K, t)$ to be in trap $i$ at time $t$ at fixed activity, and define its Laplace transform:

$$\tilde{P}_i(s, t) = \sum_{K=0}^{\infty} e^{-sK} P_i(K, t). \quad (4)$$

It can be shown that $\tilde{P}_i(s, t)$ obeys an evolution equation (similar to but different from a Master Equation) of the form $\partial_t \tilde{P}_i(s, t) = \hat{W}^K \tilde{P}_i(s, t)$ where the elements of the evolution operator take the form:

$$\hat{W}_{ij}^K = \delta_{j,i-1} \frac{e^{-s}}{B_{i-1}} - \delta_{j,i} \frac{1}{B_i}. \quad (5)$$

From $\tilde{P}_i(s, t)$ in Eq. (4), one defines the cumulant generating function of the activity $K$:

$$Z_K(s, t) = \sum_i \tilde{P}_i(s, t) = \langle e^{-sK} \rangle, \quad (6)$$

where $\langle \dots \rangle$ stands for an average over all possible histories the system can follow between time 0 and time $t$. At large time $t$ one expects

$$Z_K(s, t) \sim e^{t\psi_K(s)}, \quad (7)$$

where $\psi_K(s)$ is the largest eigenvalue of the linear operator $\hat{W}_K$ in Eq. (5) and is thus a large deviation function. At large time $t$, the derivatives of $\psi_K(s)$ will give the cumulants of the activity $K(t)$. Note that $\psi_K(s)$ depends a priori on the realization of the random variables $B_i$ (3) and will thus be itself a random variable.

The parameter $s$ allows to probe the different histories the system can follow: positive $s$ will correspond to the inactive phase of the dynamics, i.e. $K$ smaller than its average, while negative $s$ correspond to the active phase of the dynamics, i.e. $K$ larger than its average. A discontinuity in the derivatives of the large deviation function will correspond to a dynamic transition between two different phases (the active and the inactive one), reflecting the sharp interfaces between dynamical heterogeneities.

### B. Dynamic Transition

To compute the large deviation function $\psi_K(s)$ we have to find the largest eigenvalue of $\hat{W}^K$ defined in Eq. (5). The $N$ eigenvalues $\lambda_j(s)$, $1 \le j \le N$ are solutions of the equation

$$e^{-s} \frac{\tilde{P}_{i-1}(s, t)}{B_{i-1}} - \frac{\tilde{P}_i(s, t)}{B_i} = \lambda_j(s) \tilde{P}_i(s, t), \quad (8)$$

and $\psi_K(s) = \max_{1 \le j \le N} \{\lambda_j(s)\}$. The computation of the characteristic polynomial is straightforward and leads to the eigenvalue equation

$$\frac{1}{N} \sum_{i=1}^{N} \ln(1 + \lambda_j(s) B_i) = -s. \quad (9)$$

Obviously, one has $\psi_K(s=0) = 0$. Let us first consider the case $s < 0$ where there is only one
$\lambda_j > 0$ solution of Eq. (9) which thus coincides with $\psi_K(s)$. In that case, i.e. in the active phase, it is natural to assume that $\psi_K(s)$ is self-averaging so that, in the limit $N \to \infty$, one expects $\psi_K(s) \simeq \overline{\psi_K(s)}$. Using the law of large numbers to treat the sum in Eq. (9), one obtains in the $N \to \infty$ limit

$$ \mu \int_{1}^{\infty} dB B^{-1-\mu} \ln(1 + \overline{\psi_K(s)}B) = -s, \quad (10) $$

which uniquely determines $\overline{\psi_K(s)}$. To describe the fluctuations of $\psi_K(s)$ around its mean for finite $N \gg 1$ one writes $\psi_K(s) - \overline{\psi_K(s)} = \chi_K(s)N^{-1/2} + O(N^{-1})$ where $\chi_K(s)$ is a zero-mean random variable. Using the Central Limit Theorem applied to the sum in Eq. (9), one obtains that $\chi_K(s)$ is a Gaussian variable such that the distribution of $\psi_K(s)$ is given for $N \gg 1$ by

$$ p(\psi_K(s)) = \frac{\sqrt{N}}{\sigma_\alpha \sqrt{2\pi}} \exp \left[ - \frac{N(\psi_K(s) - \overline{\psi_K(s)})^2}{2\sigma_\alpha^2} \right] \quad (11) $$

where $\sigma_\alpha^2$ is given by

$$ \sigma_{\alpha}^{2} = \left( \frac{B}{1 + \overline{\psi_{K}(s)}B} \right)^{-2} \times \left[ \ln^{2}(1 + \overline{\psi_{K}(s)}B) - (\ln(1 + \overline{\psi_{K}(s)}B))^{2} \right], \quad (12) $$

where $\overline{\psi_K(s)}$ is given in Eq. (10).

Let us now consider $s > 0$ and label the $B_i$'s such that $0 < B_1 < \dots < B_N$. In that case, one obtains the bounds

$$ -\frac{1}{B_N} < \psi_K(s) < 0. \quad (13) $$

From the distribution of the variables $B_i$'s in Eq. (3) one obtains that the one of the largest one $B_N$ is given by:

$$ p(B_N) = N\mu B_N^{-1-\mu} (1 - B_N^{-\mu})^{N-1}, \quad (14) $$

from which one gets for $N \gg 1$:

$$ -\frac{1}{B_N} \approx -N^{-\frac{1}{\mu}} \Gamma(1 + \frac{1}{\mu}) e^{-\frac{1}{\mu} + 1}. \quad (15) $$

Thus, from the bounds in Eq. (13) we find that $\lim_{N \to \infty} \overline{\psi_K(s)} = 0$ for $s > 0$. On the other hand, for $s < 0$, one has from Eq. (10) $\lim_{N \to \infty} \overline{\psi_K(s)} > 0$ and in that limit $\overline{\psi_K(s)} \propto (-s)^{1/\mu}$ for small $s$. Therefore, the first derivative $\overline{\psi'_K(s)}$ is continuous but higher order derivatives $\overline{\psi'_K(s)}$ with $p \ge 1/\mu$ will display a discontinuity in $s = 0$, indicating a dynamic phase transition (of order higher than one).

We have checked the presence of this dynamic transition by solving numerically the eigenvalue

equation (9). In Fig. 1, we show a plot of $\overline{\psi_K(s)}$ as a function of $s$ for different values of $N = 20, 30, 50$ and $100$ for $\mu = 0.8$. The data were obtained by averaging over $10^6$ samples. For $s > 0$ these numerical data show that $\overline{\psi_K(s)} \to 0$ when $N \to \infty$ (one finds indeed $\overline{\psi_K(s)} \propto -N^{-1/\mu}$, consistent with the bounds in Eq. (13)). On the other hand, this plot on Fig. 1 shows that $\lim_{N \to \infty} \overline{\psi_K(s)} > 0$ for $s < 0$. The solid line in Fig. 1 is the analytical value of $\overline{\psi_K(s)}$ obtained by solving numerically Eq. (10), which is in good agreement with the numerical data and shows the presence of a dynamical phase transition.

The slope of the large deviation $\overline{\psi_K(s)}$ function is related by Eqs. (6,7) to the mean value of the activity $\langle K \rangle$ over all histories, thus in the active phase ($s < 0$) of the dynamics the activity takes a constant value $\langle K \rangle > 0$ while in the inactive phase ($s > 0$) it is found to be subextensive in the system size ($\langle K \rangle \sim N^{-1/\mu}$). This coexistence of active and inactive phases of space-time, or dynamical heterogeneity is argued to manifest itself through the glassy properties of the directed Trap Model. In the corresponding pure system, where the waiting times are no longer random variables, glassiness obviously disappears, and so does the discontinuity in the derivatives of the large deviation function. The activity is a constant $\langle K \rangle > 0$ also in the $s > 0$ phase and no dynamical heterogeneity is present.

<figure><img src="image_2.png" /><figcaption>Figure 1: Plot of $\overline{\psi_K(s)}$ as a function of $s$ for different values of $N$ and $\mu = 0.8$. The solid line is the solution of Eq. (10): for $s > 0$ the curve approaches the horizontal axis as $N^{-1/\mu}$ for $N \gg 1$.</figcaption></figure>

The numerical results in Fig. 1 were obtained assuming that the large deviation function approaches its average in the large $N$ limit and is distributed according to Eqs. (11,12). Indeed this supposition was confirmed evaluating numerically the probability distribution $p(\psi_K(s))$ of $\psi_K(s)$ for a negative $s$. In Fig. 2 we show a plot
of $p(\psi_K(s))N^{-1/2}$ as a function of $N^{1/2}(\psi_K(s) - \overline{\psi_K(s)})$ for $\mu = 0.8$ and $s = -0.4$ ($N = 100,200$). The data were obtained by averaging over $3 \cdot 10^7$ samples. The solid line is the Gaussian distribution obtained rescaling by $\sqrt{N}$ the Eqs. (11, 12). This gaussian id approached by the numerical data as the system size grows.

On the other hand, for $s > 0$, the evaluation of the average large deviation function, i.e. the solution of Eq. (10), was done limiting the values of $\psi_K(s)$ to the $N$-dependent bound (13). Indeed we found numerically that the distribution $p(\psi_K(s))$ of $\psi_K(s)$ for a positive $s$ is well fitted by the form (14):

$$p(\psi_K(s)) = N\mu\psi_K(s)^{\mu-1}(1 - \psi_K(s)^{\mu})^{N-1} \quad (16)$$

for all $N$. In the inset of Fig. 2 we show the rescaled distribution $p(\psi_K(s))N^{-1/\mu}$ as a function of $\psi_K(s)N^{1/\mu}$ for $N = 100,200$ averaged over $10^7$ samples for $s = 0.2$ and $\mu = 0.8$. These distribution are found to be well approached by the large $N$ limit of the rescaled distribution which takes the form $p(x) = \mu x^{\mu-1} \exp(-x^{\mu})$.

<figure><img src="image_3.png" /><figcaption>Figure 2: Numerical evaluation of the distribution of the large deviation function for $\mu = 0.8$ and $s = -0.4$ (N=100,200). The solid line is the rescaling of the Gaussian distribution (11). Inset: numerical evaluation of the distribution of the large deviation function for $\mu = 0.8$ and $s = 0.2$. The solid line is the large $N$ limit of the rescaling of distribution (16).</figcaption></figure>

### III. RANDOM ENERGY MODEL (REM)

We consider now a widely explored model of disordered system, the Random Energy Model introduced in [13] as a mean field spin-glass model. Its dynamics has been studied in detail by Koper and Hilhorst in Ref. [14].

#### A. Dynamics of the REM

The Random Energy Model is a system of $2^N$ energy levels $E_i$ which are independent random variables distributed according to a Gaussian distribution of zero mean and variance $\frac{1}{2}\mathcal{N}J^2$:

$$P(E) = \frac{e^{-E^2/\mathcal{N}J^2}}{\sqrt{\pi\mathcal{N}J^2}} \quad (17)$$

The model exhibits a phase transition at the critical temperature $T_c = J/2\sqrt{\ln 2}$ below which the free energy becomes a constant $E_0(\mathcal{N}) \approx \mathcal{N}J\sqrt{\ln 2}$ in the thermodynamic limit and the entropy vanishes. As a consequence, the energy levels close to $E_0(\mathcal{N})$ will dominate the low temperature phase. We thus consider a system of $N$ energy levels

$$E_i = E_0(\mathcal{N}) + \epsilon_i, \quad (18)$$

where $\epsilon_i$ is small and nonextensive [14] distributed according to

$$p(\epsilon) = \begin{cases} \rho e^{\rho(\epsilon - \epsilon_c)}, & \epsilon \ge \epsilon_c, \\ 0, & \epsilon < \epsilon_c, \end{cases}$$

where $\rho = T_c^{-1}$ is a constant and $\epsilon_c$ is a cut-off energy. In the following, we will consider the scaling limit $N \to \infty$, $\epsilon_c \to \infty$, while $Ne^{-\rho\epsilon_c} = v$ is kept fixed, where the physical properties of the REM for temperatures $T < T_c$ take finite, cut-off independent values.

The dynamics of this model is defined by a Master Equation

$$\frac{dP_i(t)}{dt} = \sum_{j \neq i} W_{ij} P_j(t) - \sum_{j \neq i} W_{ji} P_i(t), \quad (19)$$

where the transition rates $W_{ji}$ for going from level $i$ to level $j$ are taken of the form [14]:

$$W_{ji} = B_j V_j V_i \quad \text{with} \quad B_j = e^{-\beta\epsilon_j}, \quad (20)$$

where $i = 1, \dots, N$ and $V_i$ are positive random variables. One can easily check that the rates $W_{ji}$ obey the detailed balance condition. The distribution of the $B_i$'s can be easily derived from the distribution of the energies $\epsilon_i$, leading to:

$$p(B) = \begin{cases} \frac{v}{N} B^{-1-\mu}, & \text{for } \left(\frac{v}{N}\right)^{1/\mu} < B, \\ 0, & \text{otherwise} \end{cases} \quad (21)$$

where we set $\mu = \rho/\beta = T/T_c$, $\mu \in (0, 1)$. Following Ref. [14], we will consider level dependent barriers $V_i = B_i^{-q}$ with $q \in (0, 1)$. If one associates to each level $i$ an independent random magnetization $\mu_i$ with zero mean and variance $N$, the kinetic REM can be considered as a magnetic model. Within this choice of level dependent barriers it has been shown that the behavior of the equilibrium
autocorrelation function for the magnetization exhibits a non-exponential behavior [14]

$$ \overline{C(t)} \propto t^{-2\mu/q} \quad \text{for } q + \mu > 1, \qquad (22) $$

$$ \overline{C(t)} \propto t^{-\eta} e^{-Dt^{\gamma}} \quad \text{for } q + \mu < 1 \text{ and } q < \frac{1}{2}, $$

$$ \overline{C(t)} \propto t^{-2\mu/(2q-1)} \quad \text{for } q + \mu < 1 \text{ and } q > \frac{1}{2}, $$

where $D$ is a constant, $\gamma = \mu/(1-2q+\mu q/(1-q))$ and $\eta = (2\mu/(1-q)-1)\gamma$.

Here we apply the thermodynamic formalism of histories to this model and show that a dynamic transition takes place reflecting the slow dynamics in this system (22).

## B. Dynamic transition

To investigate the presence of a dynamic transition, we look again at the large deviation function $\psi_K(s)$ defined above in Eq. (7), associated to the number of configuration changes $K(t)$ between time 0 and time $t$. From the Master Equation governing the dynamics (19), one obtains the equation of evolution of $\tilde{P}_i(s,t)$ defined in Eq. (4)

$$ \begin{aligned} \frac{d\tilde{P}_i(s,t)}{dt} &= \sum_j \hat{W}_{ij}^K \tilde{P}_j(s,t) \qquad (23) \\ &= e^{-s} B_i V_i \sum_{j \neq i} V_j \tilde{P}_j(s,t) - V_i \tilde{P}_i(s,t) \sum_{j \neq i} B_j V_j. \end{aligned} $$

The large deviation function $\psi_K(s)$ is the largest eigenvalue of the evolution operator $\hat{W}^K$. The eigenvalues $\lambda_j(s)$ of $\hat{W}^K$ are solutions of the equation

$$ f(\lambda_j(s)) = 1, \qquad (24) $$

$$ f(\lambda) = \sum_i \frac{e^{-s} B_i V_i^2}{\zeta V_i + (e^{-s} - 1) B_i V_i^2 + \lambda} \qquad (25) $$

where $\psi_K(s) = \max_{1 \le j \le N} \{\lambda_j(s)\}$ and $\zeta = \sum_i B_i V_i$.

We first focus on the case $q+\mu > 1$. One notices that $f(\lambda)$ has $N$ simple poles on the negative real axis at $x_i = -(\zeta V_i + (e^{-s}-1)B_i V_i^2) < 0$, and decreases to zero for $\lambda \to \infty$. One has also $f(0) = 1$ for $s=0$, whereas $f(0) > 1$ for $s > 0$ and $f(0) < 1$ if $s < 0$. Thus in the active phase, $s < 0$, there is only one positive eigenvalue $\lambda_j$ which coincides with $\psi_K(s) > 0$. As done previously for the directed trap model, we thus suppose that $\psi_K(s)$ is self averaging in the limit $N \to \infty$ and coincides with its average. One obtains that, for $s < 0$

$$ \overline{\psi_K(s)} \sim N^{(2q+\mu-1)/\mu} \tilde{\psi}_K(s), \qquad (26) $$

where $\tilde{\psi}_K(s)$ is independent of $N$ and solution of the equation

$$ \mu\nu \int_{v^{1/\mu}}^{\infty} \frac{B^{-\mu-2q}}{\tilde{v}B^{-q} + \tilde{\psi}_K(s)} dB = e^s, \qquad (27) $$

where $\tilde{v} = \mu v^{(1-q)/\mu}/(q+\mu-1) > 0$. Performing an analysis similar to the one done for the directed trap model, one also finds that the fluctuations of $\psi_K(s)$ around its mean value $\tilde{\psi}_K(s)$ are described by a Gaussian as in Eq. (11) of variance $N^{-1/2}$.

However, the analysis in the inactive phase is more subtle. Indeed, given the singularities of $f(\lambda)$ in $x_i$, together with the fact that $f(0) > 1$ for $s > 0$, $\psi_K(s)$ satisfies the bounds

$$ -\min_i (\zeta V_i + (e^{-s} - 1) B_i V_i^2) < \psi_K(s) < 0. \qquad (28) $$

If one labels the random variables $B_i$ such that $B_1 < B_2 < \dots < B_N$ one easily sees that the minimum in the left hand side of the inequality (28) is reached for the maximum $B_N$, whose distribution is given in the large $N$ limit by

$$ p(B_N) = v\mu B_N^{-1-\mu} e^{-vB_N^{\mu}}. \qquad (29) $$

Using this distribution (29) together with the distribution of the $B_i$'s in Eq. (21) one obtains that

$$ -c_1 N^{(q+\mu-1)/\mu} < \overline{\psi_K(s)} < 0, \qquad (30) $$

where $c_1 > 0$ is a constant, independent of $N$. To find an upper bound for $\overline{\psi_K(s)}$, one writes the equation above (24) in a different form:

$$ \sum_i \frac{e^{-s} B_i V_i}{e^{-s}\zeta + (1 - e^{-s})(\zeta - B_i V_i + \lambda(s)V_i^{-1})} = 1. $$

Since the first term in the denominator is the sum of the numerators, the remaining terms $(1 - e^{-s})(\zeta - B_i V_i + \lambda(s)V_i^{-1})$ must be negative for at least one $i$. In this way we get an upper bound for $\psi_K(s)$:

$$ \psi_K(s) < (e^{-s} - 1) \min_i (\zeta V_i - B_i V_i^2) \qquad (31) $$

Again, the minimum in the right hand side of this inequality is reached for the maximum $B_N$. Performing the average, one gets

$$ \overline{\psi_K(s)} < -c_2 N^{(q+\mu-1)/\mu}. \qquad (32) $$

where $0 < c_2 < c_1$ is a constant. Combining these two bounds (30, 32) one gets for $s > 0$

$$ -\overline{\psi_K(s)} = \mathcal{O}(N^{(q+\mu-1)/\mu}). \qquad (33) $$

Thus by looking at the behavior of $\overline{\psi_K(s)}$ for large $N$ in Eq. (26, 27, 33), one shows that there is a
dynamical transition as $s$ crosses 0. Eq. (27) shows that $\tilde{\psi}_K(s) \propto (-s)^{q/(q+\mu-1)}$ for small $s$ indicating that, as in the directed trap model, the order of this dynamical transition is larger than one.

These analytical predictions for $q + \mu > 1$ in Eq. (26, 33) have been confirmed by solving numerically the eigenvalue equation (24). The average value $\tilde{\psi}_K(s)$ was computed by averaging over at least $10^6$ realizations of the disorder. In Fig. 3, we show a plot of $N^{-(\mu+2q-1)/\mu} \overline{\psi}_K(s)$ as a function of $s$. In agreement with Eq. (26, 33), one observes a dynamical transition occurring at $s=0$. Moreover it has been checked that for positive $s$ the large deviation function behaves like predicted in Eq. (33), i.e. $-\tilde{\psi}_K(s) = \mathcal{O}(N^{(q+\mu-1)/\mu})$.

We have checked numerically that this transition is also present for $q + \mu < 1$, where the equilibrium correlation function displays a non-exponential decay (22). In Fig. 4, one shows a plot of $N^{-1/2}\tilde{\psi}_K(s)$ as a function of $s$ for $q = 0.45$ and $\mu = 0.45$. In that regime, the autocorrelation function decays like a stretched exponential and our numerical data clearly shows a transition occurring for $s = 0$. Similarly, in Fig. 5, one shows a plot of $\overline{\psi}_K(s)N^{-2(q-\mu)/\mu}$ as a function of $s$ for $q = 0.6$ and $\mu = 0.35$, where the correlation function decays like a power law. Here again, our numerical data clearly shows a dynamical transition occurring for $s = 0$. Again these results were obtained averaging over at least $10^6$ realizations of the disorder. Thus in the three cases we have an anomalous behavior of the correlation function (22) and a discontinuity in the derivatives of the large deviation function $\overline{\psi}_K(s)$ at $s = 0$ (Fig. 3-5).

<figure><img src="image_4.png" /><figcaption>Figure 3: Numerical evaluation of the large deviation function $\overline{\psi}_K(s)$ for $q = 0.9$ and $\mu = 0.8$ ($q + \mu > 1$) in the Random Energy Model.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Numerical evaluation of the large deviation function for $q = 0.45$ and $\mu = 0.45$ ($q + \mu < 1$ and $q < 1/2$) in the Random Energy Model.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Numerical evaluation of the large deviation function for $q = 0.6$ and $\mu = 0.35$ ($q + \mu < 1$ and $q > 1/2$) in the Random Energy Model.</figcaption></figure>

## IV. CONNECTION WITH THE EXTREME VALUE STATISTICS

The glassy behavior of disordered systems emerges in the low temperature phase or in the long time limit. Consequently the extreme values of the disorder realizations play a dominating role with respect to the typical values. The low temperature phase of the Random Energy Model is indeed described through the distribution of the lowest energy levels (18), which is found to belong to the Gumbel universality class of the minimum of variables which are unbounded and have a distribution that decays faster than any power law to $-\infty$.

It was pointed out in [11, 12] that the Gumbel class corresponds exactly to the one step replica symmetry breaking solution in the replica language. It was also highlighted that for the Weibull
class, i.e. the extreme value statistics of bounded variables, and the Fréchet class for power-law decaying variables, the replica method cannot be used without substantial modifications.

In order to shed light on how the extreme value statistics of disorder affects the glassy properties of the system in the low temperature phase we consider two versions of the Random Energy Model. We choose for the equilibrium distribution of configurations in the low temperature phase, i.e. the Boltzmann factors (20) $B_j = e^{-\beta\epsilon_j}$, both an exponential and a uniform distribution. We find a dynamical transition in the first case but not in the second. These results are corroborated by a study of the equilibrium autocorrelation function in both cases.

### A. Exponentially distributed Boltzmann factors

We recall that the dynamics of the Random Energy Model can be described through the following Master equation:

$$ \frac{dP_i(t)}{dt} = \sum_{j \neq i} W_{ij} P_j(t) - \sum_{j \neq i} W_{ji} P_i(t), \quad (34) $$

where the transition rates for level-dependent barriers are $W_{ji} = N^{-2p} B_j^{1-q} B_i^{-q}$, with $p \ge 0$ and $q \in (0, 1)$.

We define as a version of the Random Energy Model a system whose dynamics is still described by Eq. (34), with the difference that the statistics of the energy barriers in the low temperature phase has changed. Instead of a Lévy distribution we consider a probability distribution function of the Boltzmann factors of the form:

$$ p(B_i) = p(e^{-\beta\epsilon_i}) = N\theta(B_i)e^{-NB_i}, \quad (35) $$

where the scaling in $N$ has been chosen a posteriori to find a well defined thermodynamic limit of the equilibrium autocorrelation function (see section IV C 1).

The large deviation function $\psi_K(s)$ is again the largest solution of the eigenvalue equation

$$ \sum_i \frac{e^{-s} B_i^{1-2q}}{\zeta B_i^{-q} + (e^{-s}-1)B_i^{1-2q} + \psi_K(s)} = 1. \quad (36) $$

We solved numerically the equation (36) and averaged the solution over $10^7$ realizations of the disorder. The results are plotted in fig. 6.

As in the directed Trap Model and the Random Energy Model we observe a dynamical transition. For $s < 0$ the large deviation function appears to be constant in the system size. For $s > 0$ the curve approaches the horizontal axis for $N \gg 1$.

<figure><img src="image_7.png" /><figcaption>Figure 6: Numerical evaluation of the average large deviation function $\overline{\psi_K(s)}$ for exponentially distributed Boltzmann factors, with $p = q = 0.9$. Inset: numerical evaluation of $\overline{\psi_K(s)}$ for uniformly distributed Boltzmann factors.</figcaption></figure>

Thus even if the moments of the Boltzmann factors (or the mean trapping times) are well defined, the system presents dynamical heterogeneities and we expect an anomalous slowing down of the dynamics.

### B. Uniformly distributed waiting times

Let's consider uniformly distributed Boltzmann factors:

$$ p(B_i) = N\theta\left(B_i - \frac{a}{N}\right) \theta\left(\frac{a+1}{N} - B_i\right), \quad (37) $$

where $a$ is a finite constant with $0 \notin [a, a+1]$. If we allow the $B_i$'s take the value $B = 0$, after some time the system will be frozen in a single energy level.

Evaluating numerically the large deviation function of this system we find a complete independence of the large deviation function from the system size for all $s$ (see the inset in fig. 6), and no discontinuity of the large deviation function. Since we do not observe any dynamic transition in the space of histories, we expect this system to be dynamically homogeneous.

### C. Correlation function

To compute the equilibrium autocorrelation function of the Random Energy Model with exponentially and uniformly distributed Boltzmann factors we recast the Koper and Hilhorst's study [14] of the relaxation properties of the Random Energy Model.
Starting from the Master Equation (34) with transition rates (20), the authors proved that the autocorrelation function $C(t)$ of any physical quantity that does not couple to the energies of the REM can be written as a complex integral of the form:

$$C(t) = - \oint \frac{ds}{2\pi i} \frac{e^{-st\zeta} \sum_{j=1}^{N} \sum_{k \neq j}^{N} \frac{V_j B_j V_k B_k}{(s-V_j\zeta)(s-V_k\zeta)}}{sZ \sum_{l=1}^{N} \frac{V_l B_l}{s-V_l\zeta}} \quad (38)$$

with $Z = \sum_{i=1}^{N} B_i$ and where the line integral is taken over a contour oriented counterclockwise and encloses all the poles of the integrand on the positive real axis. If we average this expression, as in [14], over a Lévy distribution we find the asymptotic behaviors (22) which characterize the glassiness of the Random Energy Model.

Here, to corroborate the presence or absence of the dynamical transition obtained previously, we average Eq. (38) over an exponential distribution and then over a uniform distribution.

### 1. Exponential distribution

We consider exponentially distributed Boltzmann factors

$$p(B) = N\theta(B)e^{-NB}, \quad (39)$$

and transition rates of the form

$$W_{ij} = N^{-2q} B_i^{1-q} B_j^{-q}, \quad (40)$$

where the $N$ dependence is chosen to give a well defined thermodynamic limit of the correlation function. After some manipulations which details are given in Appendix A, one finds

$$\overline{C}(t) = t^{\frac{4}{q+1}} \int_{0}^{\infty} dx \int_{0}^{\infty} dy \, xy e^{-t^{\frac{1}{q+1}} \tilde{g}(x,y)}, \quad (41)$$

where $\tilde{g}(x,y) = x + y + \Gamma(2-q)/(x^q + y^q)$. One easily finds that $\tilde{g}(x,y)$ has a single minimum in $\mathbb{R}^+ \times \mathbb{R}^+$ for $x^* = y^* = (q\Gamma(2-q)/4)^{1/(q+1)}$ and $\partial_x^2 g(x^*, y^*) = \partial_y^2 g(x^*, y^*) = x^{*-1}$. Thus the large time behavior $\overline{C}(t)$ is given by a saddle point calculation yielding

$$\overline{C}(t) \sim 2\pi t^{\frac{3}{q+1}} x^{*3} e^{-t^{\frac{1}{q+1}} \tilde{g}(x^*, y^*)} \quad \text{for } t \gg 1 \quad (42)$$

which shows, given that $q > 0$, that the relaxation is indeed slower than exponential.

### 2. Uniform distribution

We now consider uniformly distributed Boltzmann factors:

$$p(B) = N\theta(B - \frac{a}{N})\theta(\frac{a+1}{N} - B), \quad (43)$$

where $a$ is a constant. After some manipulations left in the Appendix A, one has

$$\overline{C}(t) = \int_{a}^{a+1} dx \int_{a}^{a+1} dy \frac{xy}{a^2} e^{-t \frac{a^{1-q}}{x^q + y^q}}. \quad (44)$$

The minimum of the integrand is reached for $x = y = a$ and the maximum for $x = y = a + 1$, thus the integral is bounded from above and below by an exponential and we conclude that:

$$\overline{C}(t) \sim e^{-tc(a)}, \quad (45)$$

where $c(a)$ is a positive function of $a$. Eq. (45) is the asymptotic behavior of the correlation function of a non-glassy system.

## V. CONCLUSION

The central statement of our work is that when a system presents glassy features (anomalous slowing down of the dynamics, non exponential decay of the correlation functions, anomalous dependence of viscosity on temperature and so on) it appears heterogeneous in time and energy landscape, and a dynamical transition in the space of time realizations the system can follow in configuration space has to be found.

Here we validated this statement for two well known disordered models of glasses, the Random Energy Model and the Directed Trap Model, detecting in both a dynamical phase transition. In the Random Energy Model the influence of the statistics of disorder over the space-time heterogeneity of the system was analyzed. Indeed we found a dynamical transition for a REM with exponentially distributed Boltzmann factors and showed that we are in presence of a glassy dynamics by computing the equilibrium autocorrelation function of this system which is shown to decay anomalously, i.e. to be a stretched exponential in time. Furthermore we considered a uniform distribution of Boltzmann factors finding no dynamical transition and showing, as expected, that the equilibrium relaxation is exponential in time and no slowing down of the dynamics is present.

## ACKNOWLEDGMENTS

We thank C. Monthus for several useful discussions. This work was supported by the French Ministry of Education through Grant No. ANR-05-JCJC-44482.
# Appendix A: COMPUTATION OF THE CORRELATION FUNCTIONS FOR THE REM

In this appendix, we present some details concerning the computation of the correlation function for the kinetic REM.

## a. General framework

We recall the Master Equation (34) with transition rates (20):

$$ \frac{dP_i(t)}{dt} = V_i B_i \sum_{j=1}^{N} V_j P_j(t) - V_i P_i(t) \zeta. \quad (A1) $$

The probability to find the system at level i at time t if it was in j at time t = 0, i.e. the Green function, can be written in terms of the eigenvalues and eigenfunctions of the Master Equation. We call $\phi_i^\lambda$ the ith component of the right-hand eigenfunction with eigenvalue -λ:

$$ \frac{d\phi_i^\lambda}{dt} = -\lambda \phi_i^\lambda. \quad (A2) $$

Using (A1) one can write:

$$ \phi_i^\lambda = \frac{V_i B_i}{V_i \zeta - \lambda} \sum_{j=1}^{N} V_j \phi_j^\lambda. \quad (A3) $$

The detailed balance condition implies that the left-hand eigenfunction is $\psi_i^\lambda = \phi_i^\lambda / B_i$, so that the Green function takes the form:

$$ G_{ij}(t) = \sum_{l=1}^{N} e^{-\lambda_l t} \frac{\phi_i^{\lambda_l} \psi_j^{\lambda_l}}{\sum_{k=1}^{N} \phi_k^{\lambda_l} \psi_k^{\lambda_l}} \\ = \sum_{l=1}^{N} e^{-\lambda_l t} \frac{\frac{V_i B_i}{V_i \zeta - \lambda_l} \frac{V_j}{V_j \zeta - \lambda_l}}{\sum_{k=1}^{N} \frac{V_k^2 B_k}{(V_k \zeta - \lambda_l)^2}}. \quad (A4) $$

The equilibrium autocorrelation function expressed in terms of these Green functions is

$$ C(t) \equiv \sum_{j=1}^{N} \{G_{jj}(t) - G_{jj}(\infty)\} P_j^{\text{eq}}, \quad (A5) $$

where $P_j^{\text{eq}} = \lim_{t \to \infty} P_j(t) \propto B_j$.

To get rid of the dependence of the autocorrelation function $C(t)$ in Eq. (A5) on the unknown eigenvalues $\lambda_l$ we consider its Laplace transform $\tilde{C}(s)$, which using (A4) can be written as:

$$ \tilde{C}(s) = \frac{1}{Z} \sum_{l=2}^{N} \frac{H(\lambda_l)}{(s + \lambda_l)F'(\lambda_l)} \quad (A6) $$

with $Z = \sum_{i=1}^{N} B_i$ and where

$$ H(z) = \sum_{j=1}^{N} \frac{V_j^2 B_j^2}{(V_j \zeta - z)^2}, \quad F(z) = \sum_{j=1}^{N} \frac{V_j^2 B_j}{V_j \zeta - z} - 1. $$

We introduce $g(z) = H(z)/((s+z)F(z))$, which has poles in $z = -s$, in $z = V_j\zeta$ and in the zeroes of $F(z)$. Looking at the eigenvalue equation (A3) we see that the sum of the residues of $g(z)$ in the zeroes of $F(z)$ will be exactly $\tilde{C}(s)$. Thus, since the contour integral of $g(z)$ along a circle centered in the origin with radius $R$ vanishes as $R \to \infty$, one finds:

$$ \tilde{C}(s) = -\frac{1}{s} \left\{ \frac{\sum_{j=1}^{N} \sum_{k \neq j}^{N} \frac{V_j B_j V_k B_k}{(V_j \zeta + s)(V_k \zeta + s)}}{Z \sum_{l=1}^{N} \frac{V_l B_l}{V_l \zeta + s}} + \frac{1}{Z^2} \sum_{j=1}^{N} B_j^2 - 1 \right\}, \quad (A7) $$

from which we see that $\tilde{C}(s)$ has poles only on the negative real axis ($s = 0$ is not a pole). Applying an Inverse Laplace transform and changing s into $-s\zeta$ we obtain the expression of the autocorrelation function $C(t)$ given in the text in Eq. (38).

## b. Exponential distribution

We consider exponentially distributed Boltzmann factors

$$ p(B) = N \theta(B) e^{-NB}, \quad (A8) $$
Using for the denominator in Eq. (38) the integral representation $\alpha^{-1} = \int_{0}^{\infty} d\lambda e^{-\alpha\lambda}$, we obtain from

Eq. (38) the expression for the averaged autocorrelation function:

$$ \begin{aligned} \bar{C}(t) = \lim_{N \to \infty} N^4 \int_0^\infty dB_i \int_0^\infty dB_j \int_0^\infty d\lambda \int_0^\infty d\mu \oint \frac{ds}{2\pi i} \frac{B_i B_j}{s(B_i^q - 1)(sB_j^q - 1)} e^{-N(B_i+B_j)} \\ \times e^{-\lambda(B_i+B_j)-\mu(\frac{B_i}{sB_i^q-1}+\frac{B_j}{sB_j^q-1})-stN^{-2q}(B_i^{1-q}+B_j^{1-q})} \\ \times I_N^{N-2} \end{aligned} \quad (A9) $$

where the integral $I_N$ is given by

$$ I_N = N \int_0^\infty dB e^{-NB-\lambda B-\mu \frac{B}{sB^q-1}-stN^{-2q}B^{1-q}} $$

Rescaling the variable $B = B'/N$ and integrating by parts we have, to leading order in $N$:

$$ I_N \approx 1 - \frac{1}{N} \int_0^\infty dB e^{-B} (\lambda - \mu + stN^{-q}(1-q)B^{-q}) \quad (A10) $$

Thus using this large $N$ expansion (A10) in the formula (A9), and integrating over $\lambda$ and $\mu$ we get:

$$ \bar{C}(t) = \lim_{N \to \infty} N^4 \int_0^\infty dB_i \int_0^\infty dB_j \oint \frac{ds}{2\pi i} \frac{B_i B_j}{s(B_i + B_j + 1)} \frac{e^{-N(B_i + B_j) - stN^{-2q}(B_i^{1-q} + B_j^{1-q} + N^q \Gamma(2-q))}}{s(B_i^q B_j^q + B_i B_j^q) - B_i - B_j - (sB_i^q - 1)(sB_j^q - 1)} $$

Changing variables according to $NB_i = x$, $NB_j = y$ and keeping only the leading orders in $N$ we have:

text (41).

$$ \begin{aligned} \bar{C}(t) = \lim_{N \to \infty} N^4 \int_0^\infty dx \int_0^\infty dy \oint \frac{ds}{2\pi i} xy \\ \times \frac{e^{-(x+y)-stN^{-q}\Gamma(2-q)}}{s(sN^{-q}(x^q + y^q) - 1)}, \end{aligned} \quad (A11) $$

c. Uniform distribution

We now consider uniformly distributed Boltzmann factors:

which has one pole in $s = N^q/(x^q + y^q)$. Complex integration thus gives:

$$ \bar{C}(t) = \int_0^\infty dx \int_0^\infty dy \, xy e^{-(x+y)} e^{-t \frac{\Gamma(2-q)}{x^q+y^q}} \quad (A12) $$

$$ p(B) = N\theta(B - \frac{a}{N})\theta(\frac{a+1}{N} - B), \quad (A13) $$

Finally we change variables $x = x't^{1/(q+1)}$, $y = y't^{1/(q+1)}$ and find the expression given in the

where $a$ is a constant. The average of the correlation function (38) takes the form:

$$ \begin{aligned} \bar{C}(t) = \lim_{N \to \infty} N^4 \int_{\frac{a}{N}}^{\frac{a+1}{N}} dB_i \int_{\frac{a}{N}}^{\frac{a+1}{N}} dB_j \int_0^\infty d\lambda \int_0^\infty d\mu \oint \frac{ds}{2\pi i} \frac{B_i B_j}{s(B_i^q - 1)(sB_j^q - 1)} \\ \times e^{-\lambda(B_i+B_j)-\mu(\frac{B_i}{sB_i^q-1}+\frac{B_j}{sB_j^q-1})-stN^{-2p}(B_i^{1-q}+B_j^{1-q})} J_N^{N-2} \end{aligned} \quad (A14) $$

where

$$ J_N = \int_{\frac{a}{N}}^{\frac{a+1}{N}} dB N e^{-\lambda B - \mu \frac{B}{sB^q-1} - stN^{-2p} B^{1-q}} \quad (A15) $$

Since the interval of integration is very small when $N \gg 1$ we can consider the integrand as a constant
over this interval and make the approximation:

$$J_N \approx e^{-\lambda a N^{-1} + \mu a N^{-1} - st N^{-2p-(1-q)a^{1-q}}} \quad (A16)$$

Thus after integration over $\lambda$ and $\mu$ and changes of variables $x = NB_i$ and $y = NB_j$, the correlation function becomes:

The integrand has one pole in $s = N^q/(x^q + y^q)$ so after complex integration and choosing $p = q$ we obtain the formula given in the text (44).

$$\begin{aligned} \bar{C}(t) = & \lim_{N \to \infty} \int_a^{a+1} dx \int_a^{a+1} dy \oint \frac{ds}{2\pi i} xy \\ & \times \frac{e^{-stN^{-2p+q}a^{1-q}}}{sa^2(sN^{-q}(x^q + y^q) - 1)} \end{aligned} \quad (A17)$$

[1] J.P. Bouchaud, L. Cugliandolo, J. Kurchan, M. Mézard, Out of Equilibrium dynamics in Spin-Glasses and other Glassy systems, cond-mat/9702070, (1997); L. Cugliandolo, Lecture notes in Slow Relaxation and non equilibrium dynamics in condensed matter, Les Houches Session 77 (2002).

[2] L. Santen, W. Krauth, Absence of Thermodynamic Phase Transition in a Model Glass Former, Nature 405, 550 (2000).

[3] J.P. Garrahan and D. Chandler, Geometrical explanation and scaling of dynamical heterogeneities in glass forming systems, Phys. Rev. Lett. 89, 035704 (2002); M. Merolle, J.P. Garrahan and D. Chandler, Space-time Thermodynamics of the Glass Transition, Proc. Natl. Acad. Sci. USA 102, 10837 (2005); R.L. Jack, J.P. Garrahan and D. Chandler, Space-time thermodynamics and subsystem observables in a kinetically constrained model of glassy systems, J. Chem. Phys. 125, 184509 (2006).

[4] D. Ruelle, Thermodynamics formalism Addison-Wesley, Reading, (1978); J.P. Eckmanns D. Ruelle, Ergodic theory of chaos and strange attractors, Rev. Mod. Phys. 57, 617 (1985).

[5] J.P. Garrahan, R.L. Jack, V. Lecomte, E. Pitard, K. van Duijvendijk, F. van Wijland, Dynamic first-order phase transition in kinetically constrained models of glasses, Phys. Rev. Lett. 98, 195702 (2007).

[6] V. Lecomte, C. Appert-Rolland and F. van Wijland, Thermodynamic formalism for systems with Markov dynamics, J. Stat. Phys. 127, 1 (2007).

[7] J.P. Bouchaud, A. Comtet and C. Monthus, On a dynamical model of glasses, J. Physique I 5, 1521 (1995); C. Monthus and J.P. Bouchaud, Models of traps and glass phenomenology, J. Phys. A 29, 3847 (1996).

[8] G. Ben Arous and J. Černý, Dynamics of Trap Models, École d'été de Physique des Houches, Session LXXXIII, 331-394 (2006).

[9] J.P. Bouchaud, Weak ergodicity breaking and aging in disordered systems, J. Phys. I (France) 2, 1705 (1992); J.P. Bouchaud and D. Dean, Aging on Parisi's tree, J. Phys. I (France) 5, 265 (1995); G. Ben Arous, A. Bovier, V. Gayrard, Aging in the Random Energy Model, Phys. Rev. Lett. 88, 087201 (2002).

[10] F. Ritort, Universal Dependence of the fluctuation-dissipation ratio on the transition rates in trap models, J. Phys. A 36, 10791 (2003).

[11] J.P. Bouchaud, M. Mézard, Universality classes for extreme value statistics, J. Phys. A: Math. Gen. 30, 7997 (1997).

[12] G. Biroli, J.P. Bouchaud, M. Potters, Extreme value problems in Random Matrix Theory and other disordered systems, J. Stat. Mech. 07, 07019 (2007).

[13] B. Derrida, Random-energy model: An exactly solvable model of disordered systems, Phys. Rev. B 24, 2613 (1981).

[14] G.J.M. Koper and H.J. Hilhorst, Nonexponential relaxation in the Random Energy Model, Physica A 160, 1 (1989).