format: text

text:
# Energy Barriers and Activated Dynamics in a Supercooled Lennard-Jones Liquid

## I Introduction.

More than thirty years ago, Goldstein Goldstein (1969) proposed to view a glass-forming system as a point moving in the high-dimensional landscape of the potential energy $`V(x)`$. In this framework he suggested to focus onto the local minima of the potential energy landscape (PEL), where the system is supposed to be trapped at low enough temperatures. Via occasional transitions to neighboring minima the system finally relaxes. Owing to the separation of time scales, one is able to describe many features of glass formers by properties of only the minima. Stillinger and Weber Stillinger and Weber (1982) formulated this idea in the language of statistical thermodynamics using the concept of basins. A basin of a given minimum is defined as the set of configurations that reach this minimum via their steepest descent path $`\dot{x}=F(x)`$. (We set $`x`$ and $`F(x)`$ as shorthands for the multidimensional vectors all particle positions and all forces, respectively.) The resulting tiling of configuration space into different basins allows one to write the free energy approximately as a function of static properties of minima, i.e. their energies and vibrational frequencies Büchner and Heuer (1999); Sciortino et al. (2000). Knowledge of the thermodynamics is in general not sufficient to predict dynamical properties like diffusion constants or relaxation times. However, experimental Richert and Angell (1998) as well as simulated Saika-Voivod et al. (2001); Sastry (2001); Scala et al. (2000) data seem to indicate that there exists a strong connection between dynamics and thermodynamics via the Adam-Gibbs relation Adam and Gibbs (1965).

Our goal is to reach a quantitative understanding of the slowing down of dynamics, as expressed by the bulk long-time diffusion constant $`D(T)`$. Mode coupling theory (MCT) Götze and Sjogren (1992) predicts a power-law behavior of the form $`1/D(T)(TT_c)^\gamma `$ above the MCT critical temperature $`T_c`$. Since $`T_c`$ is found to be higher than the glass transition temperature $`T_g`$, the MCT divergence of $`1/D(T)`$ at $`T_c`$ is not observed in practice. The common explanation for this shortcoming of MCT is that the theory neglects ‘activated processes’, or ’hopping’, which is supposed to come into play around and below $`T_c`$. Indeed it was proven by Schrøder et al. Schrøder et al. (2000) that in the vicinity of $`T_c`$ the time scale of fast local dynamics around single minima becomes well separated from the time scale of interbasin transitions. Above $`T_c`$ the common picture suggests that the dynamics are ’entropy-driven’, i.e. characterized by the search for escape directions Donati et al. (2000); Nave et al. (2000, 2002), since saddles lie far below the instantaneous potential energy of the system and thus represent no serious barriers Angelani et al. (2000a); Angelani et al. (2002). Also different observables like the average order of saddles Angelani et al. (2000a); Broderix et al. (2000) and the number of diffusion-like normal modes Donati et al. (2000); Nave et al. (2000, 2002, 2001) seem to indicate that well above $`T_c`$ the dynamics is not governed by activated transitions between adjacent minima. However, in the multi-dimensional space of particle coordinates, it is not obvious how to distinguish thermally activated from entropically limited dynamics. One possible way to do this will be discussed below. For the time being, we use the term hopping only in a formal sense, meaning that the trajectory of the system is mapped onto a sequence of jumps among minima.

Looking for a quantitative link of bulk diffusion $`D(T)`$ to PEL properties, we recently investigated hopping dynamics on the PEL in greater detail Doliwa and Heuer (2002). A priori, temporal and spatial aspects of hopping events have to be considered, the former in the shape of the waiting time distribution (WTD) of jumps, the latter by the jump lengths and directions, and correlations thereof. We found that strong backward correlations of jumps arise from the organization of minima into superstructures, which, following Stillinger (1995), we call metabasins (MB). It had already been known from previous work Büchner and Heuer (2000) that structural relaxation corresponds to jumps among MBs rather than single basins. MBs were identified with the help of a straightforward algorithm such that close-by minima between which the system performs several back- and forth jumps are identified as a single MB. Then, indeed, hopping among MBs was found to be close to a random walk with a distribution of MB waiting times. Motivated by this fact, we expressed $`D(T)`$ in the simple form,

$$D(T)=\frac{a^2}{6N\tau (T)},$$
(1)
with the mean waiting time $`\tau (T)`$ and the effective jump length $`a(T)`$. With this ansatz, we anticipated that waiting times would carry the major part of the temperature dependence. Indeed, $`a(T)`$ turned out to be constant for $`T<2T_c`$, which is why we dropped the argument of $`a(T)`$. Eq. 1 constitutes an important step towards the understanding of diffusion in supercooled liquids: it suffices to look for the physics behind MB waiting times, spatial details of hopping being expressed by a single constant.

A simple model for hopping dynamics has been discussed by Bouchaud and coworkers Monthus and Bouchaud (1996). They consider the relaxation from traps of depths $`E`$ with distribution $`\rho (E)`$, and escape rates $`\gamma (E,T)=\gamma _0\mathrm{exp}(\beta E)`$. When $`\rho (E)\mathrm{exp}(E/T_x)`$, the WTDs assume power-law tails $`\psi (\tau )\tau ^{\alpha (T)}`$ with exponents $`\alpha (T)=1+T/T_x`$. The consequence is the divergence of the mean waiting time at temperature $`T_x`$. In our recent paper, we observed that the WTDs of a binary Lennard-Jones system are in conformance with this kind of power-law decay Doliwa and Heuer (2002). As a consequence of such slowly decaying WTDs, the mean value $`\tau (T)`$ was found to be dominated by the few, very long waiting times. In other words, the temperature dependence of $`D(T)`$ follows alone from the durations of trapping in the very stable MBs. These results were obtained for small binary Lennard-Jones mixtures of $`N=65`$ particles. For a macroscopic system, which, due to its dynamic heterogeneity Sillescu (1999), contains many slow and fast subsystems in parallel, this implies the dominance of slow regions in the temperature dependence of $`D(T)`$.

The logical continuation along this line of thinking is to relate MB lifetimes to the PEL topography. The most prominent characteristics of a MB is, of course, its energy $`ϵ_{\text{MB}}`$, which is defined as the lowest energy of all its constituent minima. It is then natural to introduce the mean MB lifetime $`\tau (ϵ_{\text{MB}};T)`$ at constant $`ϵ_{\text{MB}}`$. Knowledge of $`\tau (ϵ_{\text{MB}};T)`$, together with the population of MBs, $`p(ϵ_{\text{MB}};T)`$, is sufficient to calculate $`\tau (T)`$ and thus $`D(T)`$, as we will show now. We write

$$\tau (T)=dϵ_{\text{MB}}\tau (ϵ_{\text{MB}};T)\phi (ϵ_{\text{MB}};T),$$
(2)
where $`\phi (ϵ_{\text{MB}};T)`$ is the distribution of MBs visited at temperature $`T`$. We will see that this decomposition can be achieved by a detailed analysis of the hopping dynamics. Since $`p(ϵ_{\text{MB}};T)`$ denotes the probability that at a given time the system is in a MB with energy $`ϵ_{\text{MB}}`$, it is proportional to $`\phi (ϵ_{\text{MB}};T)`$ and the time $`\tau (ϵ_{\text{MB}};T)`$ the system remains in MBs of this energy. With the appropriate normalisation one gets

$$p(ϵ_{\text{MB}};T)=\frac{\tau (ϵ_{\text{MB}};T)}{\tau (T)}\phi (ϵ_{\text{MB}};T).$$
(3)
From Eqs. 1, 2 and 3, it immediately follows the representation

$$D(T)=\frac{a^2}{6N}\frac{1}{\tau (ϵ_{\text{MB}};T)}_T.$$
(4)
Here, $`\mathrm{}_T`$ denotes the canonical time average (w.r.t. $`p(ϵ_{\text{MB}};T)`$), while $`\mathrm{}`$ is the average over MBs. Hence,

$$\{\tau (ϵ_{\text{MB}};T),p(ϵ_{\text{MB}};T)\}\tau (T)D(T),$$
(5)
where the second implication has been established in our recent paper Doliwa and Heuer (2002). The population $`p(ϵ_{\text{MB}};T)`$ is related to the single-basin population $`p(ϵ;T)`$, a purely static quantity, which has been extensively discussed in the literature Sciortino et al. (2000); Sastry (2001); Büchner and Heuer (1999). It has turned out for Lennard-Jones mixtures that the number density $`G_{\text{eff}}(ϵ)`$ of minima is approximately gaussian. Thus, the population of minima, $`p(ϵ;T)G_{\text{eff}}(ϵ)e^{\beta ϵ}`$, could be expressed by three parameters describing global PEL structure. In the present paper, we focus on $`\tau (ϵ_{\text{MB}};T)`$, our goal being to deduce it from PEL structure. If this succeeds, we have established the following connection,

$$\text{local + global PEL structure}\text{long-time dynamics},$$

which, in our opinion, pushes the understanding of diffusion in supercooled liquids a step further.

We proceed as follows. We first compute MB lifetimes from ordinary simulation, and later compare them to the prediction from PEL structure. First, we characterize the relaxation from four single, randomly selected MBs. By an exhaustive sampling of these MBs, we will be able to get some first insights into MB topology. Second, many MBs of fixed energy are considered and their lifetimes $`\tau (ϵ_{\text{MB}};T)`$ are calculated. Third, we relate MB lifetimes to PEL structure, by quantifying the MB depths, or effective barriers, which determine the temperature dependence of $`\tau (ϵ_{\text{MB}};T)`$. The physical scenario which will emerge from the results of this paper implies that MBs can be regarded as traps, surrounded by high barriers. It turned out from exhaustive explorations of PEL connectivity Doye and Wales (2001) that due to the high dimensionality of configuration space the number of escape paths from every minimum is enormous. Thus, one may anticipate that the effective barrier to leave a specific MB results as a complex superposition of individual escape paths. Therefore, enormous numerical effort is required to quantify their multitude for many different MBs.

Note that the whole analysis will be carried out in the spirit of activated barrier crossing. The extent to which this is present in supercooled liquids is quite disputed in literature. However, we will show that for temperatures in the landscape-influenced regime below $`2T_c`$, the apparent activation energy

$$E_{\text{app}}(ϵ_{\text{MB}};T)=\frac{\mathrm{d}}{\mathrm{d}\beta }\mathrm{ln}\tau (ϵ_{\text{MB}};T),$$
(6)
can indeed be identified with PEL barriers much larger than $`k_BT`$ which the system encounters when leaving a MB. Thus, together with Eq. 4, we will find that the activated escape out of deep traps is the physical mechanism behind diffusion.

To our knowledge, such a connection between dynamics and PEL barriers has never been established for a fragile glass former. In contrast, for $`\text{SiO}_2`$, the apparent activation energy of diffusion below $`T_c`$ could be related to the simple breakage of $`\text{Si}\text{O}`$ bonds Mcmillan et al. (1994); Horbach and Kob (1999).

The organization of the paper is as follows. In section II, we provide the details of our simulation, and describe the interval bisection method to identify MBs. Section III deals with the computation of apparent activation energies from relaxation dynamics. The corresponding energy barriers will be addressed in section V, after introducing our technique for finding transition states (section IV). In section VI, we independently demonstrate that barriers and associated reaction paths indeed govern relaxation. Finally, we discussion further aspects of our results in section VII and conclude in section VIII.

## II Simulation Details.

### II.1 General.

In the present work, we investigate a binary mixture of Lennard-Jones particles (BMLJ), as recently treated by two groups Broderix et al. (2000); Hernandez-Rojas and Wales (2001); see also Kob and Andersen (1995). It is characterized by the interaction potentials

$$V_{\alpha \beta }(r)=4ϵ_{\alpha \beta }[(\sigma _{\alpha \beta }/r)^{12}(\sigma _{\alpha \beta }/r)^6]$$

with the parameter set $`N=N_A+N_B=52+13=65`$, $`\sigma _{AB}=0.8\sigma _{AA}`$, $`\sigma _{BB}=0.88\sigma _{AA}`$, $`ϵ_{AB}=1.5ϵ_{AA}`$, $`ϵ_{BB}=0.5ϵ_{AA}`$, $`r_c=1.8`$. Linear functions were added to the potentials to ensure continuous forces and energies at the cutoff $`r_c`$. These modifications of the original potential by Kob and Andersen Kob and Andersen (1995) are necessary for the simulation of small systems. We use Langevin molecular dynamics simulations (MD) with fixed step size, $`\lambda ^2=0.015^2=2k_\mathrm{B}T\mathrm{\Delta }t/m\zeta `$, equal particle masses $`m`$, friction constant $`\zeta `$ set to unity, and periodic boundary conditions. Units of length, mass, energy, and time are $`\sigma _{AA}`$, $`m`$, $`ϵ_{AA}`$, and $`m\zeta \lambda ^2/2ϵ_{AA}`$, respectively. However, we will omit these units, for convenience. The mode-coupling temperature is $`T_c=0.45\pm 0.01`$ in this model system (compare Kob and Andersen (1995)). For the analysis of dynamics from the PEL perspective it is essential to use small systems, as has been stressed in the literature Büchner and Heuer (2000); Keyes and Chowdhary (2001); Grigera et al. (2002). On the other hand, naturally, the system should not be too small in order to avoid major finite-size effects. For the BMLJ, $`N60`$ turns out to be a very good compromise Broderix et al. (2000); Hernandez-Rojas and Wales (2001); Büchner and Heuer (2000), whereas $`N40`$ already causes large finite-size effects Büchner and Heuer (1999). Here we choose $`N=65`$, since the BMLJ60 system has a stronger tendency to be trapped in crystalline configurations. We stress here that the results obtained for the BMLJ65 system show no finite-size related artifacts. For example, $`D(T)`$ of the BMLJ1000 is identical to $`D(T)`$ of the BMLJ65 above $`T_c`$, see Fig. 4. In the temperature range studied, we found that the behavior of a BMLJ130 system largely resembles that of two independent copies of a BMLJ65. Thus, the generalization of the present work to larger systems should not bear any pitfalls.

Interestingly, the BMLJ65 relaxation becomes Arrhenius-like for low temperatures. Since we have no equilibrium runs of our BMLJ1000 below $`T_c`$ it is unclear whether the Arrhenius behavior of the BMLJ65 down to $`T=0.4`$ is a finite-size effect. A possible explanation is that the lower end of the PEL is reached (located at $`ϵ_{\text{min}}302`$, see Fig. 3), preventing further increase of barriers. In turn, this may be related to the fact that cooperative regions cannot grow any further in the BMLJ65, i.e. structural optimization \- which happens upon cooling - finally comes to an end. Since we work above $`T_c`$, our key results are not affected by this argument.

### II.2 Interval bisection.

By regularly quenching the MD trajectory $`x(t)`$ to the bottom of the basins visited at time $`t`$, as proposed by Stillinger and Weber, we obtain a discontinous trajectory $`\xi (t)`$. A problem from the standpoint of simulations is to resolve the elementary hopping events. Since computer time prohibits to calculate the minimum $`\xi (t)`$ for every time step $`t`$, we normally find ourselves in the situation of having equidistant quenched configurations $`\xi (t_i)`$, $`t_i=i\mathrm{\Delta }t`$, with, say, $`\mathrm{\Delta }t=10^5`$ MD steps. If the same minimum is found for times $`t_i`$ and $`t_j`$, we need not care about transitions in the meantime, because no relaxation has occurred there. If, in contrast, $`\xi (t_i)\xi (t_{i+1})`$, we must not expect $`\xi (t_{i+1})`$ to be the direct successor of $`\xi (t_i)`$, since many other minima could have been visited between $`t_i`$ and $`t_{i+1}`$. Therefore, further minimizations in this time interval are necessary. For reasons of efficiency, we apply a straightforward interval bisection method, which locates transitions to an accuracy of 1 MD step: provided $`\xi (t_{\text{start}}^{(0)})\xi (t_{\text{start}}^{(1)})`$, (a) set $`t^{(0)}t_{\text{start}}^{(0)}`$, $`t^{(1)}t_{\text{start}}^{(1)}`$, (b) reconstruct the trajectory $`x(t)`$ at time $`t^{(2)}=(t^{(0)}+t^{(1)})/2`$, (c) calculate $`\xi (t^{(2)})`$, (d) if $`\xi (t^{(2)})=\xi (t^{(0)})`$, set $`t^{(0)}t^{(2)}`$, else set $`t^{(1)}t^{(2)}`$, (e) repeat (b)-(d) until $`t^{(1)}t^{(0)}=1`$ MD step. Repeated application of the interval bisection to a simulation run $`x(t)`$ finally gives all relevant transitions. Note that the determination of all transitions including the numerous recrossings of basin borders would require minimization for every MD step! The interval bisection method thus may oversee back- and forth motions between minima which, in any event, are irrelevant for relaxation. Although computationally demanding, the above method has proved most efficient for resolving the relevant details of hopping on the PEL and is predestined for the construction of metabasins (see below).

## III Activation Energies from Metabasin Lifetimes.

### III.1 Metabasin lifetime construction.

As said above, stable configurations in the supercooled liquid are rarely due to single minima on the PEL, but mostly correspond to groups of strongly correlated minima. While the system is trapped in such a MB for a long time, a small number of minima is visited over and over again. This is well reflected by the time series of potential energies, $`ϵ(t)=V(\xi (t))`$ Büchner and Heuer (2000); Doliwa and Heuer (2002). In this section, we will dwell on the computation of mean MB lifetimes, (i) for single, selected MBs, and (ii) averaged over MBs of a given energy $`ϵ_{\text{MB}}`$, thus yielding $`\tau (ϵ_{\text{MB}};T)`$. The individual MBs of (i) correspond to long-lived MBs and thus represent typical MBs which govern the temperature dependence of $`\tau (T)`$.

For the grouping of minima into MBs and the resulting determination of their lifetimes from a regular simulation run, we use the following, straightforward algorithm Büchner and Heuer (2000). (a) determine the regions $`[t_i^{},t_i^{}]`$ where $`t_i^{}`$ is the time of the first and $`t_i^{}`$ the time of the last occurrence of minimum $`\xi (t_i)`$, (b) any two regions overlapping by less than $`\tau _{\text{mol}}`$ are cut so that $`[t_i^{},t_i^{}][t_j^{},t_j^{}]=\mathrm{}`$, where $`\tau _{\text{mol}}`$ is a small molecular time-scale, (c) any two regions overlapping by more than $`\tau _{\text{mol}}`$ are combined to $`[t_i^{},t_i^{}][t_j^{},t_j^{}]`$, (d) the lifetimes of MBs are defined by the regions after step (c), (e) the MB energy $`ϵ_{\text{MB}}`$ is defined as the lowest energy of minima visited during the MB lifetime.

A few comments on the procedure are in order. Time regions in (a) are determined by the interval bisection method which yields the time of transition from one minimum to another with an accuracy of one MD step. Step (b) is motivated by the observation that recrossings of a basin border during a transition are very probable. If we ignore this fact, i.e. set $`\tau _{\text{mol}}=0`$, step (c) would merge nearly all regions and we would end up with unphysically long MBs. Instead, we use $`\tau _{\text{mol}}=40`$, which is the basin equilibration time obtained from energy autocorrelation. Since the durations of transitions are of the same order of magnitude as $`\tau _{\text{mol}}`$, this is a sensible choice. However, we found that the results for MB lifetimes are not very susceptible to the precise value of $`\tau _{\text{mol}}`$. Step (c) itself is the realization of the MB concept. It is important to note that, different from Büchner and Heuer (2000), we will treat all MBs on the same footing here, no matter if they are short-lived or long-lived.

So far, the MB lifetime construction rests upon single trajectories, which only partially reflect the configuration space topology. In section V, the MB concept will be given a more precise, static definition, based on the return probability to the ground minimum.

### III.2 Activation Energies for Single MBs.

As noted above, the temperature dependence of $`D(T)`$ is dominated by the long-lived MBs. Generally, these are low-lying MBs, i.e. deep traps in the PEL. Since different MBs differ in their stability, a statistical treatment will be needed. As a first step, however, we restrict ourselves to the investigation of single MBs.

The relaxation times computed in this section do not stem from regular, linear simulation runs, but are obtained by artifically placing the system in a specific MB and waiting for its escape (’escape runs’). The above algorithm for the MB lifetime construction implicitly assumes that MBs finally have been left. In other words, the algorithm may not be used to determine the time where to stop the simulation due to successful escape. Fortunately, we can avoid running into this paradoxical situation by judging from an independent criterion whether an escape has been completed: if the distance of the instantaneous minimum to the starting position is greater than $`d_{\text{max}}=4`$, returning to the original basin can practically be excluded (see section V for a justification of $`d_{\text{max}}=4`$). Then, by applying the MB construction algorithm to the escape run, we obtain the lifetime of the MB.

We analyzed four low-lying ($`ϵ_{\text{MB}}<300`$), randomly selected MBs in greater detail. By repeated starts from the bottom of the MBs, we computed the mean lifetimes $`\tau _i(T)`$ as a function of temperature. From Fig. 1, we see that the relaxations from all MBs follow nicely an Arrhenius law. We note that, due to starting at minima, a short intra-basin equilibration time ($`\tau _{\text{mol}}=40`$, from energy autocorrelation) has been subtracted from the raw $`\tau _i(T)`$.

The fact that an Arrhenius form of $`\tau (ϵ_{\text{MB}};T)`$ is observed indicates that the barriers do not change any further upon lowering temperature. Put differently, MBs serve as traps surrounded by barriers with heights around $`E_{\text{app}}(i)=\mathrm{d}\mathrm{ln}\tau _i(T)/\mathrm{d}\beta `$. We will see in section V that this is indeed correct. Since $`E_{\text{app}}(i)/k_\mathrm{B}T_c>10`$, this implies a strongly activated dynamics near $`T_c`$.

### III.3 Activation Energies vs. MB energies.

As a further step, we analyze the mean relaxation time from MBs with the same energy, $`\tau (ϵ_{\text{MB}};T)`$; see Eq. 2. Clearly, the low $`ϵ_{\text{MB}}`$ are not populated at high temperatures so that regular simulation does not yield $`\tau (ϵ_{\text{MB}};T)`$ over a wide temperature range. We therefore artificially place the system in the desired MBs (in the lowest minima $`ϵ_{\text{MB}}`$ thereof) and measure the escape times as a function of temperature. Averaging over many different MBs, we obtain $`\tau (ϵ_{\text{MB}};T)`$.

Results are shown in Fig. 2 as a function of $`ϵ_{\text{MB}}`$. Below $`T=1`$, all relaxation times display Arrhenius behavior. Thus the apparent activation energies $`E_{\text{app}}(ϵ_{\text{MB}};T)`$ are temperature independent. In the following we will therefore omit the second argument. Thus, we can write

$$\tau (ϵ_{\text{MB}};T)=\tau _0(ϵ_{\text{MB}})e^{\beta E_{\text{app}}(ϵ_{\text{MB}})}.$$
(7)
As expected, the properties of MBs as expressed by $`E_{\text{app}}(ϵ_{\text{MB}})`$ depend on their ground state energy $`ϵ_{\text{MB}}`$.

We can interpret $`E_{\text{app}}(ϵ_{\text{MB}})`$ as the mean effective depth of MBs at $`ϵ_{\text{MB}}`$. Since the lower end of the energy landscape is reached at $`ϵ302`$ no deeper traps exist (compare Fig. 3, see also <sup>1</sup><sup>1</sup>1B. Doliwa and A. Heuer, in preparation). A simple statement for the depths of traps would follow if the rims of all traps were at the same level $`ϵ_{\text{th}}`$. The consequence would be $`E_{\text{app}}(ϵ_{\text{MB}})=ϵ_{\text{th}}ϵ_{\text{MB}}`$, for all $`ϵ_{\text{MB}}<ϵ_{\text{th}}`$. This simple scenario is ruled out by the data, see Fig. 2(b). Actually, a more complicated energy dependence of $`E_{\text{app}}(ϵ_{\text{MB}})`$ is expected from the very fact that the system - despite its small size - is not a completely cooperative unit, see the discussion in section VII.

The fact that we still observe Arrhenius-like relaxation in Fig. 2 indicates that the variation of trap depths at constant $`ϵ_{\text{MB}}`$ is not large, compare $`E_{\text{app}}(i)`$ from Fig. 1. Otherwise, $`E_{\text{app}}(ϵ_{\text{MB}})`$ would increase upon decreasing temperature, due to the more and more dominant, extremely deep traps. In contrast, trap depths at constant $`ϵ_{\text{MB}}`$ seem rather well defined by $`ϵ_{\text{MB}}`$, which suggests the existence of some underlying topological principle.

As seen from Fig. 2(c), the prefactor $`\tau _0(ϵ_{\text{MB}})`$ has no strong dependence on $`ϵ_{\text{MB}}`$. From high energies, it decreases at most an order of magnitude and seems to level off below $`ϵ_{\text{MB}}=297`$. Hence, for the range of energies that dominate $`\tau (T)`$ at low temperatures, it can be considered constant within error bars. In contrast to $`E_{\text{app}}(ϵ_{\text{MB}})`$, we will not be able to deduce $`\tau _0(ϵ_{\text{MB}})`$ from PEL structure. Its weak variation is therefore quite fortunate.

We will now analyze the second factor of the integrand in Eq. 2, $`\phi (ϵ_{\text{MB}};T)`$. It is shown in Fig. 3. Interestingly, the variation of $`\phi (ϵ_{\text{MB}};T)`$ is much weaker for low $`T`$ as the variation of $`p(ϵ_{\text{MB}};T)`$. From Eqs. 3 and 7, one concludes that the constancy of the distribution $`\phi (ϵ_{\text{MB}};T)`$ is equivalent to having $`E_{\text{app}}(ϵ_{\text{MB}})=ϵ_{\text{th}}ϵ_{\text{MB}}`$, with some constant $`ϵ_{\text{th}}`$. Since this simple behavior is not present, one must still have a residual temperature dependence of $`\phi (ϵ_{\text{MB}};T)`$.

Concerning $`p(ϵ_{\text{MB}};T)`$, the weak temperature dependence of its first moment for the three lowest temperatures is simply related to the probing of the lower end of the PEL. Actually, it turns out that $`p(ϵ_{\text{MB}};T)`$ is, within statistical error, identical to the corresponding distribution of minima $`p(ϵ;T)`$. One would expect this for high $`ϵ_{\text{MB}}`$, because no pronounced MBs are observed there. Considering a deep MB with many minima, this will equally effect no large difference between $`p(ϵ_{\text{MB}};T)`$ and $`p(ϵ;T)`$. The reason is that the group of minima near $`ϵ_{\text{MB}}`$ carry the largest part of the population. Since they are close to $`ϵ_{\text{MB}}`$, transferring their weight to $`ϵ_{\text{MB}}`$ when computing $`p(ϵ_{\text{MB}};T)`$ has little effect.

As a consistency check, we use the data from Fig. 2 and 3 to reproduce $`\tau (T)`$ indirectly via Eq. 2 (denoted $`\tau (T)_{\text{ind}}`$). The match with $`\tau (T)`$ is not completely trivial since the data for $`\tau (T)`$ and $`\phi (ϵ_{\text{MB}};T)`$ were gathered from a linear simulation run, while $`\tau (ϵ_{\text{MB}};T)`$ results from selected MBs of certain $`ϵ_{\text{MB}}`$, where the system has been artificially placed. As shown in Fig. 4, the agreement of $`\tau (T)`$ and $`\tau (T)_{\text{ind}}`$ is good for $`T1`$ within the possible accuracy. Note that there is no free fit parameter between them. The deviation at $`T=2`$ can be explained by the fact that $`\tau (ϵ_{\text{MB}};T)`$, above $`T=1`$, and especially for the high $`ϵ_{\text{MB}}`$, departs from Arrhenius behavior, see Fig. 2(a).

So far, all barriers or trap depths have been derived indirectly, from the temperature dependence of waiting times. A link to the PEL structure is still lacking. For instance, the activation energies $`E_{\text{app}}(i)`$ of this section are expected to reflect the local topography of the selected MBs. Indeed, they can be identified from the barriers of escape paths, as will be demonstrated in section V.

First of all, the barriers between neighboring minima are of great interest. These are known once we have in hand the corresponding transition states.

## IV Non-Local Ridge Method for finding Transition States.

### IV.1 Description of the method.

We now describe how to determine transition states (TS) from the simulation, by what we call the (non-local) ridge method. The principle idea is that TSs are local minima of basin borders. They can be pictured as the lowest points of mountain ridges on the PEL. If the system crosses a basin border at time $`t`$, the steepest descent path starting from $`x(t)`$ should end up in a TS, see Ionova and Carter (1993). In practice, however, the descent will deviate from the ridge due to numerical error, finally ending up in the minimum $`\xi _0\xi (t)`$ or $`\xi _1\xi (t+)`$. As a way out, we let the system perform two descents in parallel, on either side of the basin border, as schematically depicted in Fig. 5. More specifically, if a transition happened at time $`t`$, interval bisection yields the configurations $`x(t)y_0`$ and $`x(t+1\text{ MD step})y_1`$. From these, by further interval bisection on the straight line between $`y_0`$ and $`y_1`$, the distance to the border may be further reduced if necessary, resulting in two configurations, again called $`y_0`$ and $`y_1`$. Close as they are, they still belong to different basins. If we now let descend $`y_0`$ and $`y_1`$ in parallel, they first move along the ridge towards the transition state until they finally bend off to their respective minima. This separation is clearly not wanted, so from time to time we reduce their distance by interval bisection. After a few iterations (descents+interval bisection) the vicinity of the transition state is reached in most cases. We then use a short minimization of the auxiliary potential $`\stackrel{~}{V}=\frac{1}{2}|F(x)|^2`$ followed by a few steps of Newton-Raphson type, which bring the search for the TS to a quick convergence. Besides a vanishing force, the resulting configuration $`\zeta `$ has a Hessian matrix with one negative eigenvalue. After small displacements along the corresponding eigenvector, one reaches the adjacent minima via steepest descent. This yields the reaction path (RP) $`\zeta (s)`$, where $`s`$ is a curvilinear parameter. We set $`\zeta (0)=\zeta `$, $`\zeta (s_0)=\xi _0`$, and $`\zeta (s_1)=\xi _1`$, where either $`s_0`$ or $`s_1`$ is negative.

It can happen, though, that no saddle between $`y_0`$ and $`y_1`$ is found, but that the interval bisection locates a third minimum. The basin border splits into two at this point, and no direct saddle between the initial and final minimum is available. Thus, we also have to split the descent along the basin border into two processes and then continue separately. If the two descents are successful without further bifurcations, we are finished and have the optimum reaction path which takes a detour via a third minimum. In such a situation, the RP is clearly not very useful. It has to be stressed that bifurcations are no artifacts of the ridge method, but a topological feature of some basin borders on the PEL. Fortunately, as a signature of strong anharmonicity, they are quite rare and happen to occur only in the high-energetic regions of the PEL. For the escapes from long-lived MBs, they are of no importance.

A similar algorithm is described in the literature Ionova and Carter (1993), which, instead of minimization and interval bisection, uses local maximization between $`y_0`$ and $`y_1`$ to prevent the configurations from moving apart. Although computationally less expensive, this method is not appropriate for our purpose. As an effect of the high dimensionality, the local shape of the PEL around $`y_0`$ and $`y_1`$ gives no direct clue to the membership to basins. When descending, one may thus loose the important property of $`y_0`$ belonging to the basin of $`\xi _0`$ and $`y_1`$ belonging to that of $`\xi _1`$. This effect has indeed been reported in Ionova and Carter (1993).

In the literature, plenty of methods exist dealing with the computation of transition states. One kind of them starts from the knowledge of the initial and final minimum Elber and Karplus (1987); Matro et al. (1994); Deaven et al. (1996); Angelani et al. (2000b). After a more or less educated guess for an initial trial RP, one iteratively improves the RP according to some prescription, e.g., the minimization of an action functional. Two sources of erroneous results have to be addressed in this connection. First, the two minima in question have to be true neighbours. This can only be verified by locating two points close to the basin border, e.g. by interval bisection of the initial trial path. The numerical cost is not small; for our ridge method, for instance, about one third of the calculation time is consumed by fixing $`y_0`$ and $`y_1`$ (depending on the minimization interval of the MD run). Second, the iterative path optimization may become stuck in a local extremum, due to an unfortunate choice of the initial path.

The other kind of TS search methods start from an initial minimum and climb up to a transition state guided by the shape of the PEL. Just walking against the force, however, would be a fatal strategy, as one can see by turning the PEL upside down: ending up in a TS is numerically impossible, since one quickly runs into one of the PEL singularities (two or more identical particle positions). Eigenvector-following algorithms Wales (2001) overcome this defocussing of steepest ascent paths by walking into the direction of negative local PEL curvature. The ’activation-relaxation technique’ by Mousseau and coworkers, in contrast, steps against the force in the direction leading away from the minimum, while descending the PEL perpendicular to that direction Barkema and Mousseau (1996). A drawback of the latter methods is that the choice for the next TS to mount is not well under control. From the minimum, a starting direction is chosen, either by purely random displacements or by some hard-sphere-like particle moves Doye and Wales (2001). Unfortunately, the number of escape directions from a minimum is generally very large (at least $`O(Nd)`$ as we found in the BMLJ65), whereas the majority of these is dynamically inaccessible at low $`T`$. Hence, eigenvector-following and activation-relaxation techniques yield many TSs which only negligibly contribute to relaxation rates. Striving for the simulation of low-temperature hopping dynamics based on these methods Hernandez-Rojas and Wales (2001); Mousseau (2000); Ball and Berry (1999), one may suffer a considerable reduction of efficiency.

We finally mention two complementary means of studying energy barriers. The ’lid’ algorithm, proposed by Schön and coworkers Wevers et al. (1999), is able to find upper bounds for the depths of single basins. By performing random walks below different potential energy thresholds and by regular minimizations, one is able to compute the elevation necessary for transitions to neighboring minima. From a more theoretical perspective, Schulz has specified a relation between transition rates and the overlap of vibrations in neighbouring basins Schulz (1998).

### IV.2 Comparison to $`\stackrel{~}{V}`$-saddles.

The advantage of the ridge method is that we definitely find the relevant barrier for a transition, i.e. a first order saddle on the basin border next to the point where the border was crossed. In contrast, the method using the auxiliary potential $`\stackrel{~}{V}=\frac{1}{2}|F(x)|^2`$ as applied in recent studies Angelani et al. (2000a); Angelani et al. (2002); Broderix et al. (2000); Grigera et al. (2002) has two major drawbacks. First, the $`\stackrel{~}{V}`$ minimization locates saddles (we call them $`\stackrel{~}{V}`$-saddles), even if they are not accessible kinetically. This is because the expression $`F^{}HF`$ is not positive ($`H=H(x)`$ denotes the Hessian of $`V(x)`$), i.e. $`\stackrel{~}{V}`$-minimization can climb up to a saddle. Second, one obtains higher-order saddles and, most frequently, non-stationary points (shoulders). These configurations are of no use to us because we specifically analyze paths over the lowest barriers on basin borders, i.e., transition states.

To shed more light on the interrelation of TSs and $`\stackrel{~}{V}`$-minima, we minimized $`\stackrel{~}{V}`$ by steepest descent, starting from configurations $`x(t)`$ only if $`\xi (t)\xi (t+1\text{ MD step})`$ (like $`y_0`$ in Fig. 5). In other words, we calculated $`\stackrel{~}{V}`$-saddles exactly at transition times. If this yielded the correct TSs, our more time-consuming ridge method would be clearly useless. The difference $`\mathrm{\Delta }ϵ=ϵ_{\stackrel{~}{V}}ϵ_{\text{TS}}`$ specifies the overestimation of the true barrier by the $`\stackrel{~}{V}`$-saddle. It may also happen that the index of the $`\stackrel{~}{V}`$-saddle is different from one. The distributions of $`\mathrm{\Delta }ϵ`$ and the index are shown in Fig. 6 ($`T=0.5`$). Obviously, the $`\stackrel{~}{V}`$-saddles considerably overestimate barriers and the correct TSs are only found very rarely. Moreover, most of the $`\stackrel{~}{V}`$-saddles have an index different from one, i.e. are no TSs at all. In turn, the energy of the TS is never undersold by a $`\stackrel{~}{V}`$-saddle. In conclusion, $`\stackrel{~}{V}`$-saddles turn out to have the undesired quality of being decorrelated from the relevant TSs, i.e., from the barriers that control relaxation (see section VI).

### IV.3 Population of Basin Borders.

After Angelani and coworkers Angelani et al. (2000a); Angelani et al. (2002), the mean index of $`\stackrel{~}{V}`$-saddles vanishes at $`T_c`$. Therefore, as they have argued, dynamics above $`T_c`$ is dominated by saddles, in that there are always some unstable directions available which allow the system to relax, without traversing an additional energy barrier. Passing $`T_c`$, the mechanism suffers a drastic change, and abruptly, one is faced with an index of ca. zero, i.e., saddles have to be reached via thermal activation. Since the preceeding subsection may cast some doubts on the significance of $`\stackrel{~}{V}`$-saddles, we now want to discuss an alternative analysis of the way the population of minima versus unstable configurations evolves upon decreasing temperature. More specifically, we determine the population of basin borders,

$$p_{\text{BB}}(T)=\frac{1}{Z(T)}ddxe^{\beta V(x)}\delta (x),$$
(8)
where integration is over the non-crystalline part of configuration space, also in the partition function $`Z(T)`$, and $``$ runs over all basin borders of the PEL. This expression is impractical in numerical simulation; one may rather ask if, for some instantaneous configuration $`x`$, there is a basin border nearby. In this case, small random displacements (length $`\delta 𝐈𝐑`$, direction $`\omega 𝐈𝐑^{Nd}`$, $`|\omega |=1`$) possibly lead into another basin, i.e. $`\xi (x)\xi (x+\omega \delta )`$. This kind of PEL analysis has been recently carried out by Fabricius and Stariolo Fabricius and Stariolo (2002). One calculates

$$p_{\text{BB}}(T;\delta )=P\left(\xi (x)\xi (x+\omega \delta )\right)_{T,\omega },$$
(9)
which is the probability that random disturbances $`\omega \delta `$ will cause crossings of basin borders at temperature $`T`$. The brackets denote the canonical plus the average over the random directions $`\omega `$.

One obtains the behavior

$$p_{\text{BB}}(T;\delta )\text{const}\times p_{\text{BB}}(T)\delta ,\delta 0$$
(10)
(the constant is set to unity for convenience). The validity of Eq. 10 is demonstrated in the left inset of Fig. 7, where $`p_{\text{BB}}(T;\delta )/\delta `$ has been calculated as a function of $`\delta `$. We find that $`p_{\text{BB}}(T;\delta )/\delta `$ is constant within statistical error below $`\delta =1.2`$. As an orientation, the typical distance between neighboring minima is larger than $`2.0`$, whereas intra-MB neighbors on average are less than $`1.0`$ apart.

The main part of Fig. 7 shows results for $`p_{\text{BB}}(T)`$ in an Arrhenius plot, with $`\delta =0.7`$. Over the whole temperature range considered, $`p_{\text{BB}}(T)`$ is Arrhenius-like. The apparent activation energy is ca. 1.8, which is small in comparison with the typical values observed for MB lifetimes. However, the temperature dependence becomes stronger if we impose the constraint of a minimum distance between neighbouring minima (data not shown). In this way, we eliminate the fast intra-MB transitions, which have small barriers.

In any event, $`p_{\text{BB}}(T)`$ features no noticeable change in behavior when approaching and crossing $`T_c`$. In a different graphical representation (see right inset) one might wrongly conclude that $`p_{\text{BB}}(T)`$ disappears at some finite temperature. Stated differently, the data suggest that the increasing timescale separation upon cooling happens rather smoothly, with no distinctly new physics emerging near $`T_c`$. This is in qualitative agreement with the work of Schrøder et al. Schrøder et al. (2000), who use the incoherent scattering functions from hopping dynamics $`\xi (t)`$ to deal with the separation of intra- and interbasin dynamics. There, the initial short-time decay of scattering functions (quantified by the so-called non-ergodicity parameter) is nothing else than a measure for the population of basin borders.

## V Energy Barriers from PEL Topology.

### V.1 Return Probabilities and Metabasin Definition.

With the tools of interval bisection and TS search, we are now in the position to analyze the escapes from MBs in full detail. When a MB is left, we first resolve all minima visited during the escape. Second, all corresponding TSs and, if desired, reaction paths are calculated. An example is shown in Fig. 8. The successive RPs were spliced together to a long, multi-minima RP $`\zeta (s)`$. One might take the energy profile, $`V(\zeta (s))`$, depicted in the figure, for one of the common cartoons of a PEL. However, it rests upon real data. Berry and coworkers have produced similar charts for the relaxation of small atomic clusters towards their global minima Ball et al. (1996); Ball and Berry (1999). For $`s<5`$ one can see the typical back-and-forth hopping among the ground minima of the MB. Obviously, the corresponding barriers are not large compared to $`k_\mathrm{B}T=0.5`$. The escape starts at $`s=5`$. The first minimum reached is very unstable as expected from the small backward barrier. Indeed, if we repeatedly start in this minimum and perform a number of short simulation runs (here: $`99`$) with different random numbers, the system will return to the bottom of the MB with probability $`p_{\text{back}}=98\%`$ and leave the range of attraction only rarely. Thus, the escape is far from being complete at this stage. Going to the next minimum, the return probability decreases, but does not drop to zero. We say that the system is free if $`p_{\text{back}}`$ is smaller than $`50\%`$. As the outcome of this investigation, we obtain the energy barrier surmounted before the minimum with $`p_{\text{back}}<50\%`$ was reached, see below. The exits from other long-lived MBs mostly look the same as in the example, while the escape in one jump is not common. In other words, MBs usually have the form of a funnel with some ledges on the walls Stillinger (1995); Middleton and Wales (2001). Minima with $`p_{\text{back}}>50\%`$ are said to belong to the MB. This criterion is reminiscent of the definition of dynamic bottlenecks introduced by Chandler and coworkers Bolhuis et al. (2002).

An interesting property of a MB is its diameter $`d`$. It is defined as the maximum distance between its minima. For the MBs found in the simulation at $`T=0.5`$, the distribution of diameters is depicted in Fig. 9. The delta-peak from single-minimum MBs has been omitted. No MB with $`d>d_{\text{max}}=4`$ has been found. As a consequence, if a minimum has a distance larger than $`d_{\text{max}}`$ to some MB minimum, we can safely assume $`p_{\text{back}}50\%`$. This criterion has already been used in section III.

Based on these insights, we can now provide a more complete description of MBs (Fig. 10). First, the ground state of a MB has to be identified (kernel minimum), since the definition of $`p_{\text{back}}`$ rests upon it. At low enough temperatures, the kernel minimum will certainly be visited during the MB lifetime, due to the very low barriers among the minima on the bottom of the MB. Second, for minima beyond the distance $`d_{\text{max}}`$ from the kernel, we set $`p_{\text{back}}`$ to zero. Third, the probability $`p_{\text{back}}`$ for returning to the kernel before reaching a distance greater than $`d_{\text{max}}`$ can be assigned to every remaining minimum and -in principle- be computed by simulation. To this end, one repeatedly starts in the minimum and checks if a recurrence to the kernel occurs. Fourth, the minima with $`p_{\text{back}}>50\%`$ are defined as the MB.

Please bear in mind that $`p_{\text{back}}`$ will in general depend on temperature, since it is defined by dynamics. Correlations among minima are expected to increase towards lower temperatures, implying that MBs are no static concept but rather grow with decreasing $`T`$. In Fig. 8, e.g., the minimum at $`s6.5`$ has the ’critical’ value of $`p_{\text{back}}47\%`$ at $`T=0.5`$. Although we do not know the details of PEL connectivity around this minimum, the small backward barrier suggests that the minimum would exceed $`p_{\text{back}}=50\%`$ for still lower temperatures, thus joining the MB. However, we may also conceive some situations where a critical $`p_{\text{back}}50\%`$ is quite unsusceptible to temperature changes. This is the case if backward and forward barriers are of about the same size. We will come back to that issue later.

We further note that the explicit computation of $`p_{\text{back}}`$ can be extremely expensive. This is mainly the case when $`p_{\text{back}}`$ is small, and complete escapes beyond $`d_{\text{max}}`$ have to be awaited. However, the exact value of $`p_{\text{back}}`$ is of no great interest. In fact, it suffices to know whether $`p_{\text{back}}<50\%`$ or $`p_{\text{back}}>50\%`$. This decision can often be reached to a high confidence with few trials.

The MB lifetime algorithm in section III is based on the detection of back-and-forth jumps between minima. One mostly observes the dominant minima on the bottom of the MBs, whereas the more elevated members are only weakly populated, see Fig. 10. If MB lifetimes are to be read from a simulation run, it suffices to notice when the set of dominant MB minima has been left, since the visits to the elevated minima at the end of the MB lifetime happen quite rapidly. Thus, the algorithm of section III reduces the MB to the most populated minima, which is sufficient for the purpose of lifetime calculation from a given simulation run. In contrast, for the prediction of MB relaxation behavior as pursued in this section, the minima close to the rim of MBs are of special interest. Their elevations from the bottom of the MB give the depth of the MB.

### V.2 Barriers for Metabasin Relaxations.

In the spirit of the above remarks, we will now carry out a systematic investigation of the energy barriers overcome when escaping MBs. The goal to is to recover the apparent activation energies computed in section III from PEL topology.

The mean lifetime $`\tau _i`$ of MB $`i`$ can be expressed in terms of escape rates $`\gamma _{i,\alpha }`$ of different relaxation channels $`\alpha `$,

$$\tau _i^1=\underset{\alpha }{}\gamma _{i,\alpha }.$$
(11)
In general, each $`\gamma _{i,\alpha }`$ reflects a multi-minima escape path

$$\xi _0\stackrel{\zeta _{01}}{}\xi _1\stackrel{\zeta _{12}}{}\xi _2\mathrm{}\xi _{M1}\stackrel{\zeta _{M1,M}}{}\xi _M$$
(12)
as the one shown in Fig. 8. Here, $`\xi _0`$ is the kernel minimum and $`\zeta _{ab}`$ is the TS for $`\xi _a\xi _b`$. Suppose that the number $`M`$ of jumps in the sequence Eq. 12 is large enough to completely quit the MB’s range of attraction, i.e., $`p_{\text{back}}(M)0`$. For the escape shown in Fig. 8, e.g., $`M7`$ would be fine.

We further take for granted that the rates for single barrier crossings follow quantitatively -via transition state theory- from the height of barriers, $`E_{ab}=V(\zeta _{ab})V(\xi _a)`$ (the energy difference between the minimum $`\xi _a`$ and the TS between $`a`$ and $`b`$). Hence, rates $`g_{ab}`$ for single transitions $`\xi _a\xi _b`$ are characterized by

$$g_{ab}e^{\beta E_{ab}}.$$
(13)
A justification for this assumption, even for temperatures above $`T_c`$, will be given in section VI.

Generally, the probability of upward jumps is small at low $`T`$. Hence, climbing out of a MB in a back-and-forth fashion (e.g., $`\xi _a=\xi _{a+2}`$ and $`\xi _{a+1}=\xi _{a+3}`$) is not probable. (This is reminiscent of the fact that the activated crossing of single potential barriers happens on a short time scale, i.e. in a rather straight way.) In contrast, excursions from the main path may happen. As shown in Fig. 8, the minimum at $`s=6.5`$ is revisited at $`s=9`$ after taking a look at another minimum ($`s8`$). The latter does not appear again later on. Clearly, running into such ’dead ends’ should not contribute to the escape rate via the successful main path. We therefore eliminate such excursions from the sequence of minima, Eq. 12. From these remarks we take the liberty of assuming that no minimum appears more than once along the escape path,

$$\xi _a\xi _b,ab.$$
(14)
We are now interested in the contribution of the path Eq. 12 to the total escape rate Eq. 11. Particularly, we have to consider the question of how many single transitions are relevant for the escape process. The probability to jump from minimum $`\xi _a`$ to $`\xi _{a+1}`$ is $`g_{a,a+1}/g_a`$, where $`g_a`$ denotes the inverse lifetime of minimum $`\xi _a`$. The rate of escape via a longer pathway now is given by the rate of the first jump times the probability that the minima $`\xi _a`$ ($`a=1,\mathrm{}M`$) are visited in correct order thereafter,

$$\gamma _{i,\alpha }=g_{01}\frac{g_{12}}{g_1}\frac{g_{23}}{g_2}\mathrm{}\frac{g_{M1,M}}{g_{M1}}.$$
(15)
With the help of Eq. 13 one calculates

$$\frac{\mathrm{d}}{\mathrm{d}\beta }\mathrm{ln}\gamma _{i,\alpha }=E_{01}+\underset{a=1}{\overset{M1}{}}p_{\text{ret}}(a)(E_{a,a+1}E_{a,a1}),$$
(16)
where $`p_{\text{ret}}(a)=g_{a,a1}/g_a`$ is the probability to jump back to minimum $`a1`$ from minimum $`a`$. In the derivation of Eq. 16, we have neglected a term proportional to $`E_{a,a+1}`$ minus the average barrier when jumping from $`a`$ to a neighbouring minimum other than $`a1`$. This term strictly vanishes when performing the final summation in Eq. 11. Moreover, we made use of Eq. 14.

One possibility for calculating activation energies from Eq. 16 would be to consider the complete paths Eq. 12, where $`p_{\text{back}}(M)0`$, and determine all terms in the sum of Eq. 16. However, an accurate computation of all the desired $`p_{\text{ret}}(a)`$’s would even be more costly than the determination of the point where $`p_{\text{back}}`$ changes from above to below $`50\%`$. We therefore use the following approximation of Eq. 16, which is in conformance with our previous definition of MBs: Let $`m(T)`$ be the first minimum along the path Eq. 12, where $`p_{\text{back}}<50\%`$. Then, for all $`a<m(T)`$, we set $`p_{\text{ret}}(a)`$ to unity, while for $`am(T)`$ (i.e. outside the MB), we let $`p_{\text{ret}}(a)=0`$. Thus,

$$\begin{array}{cc}\hfill \frac{\mathrm{d}}{\mathrm{d}\beta }\mathrm{ln}\gamma _{i,\alpha }E_{i,\alpha }& E_{01}+\underset{a=1}{\overset{m1}{}}(E_{a,a+1}E_{a,a1})\hfill \\ & =ϵ_{m1}ϵ_0+E_{m,m1},\hfill \end{array}$$
(17)
where $`m=m(T)`$. In this way, the terms $`a<m(T)`$ in Eq. 16 are given higher weights, whereas those of $`am(T)`$ are neglected. We will dwell on the quality of this approximation later on.

Note that, due to the temperature dependence of $`p_{\text{back}}`$, energy barriers $`E_{i,\alpha }`$ generally increase upon cooling: At high temperatures, in contrast, correlations among minima are small, such that MBs (even the low-lying) consist of only one minimum. This effect is included in Eq. 17 by the temperature dependence of $`m(T)`$.

### V.3 Single Metabasins.

We now relate the lifetimes of single, selected MBs (cf. section III.2) to PEL barriers. By repeated starts from these MBs, the local PEL topography is sampled thoroughly, yielding sets of typical escape pathways. Whenever a MB is left, we locate the transitions by interval bisection and obtain the corresponding TSs with the help of the ridge method. Then, $`p_{\text{back}}`$ is calculated for the minima visited, until for the first time, $`p_{\text{back}}<50\%`$. Finally, the barrier $`E_{i,\alpha (k)}`$ is computed according to Eq. 17, where $`\alpha (k)`$ denotes the escape path chosen at the $`k`$th escape.

The histograms of barriers are shown in Fig. 11, for the four MBs of Fig. 1, at $`T=0.5=1.1T_c`$. Due to the slow dynamics at this temperature, the computation of $`p_{\text{back}}`$ was rather expensive. Nevertheless, the statistics should be sufficient for a reasonable estimate of the apparent activation energy. To this end, we express $`E_{\text{app}}(i)`$ of MB $`i`$ in terms of the contributions $`E_{i,\alpha }`$,

$$\frac{\mathrm{d}}{\mathrm{d}\beta }\mathrm{ln}\tau _i\tau _i\underset{\alpha }{}E_{i,\alpha }\gamma _{i,\alpha }=\underset{\alpha }{}p_{i,\alpha }E_{i,\alpha }E_{\text{app}}^{\text{est}}(i),$$
(18)
where Eqs. 11 and 17 were used. Thus, the barriers $`E_{i,\alpha }`$ are weighted by the probabilities $`p_{i,\alpha }=\gamma _{i,\alpha }/_\alpha \gamma _{i,\alpha }`$ that the escape happens via pathway $`\alpha `$. Note that the $`E_{i,\alpha (k)}`$ correspond to the pathways that were chosen by the system, i.e. they are already weighted correctly by $`p_{i,\alpha (k)}`$, compare Eq. 18. Therefore, $`E_{\text{app}}^{\text{est}}(i)`$ is just the average of the $`E_{i,\alpha (k)}`$. The values of $`E_{\text{app}}(i)`$ and $`E_{\text{app}}^{\text{est}}(i)`$, given in Fig. 11, are in good agreement. Also shown in Fig. 11 is the distribution of first barriers, $`E_{01}`$, belonging to the step $`\xi _0\xi _1`$. Evidently, the neglect of the multi-minima nature of escapes leads to a considerable underestimation of apparent activation energies.

We now continue the discussion of the temperature dependence of barriers $`E_{i,\alpha }(T)`$. At the example of MB 1 from Fig. 11, we have carried out the above program for two other temperatures, $`T=0.6`$ and $`0.8`$. The obtained distributions of barriers, $`P(E_{i,\alpha })`$, are shown in Fig. 12. We find that the estimates for the apparent activation energy ($`E_{\text{app}}^{\text{est}}(1)=6.9\pm 0.5`$, $`T=0.6`$, and $`E_{\text{app}}^{\text{est}}=6.8\pm 0.5`$, $`T=0.8`$) remain in good agreement with $`E_{\text{app}}(1)=6.7\pm 0.3`$ from section III. The distributions of barriers, however, grow narrower with decreasing temperature. Single, high barriers, contributing to the right wing of the distribution, become inaccessible at low $`T`$, i.e., the relative weights $`p_{i,\alpha }`$ of the corresponding escapes become small. This suppression of high barriers at low $`T`$ is a trivial effect.

More interesting is the vanishing of small barriers upon cooling, i.e., of the barriers $`E<5`$ in the figure. Naively, one would expect these to dominate the escape rate at low $`T`$. However, due to the stronger backward correlations (increased $`p_{\text{back}}`$), jumps over these barriers eventually do not suffice anymore to escape. As described above, the respective escape paths, $`\xi _0\mathrm{}\xi _{m(T)}`$, grow longer, and the barriers change to a different, mostly larger value.

### V.4 Average over Metabasins.

During our analysis of the escape times in section III the apparent activation energies $`E_{\text{app}}(ϵ_{\text{MB}})`$ emerged as useful quantities. Although the above results already indicate that barrier hopping is the relevant motional mechanism, a clear-cut verification requires the comparison with the average barrier the system has to cross when leaving a MB with energy $`ϵ_{\text{MB}}`$.

For this purpose we now carry out a similar program as before, with many MBs visited during an ordinary MD run. We concentrate on MBs with lifetimes of more than $`10^5`$ MD steps (179 MBs) at $`T=0.5`$. When such a MB is left, we locate the transitions by interval bisection and obtain the corresponding TSs by the ridge method. Then, we calculate $`p_{\text{back}}`$ and identify the barrier $`E_kE_{i(k),\alpha (k)}`$ according to Eq. 17. The histogram of barriers is shown as the bold line in Fig. 13. For comparison, we also show the barriers minus the contribution of the TSs $`E_{(m1)m}`$. Ignoring multi-minima correlations, we further show the histogram of first barriers $`E_{01}`$ of escapes. Evidently, the neglect of TSs or of backward correlations leads to much smaller barriers.

From the above barriers, we will now calculate estimates of apparent activation energies. When the average over lifetimes of different MBs is considered, each MB $`i`$ acquires a weight $`\phi _i`$ corresponding to its probability of occurrence,

$$\tau =\underset{i}{}\phi _i\tau _i.$$

At fixed $`ϵ_{\text{MB}}`$, the analog to Eq. 18 can then be derived

$$\frac{\mathrm{d}}{\mathrm{d}\beta }\mathrm{ln}\tau (ϵ_{\text{MB}};T)\underset{i}{}\frac{\tau _i\phi _i}{\tau (ϵ_{\text{MB}};T)}\underset{\alpha }{}p_{i,\alpha }E_{i,\alpha },$$
(19)
where summation goes over MBs of energy $`ϵ_{\text{MB}}`$. As in Eq. 18, the barriers in Eq. 19 are weighted according to their probability of occurrence, but, additionally, with the respective MB lifetimes.

In Eq. 19, we have neglected terms stemming from the variation of $`\phi _i`$’s with temperature. This is justified, since the $`\phi _i`$’s belong to the same $`ϵ_{\text{MB}}`$. Their relative weights will only vary if these MBs differ considerably in barrier heights. As already stated above, however, MBs of the same energy seem to be fairly uniform regarding this property. For the finite sample of MBs visited during an MD run, Eq. 19 then takes the form

$$E_{\text{app}}^{\text{est}}(ϵ_{\text{MB}})=\frac{\tau _kE_k}{\tau _k},$$
(20)
where summation goes over MBs of energy $`ϵ_{\text{MB}}`$. Again, the correct weighting is implicit here. This expression can be shown to converge to the right-hand side of Eq. 19 in the limit of infinitely long sampling.

In Fig. 14 we show the values of $`E_{\text{app}}^{\text{est}}(ϵ_{\text{MB}})`$, determined in this way. They perfectly agree with the apparent activation energies, derived from the analysis of relaxation times at different temperatures. Thus we have a clear-cut proof that the apparent activation energies $`E_{\text{app}}(ϵ_{\text{MB}})`$ are indeed related to barriers on the PEL and thus reflect activated behavior significantly above $`T_c`$. This again demonstrates that we not only deal with the right order of barrier sizes, but we also quantitatively link PEL topography to dynamics.

For comparison, we included the apparent activation energy which results, if only the first transitions of escapes, $`\xi _0\xi _1`$, are considered ($`E_{01}=V(\zeta _{01})ϵ_0`$). One ends up with much too small apparent activation energies. Again, multi-minima correlations turn out to be crucial for the characterization of MB depths.

In principle, the results of Fig. 14 may slightly change if all MBs rather than those with lifetimes larger $`10^5`$ MD steps were considered. However, our analysis has clearly revealed (see, e.g., Fig. 1) that the depth of the trap only mildly varies when comparing MBs with similar $`ϵ_{\text{MB}}`$. Thus inclusion of MBs with smaller values of $`\tau `$ would not significantly change the values of the apparent activation energies $`E_{\text{app}}^{\text{est}}(ϵ_{\text{MB}})`$.

Finally, we show that these results, in conjunction with $`p(ϵ_{\text{MB}};T)`$, largely explain the behavior of the diffusion constant $`D(T)`$. This is a conceptually important step, since we link $`D(T)`$ to purely structural and thermodynamical quantitites, see Eq. 5. The key is the mean lifetime $`\tau (ϵ_{\text{MB}};T)`$ of MBs at energy $`ϵ_{\text{MB}}`$, which is parametrized by $`\tau _0(ϵ_{\text{MB}})`$ and $`E_{\text{app}}(ϵ_{\text{MB}})`$ (Eq. 7). The former, $`\tau _0(ϵ_{\text{MB}})`$, however, has not been deduced from PEL properties. Its variation with MB energy is not strong (Fig. 2(c)), so we can hope that setting it to a constant will be a good approximation. Thus, Eq. 4 becomes

$$D(T)\frac{a^2}{6N\tau _0}dϵ_{\text{MB}}p(ϵ_{\text{MB}};T)e^{\beta E_{\text{app}}^{\text{est}}(ϵ_{\text{MB}})}D_{\text{est}}(T).$$
(21)
The estimated diffusion constant derived from this expression is shown in Fig. 15. The agreement of $`D(T)`$ with our estimate is satisfactory below $`T=1`$, albeit we find a slightly too strong temperature dependence for the lowest $`T`$. The deviation at $`T=2`$ is due to the depart of $`\tau (ϵ_{\text{MB}};T)`$ from Arrhenius behavior, see Fig. 2(a).

## VI Barrier Crossing.

When making use of Eq. 13, we presumed that the barriers $`V(\zeta _{ab})V(\xi _a)`$ in fact are the determinants of the temperature dependence of rates. The excellent agreement between $`E_{\text{app}}(ϵ_{\text{MB}})`$, determined from dynamics, and the $`E_{\text{app}}^{\text{est}}(ϵ_{\text{MB}})`$, from the analysis of the PEL, strongly indicates that this presumption is indeed true. We will show here in a very detailed way that at $`T=0.5=1.1T_c`$, escapes from stable MBs are perfectly activated. More precisely, two conditions are fulfilled, (i) the potential barriers are much larger than $`k_\mathrm{B}T`$, (ii) rates follow from the 1D energy profile of the RP plus corrections from perpendicular curvatures.

We will check these conditions explicitly here by an analysis of escape dynamics out of MBs. We made the observation that during every escape from a stable MB, at least one single barrier larger than $`6k_\mathrm{B}T`$ must be surmounted. Moreover, this larger jump is mostly undertaken from one of the lowest minima of the MB, compare Fig. 8. From the repeated escape runs of section III.2, we selected the most frequent ten transitions of that kind. From the respective TSs, $`\zeta _l`$, we computed the RPs, denoted $`\zeta _l(s)`$, $`l=1\mathrm{}10`$. We then investigated the motion within the MBs over a long period of the simulation where no escape had happened ($`10^7`$ MD steps each MB). The goal was to observe how the system tries to climb the different RPs. To this end, we projected the instantaneous configuration $`x(t)`$ onto each of the RPs, according to

$$s_l(t)\{s^{}:x(t)\zeta _l(s^{})=\underset{s}{\mathrm{min}}x(t)\zeta _l(s)\},$$

which means the point on the RP next to $`x(t)`$. Due to the long residences in the MBs, motion therein is largely equilibrated. Hence, if the potential energy profiles $`V\left(\zeta _l(s_l)\right)`$ along the reaction paths are of importance for the transition rates, we expect that the populations $`p_l(s_l)`$ of the RPs follow from Boltzmann’s law

$$p_l(s_l)\mathrm{exp}\left\{\beta V_l(s_l)\right\}Y_l^{}(s_l)\mathrm{exp}\left\{\beta F_l(s_l)\right\}.$$

The vibrations perpendicular to path $`l`$ are accounted for by the harmonic partition function

$$Y_l^{}(s_l)=dx\mathrm{exp}\left\{\frac{\beta }{2}\lambda _\nu x_\nu ^2\right\}\delta (x\widehat{t}(s_l)),$$

where $`\lambda _\nu `$ are the eigenvalues of the Hessian matrix, $`H(s_l)`$, $`x_\nu `$ the components of $`x`$ along the eigenvectors, and $`\widehat{t}(s_l)`$ is the tangent to the reaction path.

The upper inset of Fig. 16 shows an example of $`p_l(s_l)`$ vs. $`F_l(s_l)/k_\mathrm{B}T`$. The population of the reaction path follows nicely the prediction from its energy profile. For RPs with complicated shapes, this correspondence can be disturbed. The worst agreement of the considered RPs is shown in the second inset. Still, a clear correlation of RP population with energy is present. We compiled the results for all 10 RPs in Fig. 16 as a parametric plot of $`\mathrm{ln}p_l`$ vs. $`F_l/k_\mathrm{B}T`$. Curves of slope one result from a perfect equivalence of $`p_l`$ to $`F_l/k_\mathrm{B}T`$. Here, we find an average slope of $`0.92`$. Since transition rates are proportional to the population of TSs, the implications of these results are obvious: MB jump rates follow from energy barriers. We finally note that the vibrational terms $`\mathrm{ln}Y_l^{}(s_l)`$ are minor as compared to $`\beta V(s_l)`$.

In view of these results, it is a little surprising that the TS location with the help of the auxiliary potential $`\stackrel{~}{V}`$ was that unsuccessful (cf. section IV). Since the RP population suits well the harmonic description of the RP, one expects that motion near the TS is quite harmonic, too. Minimizing $`\stackrel{~}{V}`$ in a harmonic potential directly yields the stationary state. Consequently, one should easily find the TS when starting from a configuration at $`s0`$. After section IV, this is not the case, so at least minor anharmonicities must be present.

## VII Discussion.

The metabasin concept is at the heart of the present study. The important insight is that, upon cooling, not only the time scale of inter-basin transitions becomes well separated from intra-basin vibrations, but also that a similar separation occurs between MB hopping and intra-MB transitions. Recently, Biroli and Kurchan have analysed the general problem of defining metastable states in glassy systems Biroli and Kurchan (2001). They conclude that one has no absolute notion of a state without making reference to a time scale and hence to dynamics. Also the present definition of MBs relies on the dynamics of the system. It is, however, independent of time scale and exclusively depends on the $`p_{\text{back}}`$ values which directly reflect the topological properties of the PEL.

Our MB definition (see Fig. 10) is devised to eliminate the information on trivial back-and-forth jumps within MBs. This strongly correlated type of motion is reminiscent of the particles’ rattling in the cages formed by their neighbors. Similarly, escaping from MBs seems to be equivalent to the breaking of cages and thus to structural relaxation. Guided by this idea, we have examined MB relaxation in great detail:

First, for repeated relaxation from the same MB, we calculated the mean relaxation time $`\tau _i`$ and found Arrhenius behavior in all cases. The simplest view is that the apparent activation energies $`E_{\text{app}}(i)`$ from the Arrhenius-like $`\tau _i`$ (Fig. 1) should correspond to the depths of these MBs, i.e. to the typical heights of barriers that surround the MBs. Indeed, this has been quantitatively confirmed for the four randomly selected, low-lying MBs (see Fig. 11). A direct conclusion from the constancy of $`E_{\text{app}}(i)`$ is that the system does not find smaller and ever smaller barriers upon decreasing $`T`$.

Although not of statistical relevance for the whole PEL, the results for the four single MBs give us a detailed picture of the local PEL topography. An important outcome is the variation of barrier heights with temperature, see Fig. 12. We have already discussed that low barriers increase upon cooling, due to enhanced multi-minima correlations (growing MBs), while unnecessarily high barriers are suppressed. Both effects seem to cancel, so that the mean barrier, $`E_{\text{app}}^{\text{est}}`$, remains constant, leading to Arrhenius behavior below $`T=1`$. This cancellation may be fortunate, at least we can offer no explanation for it, here. As depicted in Fig. 12, the distribution of barriers becomes more and more narrow when going from $`T=0.8`$ to $`T=0.5`$, but the mean value, i.e. $`E_{\text{app}}^{\text{est}}(1)`$, remains constant. The constant apparent activation energy of MB 1 down to $`T=0.45`$ implies that the mean value of the distribution of barriers has not increased. We thus speculate that the growth of barriers due to increasing multi-minima correlations has essentially come to an end at $`T0.5`$. Although the temperature dependence of the barrier distribution has only been analyzed for a single MB, the constancy of apparent activation energies of the other three MBs and the temperature independence of $`E_{\text{app}}(ϵ_{\text{MB}})`$ support this idea. Stated differently, the development of superstructures of minima seems to cease at some temperature above $`T_g`$. Expressed by $`p_{\text{back}}`$, this means that no minimum with $`p_{\text{back}}<50\%`$ will surpass $`p_{\text{back}}=50\%`$ upon further cooling, thus being unable to join the MB in question. Hence, an escape sequence found at one temperature $`T0.5`$ has the same length at another one, i.e., from some temperature on, the minimum $`\xi _{m(T)}`$ remains at $`p_{\text{back}}<50\%`$ for $`T0`$; we then say it terminates the sequence. It is an interesting question under what circumstances such termination happens. A trivial example would be a ’transit’ minimum with one backward and one forward exit, where taking the forward leads to a minimum with $`p_{\text{back}}0`$. If the backward barrier was higher than the forward one, $`p_{\text{back}}`$ would go to zero for $`T0`$. On the other hand, the minima inside MBs generally feature growing $`p_{\text{back}}`$’s upon cooling, because the energetic gain of returning becomes more and more attractive. Ideally, thus, for $`T0`$, we would have $`p_{\text{back}}1`$ within MBs, and $`p_{\text{back}}0`$ outside. This provides a plausible, physical basis for computing barrier heights according to Eq. 17, at least in the limit $`T0`$. Clearly, a more detailed investigation of the temperature dependence of $`p_{\text{back}}`$ is necessary to back these conclusions.

Second, we analyzed the average relaxation times $`\tau (ϵ_{\text{MB}};T)`$ from MBs at fixed energy $`ϵ_{\text{MB}}`$. Again, they displayed Arrhenius behavior, with apparent activation energy $`E_{\text{app}}(ϵ_{\text{MB}})`$ (see Fig. 2), which compared well with the prediction from PEL barriers (Fig. 14). In this connection, a recent paper by Grigera et al. Grigera et al. (2002) is of interest. The authors use the $`\stackrel{~}{V}`$-potential to compute saddles in a binary soft-sphere mixture ($`N=70`$). From the TSs among these saddles (index one, no shoulder), they perform steepest descents to obtain the connected minima. They define barriers as the energy difference $`\mathrm{\Delta }U`$ from the TSs to the lower one of the connected minima $`ϵ=\mathrm{min}(ϵ_0,ϵ_1)`$. Plotting the average $`\overline{\mathrm{\Delta }U}(ϵ)`$, they find a similar curve to our $`E_{\text{app}}(ϵ_{\text{MB}})`$, Fig. 14, i.e., a strong increase of barriers towards lower energies. In contrast, when carrying out the same analysis for our BMLJ65, we found a nearly constant $`\overline{\mathrm{\Delta }U}(ϵ)`$, a curve close to the first barriers of escapes $`E_{01}`$ shown in Fig. 14. We would have expected this result, since the multi-step nature of escapes in the BMLJ65 has clearly been demonstrated. On the other hand, the contrasting result of Grigera et al. indicates that the soft-sphere PEL is not organized in multi-minima superstructures. A clarification of this point would be very useful.

Note that $`E_{\text{app}}(ϵ_{\text{MB}})`$ is of special importance since it bridges the separation between dynamics (diffusion constant $`D(T)`$) and thermodynamics (population of $`ϵ_{\text{MB}}`$). Clearly, an understanding of $`E_{\text{app}}(ϵ_{\text{MB}})`$ from basic principles is highly desirable. It is plausible that the simple form $`E_{\text{app}}(ϵ_{\text{MB}})=ϵ_{\text{th}}ϵ_{\text{MB}}`$ can only be expected for a system acting as a completely correlated entity. In contrast, two independently relaxing subsystems should generally show a weaker dependence of $`E_{\text{app}}(ϵ_{\text{MB}})`$ on $`ϵ_{\text{MB}}`$. This can be seen by a very simple argument. Consider two independent, identical systems, with MB energies $`ϵ_{\text{MB}}^{\text{(1)}},ϵ_{\text{MB}}^{\text{(2)}}`$ and activation energy $`\stackrel{~}{E}_{\text{app}}(ϵ_{\text{MB}}^{\text{(1,2)}})`$. What can be said about $`E_{\text{app}}(ϵ_{\text{MB}})`$ of the union of these systems, at MB energy $`ϵ_{\text{MB}}=ϵ_{\text{MB}}^{\text{(1)}}+ϵ_{\text{MB}}^{\text{(2)}}`$? In the limit of low temperatures, the apparent activation energy is given by $`\mathrm{min}(\stackrel{~}{E}_{\text{app}}(ϵ_{\text{MB}}^{\text{(1)}}),\stackrel{~}{E}_{\text{app}}(ϵ_{\text{MB}}^{\text{(2)}}))`$. A proper average over different realizations $`ϵ_{\text{MB}}^{\text{(1)}}`$, $`ϵ_{\text{MB}}^{\text{(2)}}=ϵ_{\text{MB}}ϵ_{\text{MB}}^{\text{(1)}}`$ yields $`E_{\text{app}}(ϵ_{\text{MB}})`$ of the combined system. Instead of carrying out this average, we use the fact that $`\stackrel{~}{E}_{\text{app}}(ϵ_{\text{MB}})`$ is a monotonous function and estimate

$$0\mathrm{min}(\stackrel{~}{E}_{\text{app}}(ϵ_{\text{MB}}^{\text{(1)}}),\stackrel{~}{E}_{\text{app}}(ϵ_{\text{MB}}^{\text{(2)}}))\stackrel{~}{E}_{\text{app}}(ϵ_{\text{MB}}/2).$$

Thus, $`0E_{\text{app}}(ϵ_{\text{MB}})\stackrel{~}{E}_{\text{app}}(ϵ_{\text{MB}}/2)`$, which means that the combined system shows a weaker dependence on $`ϵ_{\text{MB}}`$ than a single copy. For a reasonable PEL topology, one would expect $`|\mathrm{d}E_{\text{app}}(ϵ_{\text{MB}})/\mathrm{d}ϵ_{\text{MB}}|1`$, because barriers should not mount up more than one descents in the PEL. Since the $`ϵ_{\text{MB}}`$-dependence of $`E_{\text{app}}`$ becomes weaker for larger systems, it in turn should increases towards smaller $`N`$. As a speculation, this might open a way of estimating the size of cooperative regions.

The results shown in Fig. 15, obtained via Eqs. 4 and 7, demonstrate the use of the present work. From PEL barriers ($`E_{\text{app}}^{\text{est}}(ϵ_{\text{MB}})`$) and thermodynamics ($`p(ϵ_{\text{MB}};T)`$) we are able to produce a reasonable estimate of dynamics. An overall proportionality factor $`1/\tau _0`$ remains as an adjustable parameter, since it could not be predicted from PEL structure. As discussed in section III, one may use $`p(ϵ;T)`$ instead of $`p(ϵ_{\text{MB}};T)`$, since they are nearly identical. This is very convenient, because upon contructing $`p(ϵ;T)`$, no information about dynamics is needed. The breakdown of the Arrhenius form of $`\tau (ϵ_{\text{MB}};T)`$ above $`2T_c`$ limits our description to the temperatures $`T2T_c`$. In any event, we would not have dared to make quantitative statements on the basis of the hopping picture above the landscape-influenced temperature regime. At $`T=1`$, for instance, we have $`ϵ_{\text{MB}}_T289`$ (Fig. 3), where we already find PEL barriers $`E_{\text{app}}(289)1=k_\mathrm{B}T`$ (Fig. 2).

From the fact that we could quantitatively relate MB lifetimes to PEL barriers below $`2T_c`$ and the results from section VI, we see that there exist activated barrier crossing events significantly above $`T_c`$. As shown before Doliwa and Heuer (2002), these escape processes from stable MBs determine the temperature dependence of the diffusion constant also above $`T_c`$. Thus, the general statement that hopping events are not relevant there (see, e.g. Franosch and Götze (1999)) is not correct for the BMLJ system. This implies that the ideal MCT can be applied to systems for which activated processes determine the time-scale of relaxation. Thus it seems that the theoretical description of the cage effect in terms of structural quantities, as done in MCT, works independent of whether the cage effect is purely entropic (like in hard-sphere systems) or is to a large degree based on barrier-crossing events.

Moreover, with the help of the unbiased quantity $`p_{\text{BB}}(T)`$, we were able to measure the population of basin borders. No indication for an abrupt change of relaxation mechanism could be observed in $`p_{\text{BB}}(T)`$; in contrast, the separation of intra- and inter-basin motion seems to happen rather smoothly (see Fig. 7). Thus, there is no qualitative change around $`T_c`$.

We finally discuss the relation of our work to the instantenous normal mode approach (INM) which considers the number of ’diffusive modes’, $`f_{\text{diff}}`$, to be at the physical basis of diffusion Nave et al. (2000, 2001, 2002). From the directions corresponding to negative eigenvalues of the Hessian $`H(x(t))`$ (unstable directions), one filters out the ’diffusive’ directions. Considering the energy profile on the straight lines along the unstable directions, La Nave et al. observed extremely small barriers, indicating completely ’entropic’ dynamics at the considered temperatures Nave et al. (2001). This conclusion, though reached for a model of supercooled water, is in contrast to our findings of the relevance of energetic barriers. A possible key to this apparent contradiction is that $`f_{\text{diff}}`$ is directly related to the fraction of time spent in ’mobile’ regions of configuration space. In contrast, we have concentrated on the durations of the stable, immobile structures. As the consequence of longer and longer residences in deep MBs, the mobile fraction becomes smaller and smaller. Thus, one observes a relation between $`D(T)`$ and $`f_{\text{diff}}`$, although it is the long trapping times which are the reason for the slowing down of dynamics.

We further note that the MB concept is not implemented in the INM approach. Supercooled water, e.g., exibits very pronounced MB correlations in the time series of minima, even for a ’large’ system of 216 particles Giovambattista et al. (2002). Generally, fragile glass formers are expected to have a ’rugged’ PEL, i.e. exibit extensive superstructures of minima Stillinger (1995). In view of this insight, the success of INM analyses for the latter type of systems is quite surprising.

## VIII Conclusion.

Our goal in this paper was to shed some light on the temperature dependence of the diffusion constant. In our previous work Doliwa and Heuer (2002), metabasins turned out as a useful concept that reduces correlations of subsequent PEL-hopping events. Taking seriously these correlations, the present investigation went a step further into this direction, by relating the temperature dependence of relaxation to the depths of these multi-minima superstructures. We have shown in this paper that a quantitative link between PEL structure and dynamics is possible above $`T_c`$. However, our approach is still phenomenological at this stage: we are far from predicting $`E_{\text{app}}(ϵ_{\text{MB}})`$ from more general PEL properties or even the interaction potentials. To achieve this is a major challenge, implying a deep understanding of energy landscape topology.

Further insights might be obtained by unveiling the real-space aspects of MB relaxation. Here, the correspondance of MBs to the cage effect should serve as the guiding principle. An interesting question along this line would be if some of the real-space phenomena found in supercooled liquids (e.g. the string-like motion discovered by Donati et al. Donati et al. (1998)) can be traced back to energy landscape features.

We are pleased to thank M. Fuchs, D.R. Reichman, and H.W. Spiess for helpful discussions. This work has been supported by the DFG, Sonderforschungsbereich 262.


text:
# What can exoplanets tell us about our Solar System?

## 1. Identifiable Trends in Exoplanet Data

Despite the fact that massive planets are easier to detect, the mass distribution of detected planets is strongly peaked toward the lowest detectable masses. And despite the fact that short period planets are easier to detect, the period distribution is strongly peaked toward the longest detectable periods. In Lineweaver & Grether (2002, hereafter LG) we quantified these trends as accurately as possible. Here we update this analysis by including the 27 exoplanets detected between January and August 2002. As in LG, we identify a less-biased subsample of exoplanets (thick rectangle of Fig. 1). Within this subsample, we correct for completeness and then quantify trends in mass and period (Fig. 2) that are less biased than trends based on the full sample of exoplanets. Straightforward extrapolations of these trends, into the area of parameter space occupied by Jupiter, indicates that Jupiter lies in a region densely occupied by exoplanets.

Naef et al. (2001) point out that none of the planetary companions detected so far resembles the giants of the Solar System. However, this observational fact is consistent with the idea that our Solar System is a typical planetary system. Fig. 1 shows that selection effects can easily explain the lack of detections of Jupiter-like planets. Exoplanets detected to date can not resemble the planets of our Solar System because the Doppler technique used to detect exoplanets has not been sensitive enough to detect Jupiter-like planets. We may be sampling the tail of a distribution – the only part that we are capable of sampling. If the Sun were a target star in one of the Doppler surveys, no planet would have been detected around it. This situation is about to change.

Our analysis suggests that Jupiter is more typical than indicated by previous analyses, including our own (LG). For example, in Fig. 2, our $`\alpha =1.6`$ slope is slightly steeper than the $`\alpha =1.5`$ found in LG and is steeper than the $`\alpha 1.0`$

of other previous analyses (Jorissen, Mayor & Udry 2001, Zucker & Mazeh 2002, Tabachnik & Tremaine 2001). This means that (within the same period range) instead of $`M_{Jup}`$ mass exoplanets being twice as common as $`2M_{Jup}`$ exoplanets, we find they are slightly more than three times as common. Similarly we find there are $`3`$ times as many $`0.5M_{Jup}`$ as $`M_{Jup}`$ exoplanets. When the histogram of all 98 (= 101- 3 brown dwarfs) exoplanets is fit, including the highly under-sampled lowest $`Msin(i)`$ bin, the result is $`\alpha =1.2\pm 0.2`$. This should be compared to the $`0.8`$ of earlier fits (LG) and can be found in Marcy et al. 2002. This value of $`\alpha `$ is close to the $`0.8\pm 0.2`$ found for very low mass stars (Bejar et al. 2001). When the lowest exoplanet $`Msin(i)`$ bin is ignored because of known incompleteness we obtain $`\alpha =1.5`$ (LG reported $`1.1`$ for this case). Fitting the $`Msin(i)`$ histogram of the less-biased sample of 52 exoplanets, uncorrected for under-sampling, yields $`\alpha =1.5`$ (LG reported $`1.3`$). After correcting for under-sampling (with a 4 planet correction) we obtain our final result: $`\alpha =1.6\pm 0.2`$ (LG reported $`1.5`$). The 27 exoplanets detected between January and August 2002 fit the trends quantified in LG but also indicate that the slopes are even slightly steeper and thus that Jupiters are slightly more common than indicated by LG.

The biggest uncertainty in this analysis is not the linear approximation used to make the necessary completeness correction – it is probably not knowing how far one can reasonably extrapolate the trends identified. The surface density of the material in the protoplanetary disk available to make planets has to decrease and then drop off at some point. Thus, there is a danger of extrapolating a trend into this region – but where is it? This uncertainty is why we did not extend our analysis to Saturn’s orbit and declare more speculatively that Saturns are typical gaseous planets. However, since the exoplanet data indicates that Jupiters are common and we know of no models in which Jupiters are formed readily but Saturns are not, we see no reason to believe that this eventual drop off is near the region (inner edge of the ice zone between 4 and 10 AU) where abundant material is expected. Thus, since the data indicates that Jupiters are abundant, the most reasonable hypothesis is that Saturns probably are too. Exoplanet data on Saturns ( $`29`$ year orbital period) may be available in a decade or two.

## References

Bejar, V.J.S. et al. 2001, ApJ, 556, 830-836

Jorissen, A., Mayor, M. & Udry, S. 2001, A&A, 379, 992, astro-ph/0105301

Lineweaver, C.H. and Grether, D. 2002, Astrobiology, in press, http://xxx.lanl.gov/abs/astro-ph/0201003

Marcy, G.W. et al. 2002, http://exoplanets.org/science.html

Naef, D. et al. 2001, A&A 375, 205, astro-ph/0196255

Tabachnik, S. & Tremaine, S. 2001, astro-ph/0107482

Zucker, S. & Mazeh, T. 2002, Astrophys. Jour. 562, 1038, astro-ph/0106042


text:
# The NuTeV Anomaly, Neutrino Mixing, and a Heavy Higgs

The NuTeV experiment Zeller:2001hh  has measured the ratios of neutral to charged current events in muon (anti)neutrino – nucleon scattering:

$`R_\nu `$ $`=`$ $`{\displaystyle \frac{\sigma (\nu _\mu N\nu _\mu X)}{\sigma (\nu _\mu N\mu ^{}X)}}=g_L^2+rg_R^2,`$ (1)
$`R_{\overline{\nu }}`$ $`=`$ $`{\displaystyle \frac{\sigma (\overline{\nu }_\mu N\overline{\nu }_\mu X)}{\sigma (\overline{\nu }_\mu N\mu ^+X)}}=g_L^2+{\displaystyle \frac{g_R^2}{r}},`$ (2)
where

$$r=\frac{\sigma (\overline{\nu }_\mu N\mu ^+X)}{\sigma (\nu _\mu N\mu ^{}X)}\frac{1}{2},$$
(3)
and has determined the parameters $`g_L^2`$ and $`g_R^2`$ to be

$`g_L^2`$ $`=`$ $`0.3005\pm 0.0014,`$ (4)
$`g_R^2`$ $`=`$ $`0.0310\pm 0.0011.`$ (5)
The Standard Model (SM) predictions based on a fit to LEP/SLD data are cited as $`[g_L^2]_{\mathrm{SM}}=0.3042`$ and $`[g_R^2]_{\mathrm{SM}}=0.0301`$ in Ref. Zeller:2001hh  and we see a $`2.6\sigma `$ disagreement in $`g_L^2`$. In terms of the ratios $`R_\nu `$ and $`R_{\overline{\nu }}`$, this means that the neutral current events are not as numerous as the SM predicts when compared to the charged current events. On the LEP/SLD side, the measured invisible width is also known to be $`2\sigma `$ below the SM prediction. Both these facts seem to suggest that the $`Z\nu \nu `$ coupling is somehow suppressed.

Suppression of the $`Z\nu \nu `$ coupling can be arranged in a variety of models by mixing the neutrino with a heavy gauge singlet state as discussed in Ref. Davidson:2001ji . For the sake of simplicity, consider the case where all three generations of light neutrinos mix with a universal angle with such states. In such theories, if the $`Z\nu \nu `$ couplings are suppressed by a universal factor of $`(1\epsilon )`$, then the $`W\mathrm{}\nu `$ couplings are also suppressed by a factor of $`(1\epsilon /2)`$. Then the numerators of $`R_\nu `$ and $`R_{\overline{\nu }}`$ will be suppressed over their denominators, so such a mixing could in principle explain the NuTeV anomaly.

However, one must recall that one of the inputs used in calculating the SM predictions is the Fermi constant $`G_F`$ which is extracted from the muon decay constant $`G_\mu `$. The suppression of the $`W\mathrm{}\nu `$ couplings would lead to the correction

$$G_F=G_\mu (1+\epsilon ),$$
(6)
which will affect all SM predictions. The authors of Ref. Davidson:2001ji  conclude that $`\epsilon `$ cannot be chosen to explain the NuTeV data without destroying the excellent agreement between LEP/SLD and the SM.

But what if we include oblique corrections from new physics Peskin:1990zt  to compensate for the shift in $`G_F`$? Note that the Fermi constant always appears multiplied by the $`\rho `$-parameter in neutral current amplitudes so a shift in $`G_F`$ should be absorbable into the $`T`$-parameter.

In Ref. LOTW , we perform a global fit to the NuTeV and LEP/SLD data to determine what values of $`\epsilon `$, $`S`$, and $`T`$ are preferred. Using $`m_H=115\mathrm{GeV}`$, $`m_t=174.3\mathrm{GeV}`$ as the reference SM, we obtain excellent agreement between theory and experiment at

$`S`$ $`=`$ $`0.05\pm 0.10,`$ (7)
$`T`$ $`=`$ $`0.40\pm 0.16,`$ (8)
$`\epsilon `$ $`=`$ $`0.0028\pm 0.0010.`$ (9)
What kind of ‘new physics’ is compatible with these values of $`S`$ and $`T`$? It turns out that the simplest solution is the SM itself with a larger Higgs mass. Thus neutrino mixing together with a heavy Higgs will reconcile the LEP/SLD and NuTeV results.

Details of this analysis will be presented in Ref. LOTW .

Acknowledgments

I would like to thank Will Loinaz, Naotoshi Okamura, and L. C. R. Wijewardhana for their collaboration on this project, and Michio Hashimoto and Mihoko Nojiri for helpful discussions. This research was supported by the U.S. Department of Energy, grant DE–FG05–92ER40709, Task A.


text:
# Lyapunov exponents in Minkowskian U(1) gauge theory Support by FWF P14435-TPH and computer code for Lyapunov exponents by T.S. Biró are gratefully acknowledged.

## 1 Monopole and photon part of U(1)

We begin with a $`4d`$ U(1) gauge theory described by the Euclidean action $`S\{U_l\}=\beta _p(1\mathrm{cos}\theta _p)`$, where $`U_l=U_{x,\mu }=\mathrm{exp}(i\theta _{x,\mu })`$ and $`\theta _p=\theta _{x,\mu }+\theta _{x+\widehat{\mu },\nu }\theta _{x+\widehat{\nu },\mu }\theta _{x,\nu }.`$ We are interested in the relationship between monopoles and classical chaos across the phase transition at $`\beta _c1.01`$. Following Ref. , we have factorized our gauge configurations into monopole and photon fields. The U(1) plaquette angles $`\theta _{x,\mu \nu }`$ are decomposed into the “physical” electromagnetic flux through the plaquette $`\overline{\theta }_{x,\mu \nu }`$ and a number $`m_{x,\mu \nu }`$ of Dirac strings through the plaquette

$$\theta _{x,\mu \nu }=\overline{\theta }_{x,\mu \nu }+2\pi m_{x,\mu \nu },$$
(1)
where $`\overline{\theta }_{x,\mu \nu }(\pi ,+\pi ]`$ and $`m_{x,\mu \nu }0`$ is called a Dirac plaquette.

## 2 Classical chaotic dynamics from quantum Monte Carlo initial states

Chaotic dynamics in general is characterized by the spectrum of Lyapunov exponents. These exponents, if they are positive, reflect an exponential divergence of initially adjacent configurations. In case of symmetries inherent in the Hamiltonian of the system there are corresponding zero values of these exponents. Finally negative exponents belong to irrelevant directions in the phase space: perturbation components in these directions die out exponentially. Pure gauge fields on the lattice show a characteristic Lyapunov spectrum consisting of one third of each kind of exponents . Assuming this general structure of the Lyapunov spectrum we investigate presently its magnitude only, namely the maximal value of the Lyapunov exponent, $`L_{\mathrm{max}}`$.

The general definition of the Lyapunov exponent is based on a distance measure $`d(t)`$ in phase space,

$$L:=\underset{t\mathrm{}}{lim}\underset{d(0)0}{lim}\frac{1}{t}\mathrm{ln}\frac{d(t)}{d(0)}.$$
(2)
In case of conservative dynamics the sum of all Lyapunov exponents is zero according to Liouville’s theorem, $`L_i=0`$. We utilize the gauge invariant distance measure consisting of the local differences of energy densities between two $`3d`$ field configurations on the lattice:

$$d:=\frac{1}{N_P}_P\left|\mathrm{tr}U_P\mathrm{tr}U_P^{}\right|.$$
(3)
Here the symbol $`_P`$ stands for the sum over all $`N_P`$ plaquettes, so this distance is bound in the interval $`(0,2N)`$ for the group SU(N). $`U_P`$ and $`U_P^{}`$ are the familiar plaquette variables, constructed from the basic link variables $`U_{x,i}`$,

$$U_{x,i}=\mathrm{exp}\left(aA_{x,i}^cT^c\right),$$
(4)
located on lattice links pointing from the position $`x=(x_1,x_2,x_3)`$ to $`x+ae_i`$. The generators of the group are $`T^c=ig\tau ^c/2`$ with $`\tau ^c`$ being the Pauli matrices in case of SU(2) and $`A_{x,i}^c`$ is the vector potential. The elementary plaquette variable is constructed for a plaquette with a corner at $`x`$ and lying in the $`ij`$-plane as $`U_{x,ij}=U_{x,i}U_{x+i,j}U_{x+j,i}^{}U_{x,j}^{}`$. It is related to the magnetic field strength $`B_{x,k}^c`$:

$$U_{x,ij}=\mathrm{exp}\left(\epsilon _{ijk}aB_{x,k}^cT^c\right).$$
(5)
The electric field strength $`E_{x,i}^c`$ is related to the canonically conjugate momentum $`P_{x,i}=\dot{U}_{x,i}`$ via

$$E_{x,i}^c=\frac{2a}{g^3}\mathrm{tr}\left(T^c\dot{U}_{x,i}U_{x,i}^{}\right).$$
(6)
The Hamiltonian of the lattice gauge field system can be casted into the form

$$H=\left[\frac{1}{2}P,P+\mathrm{\hspace{0.17em}1}\frac{1}{4}U,V\right].$$
(7)
Here the scalar product stands for $`A,B=\frac{1}{2}\mathrm{tr}(AB^{})`$. The staple variable $`V`$ is a sum of triple products of elementary link variables closing a plaquette with the chosen link $`U`$. This way the Hamiltonian is formally written as a sum over link contributions and $`V`$ plays the role of the classical force acting on the link variable $`U`$.

We prepare the initial field configurations from a standard four dimensional Euclidean Monte Carlo program on a $`12^3\times 4`$ lattice varying the inverse gauge coupling $`\beta `$ . We relate such four dimensional Euclidean lattice field configurations to Minkowskian momenta and fields for the three dimensional Hamiltonian simulation by identifying a fixed time slice of the four dimensional lattice.

## 3 Chaos, confinement and continuum

We start the presentation of our results with a characteristic example of the time evolution of the distance between initially adjacent configurations. An initial state prepared by a standard four dimensional Monte Carlo simulation is evolved according to the classical Hamiltonian dynamics in real time. Afterwards this initial state is rotated locally by group elements which are chosen randomly near to the unity. The time evolution of this slightly rotated configuration is then pursued and finally the distance between these two evolutions is calculated at the corresponding times. A typical exponential rise of this distance followed by a saturation can be inspected in Fig. 1 from an example of U(1) gauge theory in the confinement phase and in the Coulomb phase. While the saturation is an artifact of the compact distance measure of the lattice, the exponential rise (the linear rise of the logarithm) can be used for the determination of the leading Lyapunov exponent. The left plot exhibits that in the confinement phase the original field and its monopole part have similar Lyapunov exponents whereas the photon part has a smaller $`L_{max}`$. The right plot in the Coulomb phase suggests that all slopes and consequently the Lyapunov exponents of all fields decrease substantially.

The main result of the present study is the dependence of the leading Lyapunov exponent $`L_{\mathrm{max}}`$ on the inverse coupling strength $`\beta `$, displayed in Fig. 2. As expected the strong coupling phase is more chaotic. The transition reflects the critical coupling to the Coulomb phase. The plot shows that the monopole fields carry Lyapunov exponents of nearly the same size as the full U(1) fields. The photon fields yield a non-vanishing value in the confinement ascending toward $`\beta =0`$ for randomized fields which indicates that the decomposition works well only around the phase transition.

An interesting result concerning the continuum limit can be viewed from Fig. 3 which shows the energy dependence of the Lyapunov exponents for the U(1) theory and its components. One observes an approximately linear relation for the monopole part while a quadratic relation is suggested for the photon part in the weak coupling regime. From scaling arguments one expects a functional relationship between the Lyapunov exponent and the energy

$$L(a)a^{k1}E^k(a),$$
(8)
with the exponent $`k`$ being crucial for the continuum limit of the classical field theory. A value of $`k<1`$ leads to a divergent Lyapunov exponent, while $`k>1`$ yields a vanishing $`L`$ in the continuum. The case $`k=1`$ is special leading to a finite non-zero Lyapunov exponent. Our analysis of the scaling relation (8) gives evidence, that the classical compact U(1) lattice gauge theory and especially the photon field have $`k2`$ and with $`L(a)0`$ a regular continuum theory. The monopole field signals $`k1`$ and stays chaotic approaching the continuum.