Free energy landscape in a dense hard sphere system

Free energy landscape of a dense hard-sphere system Chandan Dasgupta∗ arXiv:cond-mat/9808142v1 [cond-mat.dis-nn] 13 Aug 1998 Department of Physics, Indian Institute of Science, Bangalore 560012, India Oriol T. Valls School of Physics and Astronomy and Minnesota Supercomputer Institute, University of Minnesota, Minneapolis, Minnesota 55455 (October 9, 2018) The topography of the free energy landscape in phase space of a dense hard sphere system characterized by a discretized free energy functional of the Ramakrishnan-Yussouff form is investigated numerically using a “microcanonical” Monte Carlo procedure. We locate a considerable number of glassy local minima of the free energy and analyze the distributions of the free energy at a minimum and an appropriately defined phase-space “distance” between different minima. We find evidence for the existence of pairs of closely related glassy minima (“two-level systems”). We also investigate the way the system makes transitions as it moves from the basin of attraction of a minimum to that of another one after a start under nonequilibrium conditions. This allows us to determine the effective height of free energy barriers that separate a glassy minimum from the others. The dependence of the height of free energy barriers on the density is investigated in detail. The general appearance of the free energy landscape resembles that of a putting green: relatively deep minima separated by a fairly flat structure. We discuss the connection of our results with the Vogel-Fulcher law and relate our observations to other work on the glass transition. 64.70.Pf, 64.60.Ak, 64.60.Cn I. INTRODUCTION When a liquid is cooled to temperatures below the equilibrium freezing temperature at a sufficiently fast rate to prevent crystallization, it enters a metastable supercooled state. As the temperature is lowered further, the supercooled liquid undergoes a glass transition to a state in which it behaves in most ways like a disordered solid. The dynamic behavior of supercooled liquids near the glass transition exhibits many peculiar features [1–3], such as multi-stage, non-exponential decay of fluctuations and a rapid growth of relaxation times, which are not fully understood theoretically. An intuitively appealing description that is often used [4,5] for qualitative explanations of the observed behavior near the glass transition is based on the so-called “free energy landscape” paradigm. The starting point of this description is a free energy functional that expresses the free energy of a liquid as a functional of the time-averaged local number density. At high temperatures (or at low densities in systems, such as those consisting of hard spheres, where the density is the control parameter), this free energy functional is believed to have only one minimum that represents the uniform liquid state. As the temperature is decreased to values near the equilibrium crystallization temperature (or mutatis mutandi the density is increased), a new minimum representing the crystalline solid, characterized by a periodic modulation of the local density, should also develop. In the “free energy landscape” paradigm, it is assumed that a large number of “glassy” local minima of the free energy, characterized by inhomogeneous but aperiodic density distributions, also come into existence at temperatures close to the equilibrium crystallization temperature. If the system gets trapped in one of these glassy local minima as it is cooled rapidly from a high temperature, crystallization does not occur and the subsequent dynamics of the system is governed by thermally activated transitions among a subset of the large number of metastable glassy minima. If the system visits a large number of these minima during its evolution over a relatively long observation time, it behaves like a liquid over such time scales, in the sense that the time-averaged local density remains uniform. However, the dynamic behavior in this regime, being governed by thermally activated transitions over free energy barriers of varying height, is expected to be slow and complex. In this picture, the glass transition occurs when the time scale of transitions among the glassy minima becomes so long that the system remains confined in a single “valley” of the landscape over experimentally accessible time scales. The general features of the free energy landscape described above are very similar to those found in analytic studies [6–8] of certain generalized spin glass models with infinite-range interactions. Similar features have also been found in recent studies [9,10] of spin models with complicated infinite-range interactions, but no quenched disorder. The equilibrium and dynamic behavior of these mean-field models exhibit a striking similarity with the phenomenology of the glass transition. These results suggest that the free energy landscape paradigm may indeed provide a fitting 1 framework for the development of a theoretical understanding of the behavior of supercooled liquids near the glass transition. The development of such a description would obviously require detailed information about the topography of the free energy landscape of dense supercooled liquids. Since the analytic methods used in the aforementioned studies of mean-field models with infinite-range interactions can not be readily generalized to study physical systems with short-range interactions, investigations of the properties of the free energy landscape of simple glass-forming liquids require the use of appropriate numerical methods. We have carried out a number of numerical studies aimed at elucidating the relation between the dynamic behavior of simple model liquids and the structure of the free energy surface in phase space. These studies, carried out for a dense hard-sphere system, are based on a model free energy functional proposed by Ramakrishnan and Yussouff(RY) [11]. A discretized version of this free energy functional was found [12,13] to exhibit a large number of glassy local minima at densities close to or higher than the value at which equilibrium crystallization occurs. [The control parameter for a hard-sphere system is the dimensionless density n∗ ≡ ρ0 σ 3 , where ρ0 is the average number density in the fluid phase and σ is the hard-sphere diameter; increasing (decreasing) n∗ has the same effect as decreasing (increasing) the temperature of systems for which the temperature is the relevant control parameter.] From numerical studies [14–16] of a set of Langevin equations appropriate for this system, we found that the nature of the dynamics changes qualitatively at a “crossover” density near n∗x = 0.95. The dynamics of a system initially prepared in the uniform liquid state continues to be governed by small fluctuations near the uniform liquid minimum of the free energy as long as the density is lower than this crossover value. For values of n∗ higher than the crossover density, the dynamic behavior is governed by transitions among the glassy minima. The time scales for such transitions were estimated from a standard Monte Carlo (MC) method in Ref. [17] and found to increase rapidly with increasing density. In this paper, we present the results of a numerical study in which a new approach is used for further investigations of the properties of the free energy landscape of a dense hard-sphere system. This study is based on the discretized free energy functional [11] used in our previous work. The development of an understanding of the dynamics of the system in the regime where it is governed by transitions among the glassy minima of the free energy requires information about properties of the free energy landscape such as the number of glassy minima, the distribution of their free energies and overlaps, the heights of the saddle points that connect different glassy minima, and how the system evolves from one minimum to another through these saddle points. One also needs to determine the dependence of these quantities on the average density which, as mentioned above, is the relevant control parameter for the hard-sphere system. In the present study, we have developed and used a “microcanonical” MC procedure to obtain quantitative information about some of these features of the free energy landscape. As described in section II below, this MC procedure enables us to study in detail the process of transition between different glassy minima of the free energy and thus provides valuable information about the topography of the free energy surface in phase space. We have also located a large number of glassy minima of the free energy in the course of this study. This gives us useful information about some of the relevant statistical properties of the collection of glassy minima and the dependence of these properties on the density. The main results obtained from this study are summarized below. By performing a study of the probability of transition from a particular glassy minimum to any other as a function of the free energy increment (the excess free energy measured from that at the original minimum), we find that the value of the free energy at which transitions to other minima begin to occur with a high probability is nearly the same for different glassy minima. This observation suggests that the free energy surface in phase space has a “putting green like” topography in which the glassy minima are like ”holes” of varying depth embedded in a relatively flat background. The total number of glassy minima is a sensitive function of the discretization scale and the sample size. For “commensurate” (as defined below) values of these quantities, which allow the existence of a crystalline minimum (which, when present, is the global minimum of the free energy at high densities), the number of glassy minima is relatively large. Systems with incommensurate values of the discretization scale and the sample size exhibit no crystalline minimum and a substantially smaller number of glassy minima. For this reason, we have carried out all our studies of the statistical properties of glassy minima for a commensurate system. We find that the total number of glassy minima for such a system remains nearly constant as the density is varied in the range 0.94 ≤ n∗ ≤ 1.06. The free energies of the glassy minima are distributed over a wide range between the free energy of the uniform liquid and that of the crystalline solid. The width of this range increases as the density is increased. This observation, together with the result that the number of minima is nearly independent of the density, implies that the number of minima per unit interval of the free energy (the “density of states” of glassy minima) decreases with increasing density. A suitably defined “phase space distance” between two different glassy minima also shows a broad distribution. Our study shows the existence of pairs of glassy minima that differ from each other in the rearrangement of a very small number of particles. The height of the free energy barrier that separates two minima belonging to such a pair is found to be quite small. Such pairs may be identified as “two-level systems” which are believed [18] to exist in all glassy systems. Our study of the probability of transition from a particular glassy minimum to the others as a function of the free energy increment and the MC “time” t allows us to define an effective barrier height that depends rather weakly on 2 t. Some of our results for the dependence of this barrier height on the density have been briefly reported [19] in a recent paper. As described there, we found that the growth of this effective barrier height with increasing density is consistent with a Vogel-Fulcher form [20] appropriate for a hard-sphere system [21]. From our numerical results about how the dependence of the effective barrier height on t changes as the density is increased, we were able to conclude that the growth of the barrier height (and the consequent growth of the relaxation time) is primarily due to entropic effects arising from an increase in the difficulty of finding low free-energy paths (saddle points) that connect one glassy local minimum with the others. Some of the details not included in Ref. [19] are provided in the present paper. We also relate the new results described above with the conclusions reached in Ref. [19]. The rest of this paper is organized as follows. In section II, we define the model system studied and describe the numerical methods used. Section III contains a detailed description of the results obtained in our study. Finally, in section IV, we summarize the main conclusions and discuss them in the context of other related work on the glass transition. II. MODEL AND METHODS In this section we define the model free energy used in our study, and introduce the “microcanonical” MC method that we have developed as a means of studying the topography of the free energy surface of the model in phase space. We also discuss in detail the initial conditions and parameters used. A. The free energy As explained in the Introduction, our system is characterized by a free energy functional F [ρ] which is of the following form [11]: Z F [ρ] = Fl (ρ0 ) + kB T dr{ρ(r) ln(ρ(r)/ρ0 ) − δρ(r)}  Z Z ′ ′ ′ − (1/2) dr dr C(|r − r |)δρ(r)δρ(r ) , (1) where Fl (ρ0 ) is the free energy of the uniform liquid at density ρ0 , and δρ(r) ≡ ρ(r) − ρ0 is the deviation of the density ρ at point r from ρ0 . We take our zero of the free energy at the uniform liquid value, i.e. we set Fl (ρ0 ) equal to zero. In Eq.(1), T is the temperature and the function C(r) the direct pair correlation function [22] of the uniform liquid at density ρ0 , which we express in terms of the dimensionless density n∗ ≡ ρ0 σ 3 by making use of the Percus-Yevick [22] approximation. This approximation is known to be quite accurate if the value of ρ0 is not very high, and should be adequate for all the densities (n∗ ≤ 1.06) considered in this study. It is well-known [22] that the direct pair correlation function of simple model liquids characterized by an isotropic, short-range pair-potential with a strongly repulsive core (such as the Lennard-Jones liquid) is very similar to that of the hard-sphere system at high densities. Therefore, we expect the results obtained from this study to apply, at least qualitatively, to other dense model liquids. To perform the numerical calculations, we discretize our system by introducing a three dimensional cubic lattice of size L3 and mesh constant h in which a discrete set of variables, ρi , i = 1, L3 , are defined as ρi ≡ ρ(ri )h3 , where ρ(ri ) is the density at mesh point i. It is often convenient, in performing and describing the calculations, to deal with a dimensionless, normalized free energy per particle f [ρ] defined as: f [ρ] = βF [ρ]/N (2) where N = ρ0 (Lh)3 = n∗ L3 a3 is the total number of particles in the simulation box, β ≡ 1/(kB T ) and a is the ratio h/σ. B. The microcanonical Monte Carlo method Our main objective in this work is to find an efficient way to investigate the topography of the free energy landscape of the hard-sphere system described by the discretized form of the free energy functional defined in Eq.(1). Basically, what one would like to do is to start the system in a known free energy state (e.g. a glassy local minimum of the free energy), and then investigate the topography of the free energy surface near the starting point by allowing the 3 system to evolve in time, finding out which configurations it subsequently visits and where it ends up. A conventional Metropolis algorithm MC procedure, as performed at lower densities in our previous work [17], is not the most efficient way of doing this: From a computational point of view, a certain amount of computer time is spent at every step of a conventional MC simulation in evaluating the exponential of the free energy change. More important, in a conventional MC simulation carried out at the rather high densities we will consider here, it would take a very long time for the system to move out of the basin of attraction of the minimum in which it is initially placed. This makes a conventional MC study of the process of transitions among free-energy minima prohibitively expensive in the density range we consider. In order to overcome these difficulties of a standard MC simulation, we have devised another procedure which we call the “microcanonical” MC method. The reason for the quotes surrounding the word “microcanonical” is that in our method the quantity that is constrained to be lower than a specified value is not the energy but the free energy F defined in Eq.(1). Specifically, the procedure works as follows: we choose a trial value of what we call the free energy increment, which we denote by ∆F , or alternatively by ∆f if we are dealing with the dimensionless version of Eq.(2). Then, starting with initial conditions which, as discussed below, correspond to a configuration where the free energy is at a local minimum, we perform a MC simulation in which we sweep the sites i of the lattice sequentially. At each step and site, we pick another site j at random from the ones that lie within a distance σ from the site i. We then attempt to change the values of ρi and ρj to p(ρi + ρj ) and (1 − p)(ρi + ρj ), where p is a random number distributed uniformly in [0, 1]. The attempted change is accepted, and this is the crucial point, if and only if the free energy after the change is less than Fmax ≡ F0 + ∆F where F0 is the initial value, that is, the value of the free energy at the minimum where we start the computation. The simulation proceeds up to a maximum “time”, tm , measured in MC steps per site (MCS). In implementing this procedure, they key point is that we perform a sweep over a range of values of ∆F , with the same initial conditions. Obviously, if ∆F is “too small”, meaning that it is smaller than the height of the lowest free energy barrier between the starting minimum and any other minima that may be “nearby”, the system is going to remain in the basin of attraction of the starting minimum. As we increase ∆F , there will eventually be one or more accessible minima, that is, minima that the system can find within a “time” t < tm . These minima are separated from the initial minimum by free energy barriers of height less than ∆F . The system may then move to a region of phase space in the basin of attraction of this new minimum, or one of several newly accessible minima. As ∆F is further increased, additional minima will be made accessible, and furthermore, since additional paths will become available between the initial minimum and the minima already accessible at smaller values of ∆F , these minima may be reached in fewer MC steps. Clearly, if one obtains the information of near which minima the system is, and how long it takes to get there, one can begin to map out the free energy landscape. To find out which basin of attraction the system is in at a certain time t, we save the configurations (i.e. the values of the variables ρi ) at suitable, relatively frequent, time intervals ∆t. These configurations are then used as the inputs in a minimization procedure [12] that determines which basin of attraction the system is in. This is done by finding the local free energy minimum the system moves to when one tries to minimize the free energy by making small changes in the variables ρi in such a way that F is always lowered. The entire procedure, that is, running the simulation up to a certain time tm for a set of values of ∆F , saving the configurations at intervals ∆t, and analyzing them, must obviously be repeated a certain number of times (the “number of runs”) and averaged over. Furthermore, suitable ranges of values for tm , ∆F and ∆t must be fixed at the beginning from theoretical considerations coupled to trial runs. Nevertheless, these steps would also be required for a standard MC procedure and therefore the need for them does not in any way detract from the efficiency of the microcanonical method. We have carried out the numerical procedure outlined above at a number of densities in the range 0.94 ≤ n∗ ≤ 1.06. We did not consider densities lower than 0.94 because our earlier work [15,16] has shown that the dynamics of the system is governed by transitions among glassy local minima only at higher densities. Since, as mentioned above, the Percus-Yevick approximation used for the direct correlation function C(r) appearing in Eq.(1) becomes less accurate at relatively high densities [22], values of n∗ > 1.06 were not considered. C. Initial states and system parameters Our computations were performed for two different sets of the two computational system parameters – the sample size L and the mesh size h. In one case we took these two parameters to be commensurate with a close-packed lattice, and in the other incommensurate. This was done chiefly in order to study the dependence of the structure of the free energy landscape on the commensurability properties of the computational system parameters, as well as on their values. We also considered two different kinds of initial conditions, so that we could investigate the topography of the 4 free energy surface in different regions of phase space. The computationally more intensive part of our simulations was carried out for systems of size L = 15 with periodic boundary conditions and mesh size h = σ/4.6. These values of L and σ are incommensurate with a close-packed lattice, and as a result no crystalline minimum was found for these samples. Two kinds of initial conditions were used for such systems. The first kind is the same as that used in Ref. [17]. These are configurations obtained by first allowing the system to evolve from a uniform initial state under Langevin dynamics [14,15] until its free energy (which, we recall, includes a current dependent term in the Langevin model) reaches zero (indicating the departure of the system from the basin of attraction of the uniform liquid minimum of the free energy), and then using the minimization procedure to reach the minimum whose basin of attraction the system is in at that point. That minimum configuration is then the starting point of the present work. All the minima found this way exhibit glassy structure, as determined by the form of the two-point correlation function (see below) of the local density variables ρi . At higher densities, where the Langevin computation is inappropriate, the minima found at lower densities were scaled up by running the minimization program at the higher density using the lower density configuration (which, of course, is not a minimum at the higher density) as the starting point. The other portion of the computations was performed for systems with L = 12 and h = 0.25σ. These values are commensurate with a close-packed (fcc) structure, so that a crystalline minimum is found at sufficiently high densities. Starting configurations used for simulations carried out for such samples were obtained by using the minimization procedure discussed above on randomly inhomogeneous initial configurations. Out of several glassy minima found this way, we selected a few with structures similar to that of the minima used in simulations of the L = 15 sample. Because of the smaller size of these samples, we were able to explore more extensively several aspects of the problem under consideration. Our computations for the L = 15 sample were carried out for a time range tm = 15000 MCS. Computations for the system with L = 12 were usually carried out to tm = 8000 MCS. In both cases the density range 0.94 ≤ n∗ ≤ 1.06 was covered. For the larger size and longer maximum time, an interval ∆t = 5000 MCS was used, while a closer spacing, ∆t = 2000 MCS, was chosen for L = 12. The structure of a local minimum of the free energy may be characterized by the two-point correlation function g(r) of the frozen local density variables ρi at the minimum. This function is defined as X X fij (r)], (3) ρi ρj fij (r)/[ρ2av g(r) = i>j i>j where ρav ≡ i ρi /L3 is the average value of the ρi at the minimum (values of ρav vary from one glassy minimum to another, but are always slightly higher than ρ0 h3 , the value of ρi at the uniform liquid minimum), and fij (r) = 1 if the separation between mesh points i and j lies between r and r + ∆r (∆r is a suitably chosen bin size), and fij (r) = 0 otherwise. In Fig. 1, we have shown the pair correlation functions for two typical minima used as initial states in our simulations. From the structure of g(r) shown in this figure, it is clear that both these minima are glassy. It is also apparent that the structure of the L = 12 minimum is quite similar to that of the L = 15 one. Other minima used in our simulations have a similar structure. P III. RESULTS In this section, we describe in detail the numerical results obtained from our study, and present our analysis of the numerical data. A. Transitions between minima First, we discuss the qualitative behavior of the system as it evolves in “time” under the microcanonical MC “dynamics” from the initial state at t = 0 to t = tm as described in subsection II B. During the evolution of the system, we monitor the dimensionless free energy βF and the maximum and minimum values of the discretized density variables ρi , i = 1, L3 . The maximum value is useful for detecting possible transitions to the neighborhood of the uniform liquid minimum. If the system fluctuates near one of the inhomogeneous minima of the free energy, then the maximum value of ρi would be much higher than the value (close to ρ0 h3 ) it would have if the system were in the vicinity of the uniform liquid minimum. We find that the system does not move to the neighborhood of the liquid minimum for the values of ∆F considered here. The total free energy is found to remain nearly constant at a value slightly lower than the maximum allowed value, Fmax = F0 + ∆F . In some of the runs, we have also monitored, at much more frequent time intervals, a quantity d(t) that measures the “phase space distance” of the system point at time t from the starting point at t = 0. This quantity is defined as 5 d2 (t) = X [ρi (t) − ρi (0)]2 , (4) i where ρi (0) are the values of the density variables at the minimum from where the simulation is started. By monitoring the time-dependence of this quantity, we obtain useful information about how the system explores the free energy landscape as it evolves in time. We find that if the value of the free energy increment ∆F is small enough so that the system remains confined in the basin of attraction of the original minimum over the duration of the simulation, then the phase space distance d(t) saturates (or continues to increase very slowly) after a rapid initial increase. The value at which d(t) levels off increases as ∆F is increased. For values of ∆F that are sufficiently large for the system to be able to move to the basins of attraction of other minima, the transitions to other basins of attraction are usually (but not always) indicated by sudden increases in the value of d(t). Typical results for the time-dependence of d2 (t) for three different values of ∆F are shown in Fig. 2. The data shown were obtained for a L = 12 system at n∗ = 1.02. For ∆f = 1.0 and ∆f = 1.4, the system was found to remain in the basin of attraction of the initial minimum during the time scale (8000 MCS) of the simulation. In the run with ∆f = 1.9, the system was found to have moved to the basin of attraction of a different minimum at t = 2000 MCS. It moved to the basin of attraction of the crystalline minimum between times t = 2000 MCS and t = 4000 MCS, and stayed there for the remaining part of the run. While any signature of the first transition from the initial minimum to the intermediate one is not clearly visible in the time-dependence of d(t) (possibly due to the overlap of any such signature with the initial rapid increase of d(t)), the subsequent transition to the crystalline minimum is clearly indicated by a rapid rise (and eventual saturation) of d(t). We now turn to our results for the process of transition of the system from the basin of attraction of the initial glassy minimum to the basins of other minima of the free energy. The main quantity that we use to present our analysis of these results is the “critical” value of ∆f , that is, the value of this quantity at which the system begins to find other minima with a high probability. This value, which we will denote as ∆fc (t) (or ∆Fc (t) as the case may be), is algorithmically defined as follows: at every time investigated (i.e. times 5000, 10000 and 15000 MCS for L = 15 and times 2000, 4000, 6000 and 8000 MCS for the L = 12 samples), we test for progressively increasing values of ∆f what is the probability, P (∆f, t), that the system has moved to the basin of attraction of a free energy minimum distinct from the one in which it was started. This probability, which we obtain by averaging over a sufficient number of runs, is (at constant time) obviously zero for very small ∆f and rises toward unity as ∆f increases. At a constant ∆f , it increases somewhat with MC time, as the system is allowed to explore further regions of phase space. We define ∆fc (t) as the value of ∆f for which, at that time, the switching probability reaches 1/2. Note that the quantities P and ∆fc are also functions of n∗ . We find that the number of runs required to obtain ∆fc reliably and reproducibly is about ten to fifteen. The appropriate range of ∆f to be investigated is then determined by the need to hit the middle range of the switching probability values. This requires a bit of trial and error initially, but once the range is found, then running simulations at values of ∆f successively stepped up by increments of 0.05 is sufficient. The whole procedure is illustrated in Fig. 3 and Fig. 4, where we have shown the results for the transition probability as a function of the free energy increment ∆f for a L = 12 minimum at n∗ = 1.02 (Fig. 3), and a L = 15 minimum at n∗ = 0.99 (Fig. 4). The validity of our procedure for estimating the values of ∆fc from the numerical data (the estimated values of ∆fc are indicated in these figures) can be judged from the plots. It is clear from the data shown in these figures that the uncertainty in the estimated values of ∆fc is ∼ 0.05, the spacing between successive values of ∆f used in the simulation. As mentioned above, the determination of the probability of transition as a function of ∆f requires repeating the numerical procedure a number of times for a fixed set of values of ∆f . We find that the minima to which the system moves for values of ∆f close to or higher than ∆fc are, in general, different for different runs. This is more obvious for L = 12 samples which, as discussed below, exhibit a larger number of distinct glassy minima. This observation suggests that ∆fc represents a measure of the free energy increment for which a relatively large region of phase space becomes accessible to the system. Another observation that supports this interpretation is that the system almost never returns to the basin of attraction of the initial minimum after making a transition to the basin of attraction of a different one: after having left the initial minimum, the system cannot find its way back. In a few runs, we found transitions at relatively small values of ∆f which are always to the basin of attraction of the same minimum. In most of these cases, the new minimum was found to be very “close” in phase space (as measured by the quantity defined in Eq. (4)), to the initial one. These are examples of so-called “two-level systems” discussed in more detail in the next subsection. In a few cases, we found that the new minimum to which all the transitions occur at low values of ∆f is not close to the initial one. These are examples of “special” paths with low barrier heights which connect the initial minimum with another specific minimum. Since such transitions and the ones between minima which are very close to each other do not correspond to the opening up of large regions of phase space, we did not include such transitions in the calculation of the transition probability. In this way, we determined ∆fc as a function of the density n∗ and Monte Carlo time t for the ranges of time and density mentioned above. Some of our results for ∆fc are shown in Fig. 5 for a L = 12 minimum and four different 6 values of n∗ , and in Fig. 6 for a L = 15 minimum and also four values of n∗ . We find, as implied above, that ∆fc is a weak function of time, and also of course a stronger function of n∗ . The dependence of ∆fc on t for a fixed value of n∗ becomes more pronounced as n∗ is increased. These dependences were studied in detail in Ref. [19], where analytic fitted forms are given. Part of this analysis and some of the conclusions drawn from it are summarized in subsection III C for completeness. Finally, we note that although our model and numerical procedure are different from those used in most existing numerical studies of dense liquids (such studies use conventional MC or molecular dynamics to simulate the behavior of models defined by a microscopic Hamiltonian), some of the general features found in existing simulations (and also in experiments) are reproduced in our work. We find that if the value of ∆F is such that βFmax exceeds an upper threshold, then the system moves within a few hundred MC steps to the vicinity of the uniform liquid minimum. The value of this upper threshold is found to be close to βF = 5.0. This is the microcanonical analog of the melting transition. As mentioned above, this threshold value is not crossed in the simulations from which the results described here were obtained. We also find that in runs with βFmax close to, but lower than the upper threshold, the system moves to the basin of attraction of the crystalline minimum (for L = 12) with a high probability. This is nothing but the process of annealing: it is well-known from experiments and simulations that crystallization may be induced by heating a glassy system to a temperature close to (but lower than) its melting temperature and then cooling it down. B. Properties of glassy minima In the course of our computations, we have located many of the glassy minima of the free energy. As mentioned above, for the “incommensurate” L = 15 sample used in our work, the number of minima we have located at each density is not large. The total number of minima found for this sample varies in the range of four to six, with some tendency to higher values in the lower part of the n∗ range considered here. The “commensurate” L = 12 sample exhibits, as we shall see below, a substantially larger number of minima, one of which is crystalline (fcc). For this reason, we consider chiefly the results obtained for L = 12, for which we can produce significant statistics, in this subsection. While studying the process of transitions among the minima, we carried out a large number of minimization runs with many different initial states. The total number of such runs is of the order of 103 for each of the values of n∗ studied. While our procedure does not correspond to an exhaustive search for all the local minima of the system, the fairly large number of initial states considered for each value of n∗ ensures that, for L = 12 at least, we did locate a large fraction of the local free energy minima of the system. So, the statistical information obtained from our study can be expected to be representative of the full collection of local minima. The total number of local minima of the L = 12 system remains nearly constant as the density is varied in the range 0.96 ≤ n∗ ≤ 1.06. This number is close to 25. The numbers for different values of n∗ show small variations, but there is no clear systematic trend in the dependence of this number on the density. In most cases, a minimum found at a particular density may be “followed” to higher or lower densities by using the values of ρi at the minimum at the first density as inputs to the minimization program at the new density. In a few cases, we find that a minimum disappears as the density is increased or lowered, but such occurrences are rare. From these observations, we conclude that the total number of glassy minima does not exhibit any strong dependence on the density. Our limited investigation of the variation of the free energies of the glassy minima with density suggests that the ordering of the free energies remains the same (i.e. free energies of different minima do not cross) as the density is changed. The free energies of these minima are distributed in a band that lies between the free energy of the uniform liquid (which, we recall, is taken to be the zero of the free energy scale) and that of the crystalline solid. The width of this band increases with increasing n∗ . Since the number of minima is approximately independent of the density, this implies that the “density of states” of the glassy minima decreases as n∗ is increased. Specifically, let p(βF )δ be the probability of finding a glassy minimum with dimensionless free energy between βF − δ/2 and βF + δ/2. We have calculated this quantity from our data at different values of n∗ . Representative results at two densities, n∗ = 0.96 and n∗ = 1.02, are shown in Fig. 7. The values of δ used are 4.0 and 8.0 for n∗ = 0.96 and n∗ = 1.02, respectively. While the distributions for the two densities are qualitatively similar, the range of βF over which p(βF ) is nonzero is clearly wider at the higher density. The consequent decrease in the values of p(βF ) with increasing density is also clearly seen. Both distributions show peaks near the upper end, and tails extending to substantially lower values. However, the lowest free energy of the glassy minima is substantially higher than the free energy of the crystalline minimum (for the crystalline minimum, βF = −102.4 for n∗ = 0.96 and βF = −167.4 for n∗ = 1.02). If the probability of finding the system in a glassy minimum is assumed to be proportional to the Boltzmann factor e−βF , then only those minima with free energies lying near the lower end of the band would be relevant in determining the equilibrium and dynamic properties of the system. Our results indicate that the number of such “relevant” minima decreases with 7 increasing n∗ . In the present study, we find certain correlations between the free energy of a glassy minimum and its structure. Similar correlations were also found and described in some detail in Ref. [15]. Specifically, we find that minima with lower free energies have more “structure” (as indicated by e.g. the heights of the first and second peaks of the correlation function g(r) defined in Eq.(3)) and higher density that those with lower free energies. We have also studied how the distributions of the local density variables in two distinct glassy minima differ from one another. To do this, we need a measure of the difference between the distributions of ρi in two glassy minima. This measure should satisfy the requirement that it yield a zero value for the difference between two configurations if one of them can be mapped to the other by a symmetry operation of the computational mesh. The symmetries of the cubic mesh used in our computation include the 48 symmetry operations of a simple cubic lattice and all translations, taking into account the periodic boundary conditions. The quantity dm (1, 2) that we have used to measure the difference in the density distributions at two minima labeled “1” and “2” is defined as follows: X (1) 1 (2) (5) [ρi − ρR(i) ]2 , dm (1, 2) = min{R} 2 i (1) (2) where ρi and ρi are the discretized densities at two minima, R represents one of the symmetry operations mentioned above, R(i) is the mesh point to which mesh point i is transformed under R, and min{R} means that the R that minimizes the quantity on the right is to be taken. Since the variables ρi in an inhomogeneous minimum are close to one at the mesh points corresponding to the locations of the “particles”, and close to zero at the other mesh points, the quantity dm basically measures the number of particles whose positions are different in the two minima being compared. In Fig. 8, we display in histogram form the results for the distribution of dm at two values of the density. The two distributions are qualitatively similar. Both are small at small values of dm and exhibit peaks near dm = 15, which corresponds to about half of the total number of particles having different locations in the two minima. From these results, we conclude that most of the glassy minima are rather different from one another. The arrangement of the particles in the glassy minima is also very different from that in a crystalline minimum, as indicated by the observation that the value of dm almost always lies above 15 if one of the two minima being compared is glassy and the other one is crystalline. The degree of similarity between two different minima may also be quantified in terms of their “overlap” [8]. For the discretized system considered here, the dimensionless overlap q(1, 2) between two minima labeled “1” and “2” may be defined in the following way: X (1) 1 (2) q(1, 2) = (6) [ρi − ρav ][ρR(i) − ρav ], max{R} ρav L3 i where ρav is the average value of the ρi , which is assumed to be the same in the two minima, and max{R} means that the R that maximizes the quantity on the right is to be taken. Using the aforementioned fact that at the glassy minima found in the density range considered here, the values of ρi are close to one at a small number of mesh points and close to zero at others, the following approximate relation between q and dm may be derived easily: q(1, 2) ≃ 1 − dm (1, 2)/N − ρav , (7) where N ≡ ρav L3 is the total number of particles in the simulation box. The observation that the distribution of dm has a peak near dm = N/2 then implies that the distribution of q peaks near the value 0.5. The distributions shown in Fig.8 extend down to values of dm as small as 2 or 3, indicating that there are a few pairs of glassy minima which are very similar to each other. For each value of n∗ , we find a small number (3-5) of such pairs of minima. To take an example, for n∗ = 0.96, we have found two minima, with free energies βF = −47.7 and −47.0, for which the value of dm is 2.4. A detailed examination of the density distributions at these two minima reveals that the main difference between their structures comes from small displacements of just two particles. Of course, these displacements also produce small changes in the values of ρi at neighboring mesh points. We believe that these pairs of minima are examples of “two-level systems” whose existence in glassy materials was postulated [18] many years ago in order to account for some of the experimentally observed low-temperature properties. The height of the free energy barrier that separates two members of a two-level system is expected to be low. Our observations are consistent with this expectation. For the pair of minima mentioned above, we find that if we start the system from the minimum with βF0 = −47.7 and carry out our numerical procedure for finding transitions to other minima, the system begins to show transitions to the minimum with βF0 = −47.0 as the value of ∆f is increased above 0.7. For 0.7 ≤ ∆f ≤ 1.4, all the transitions are to the other member of the two-level system. Transitions to other minima begin to appear only for higher values of ∆f . (As noted above, we did not include transitions between the members of a two-level system in our calculation of ∆fc .) 8 We have also looked at how the quantity Fc = F0 + ∆Fc , which measures the value of the total free energy at which transitions to other minima begin to occur with a high probability, varies from one minimum to another. As exemplified by Fig.9, where we present the results for βFc for four minima at n∗ = 0.96 and for three minima at n∗ = 1.02, the value of this quantity is nearly constant for each value of n∗ . While the values of βF0 vary over a range of about 15 at n∗ = 0.96, and over a range of about 40 at n∗ = 1.02, the calculated values of βFc are nearly the same (within the error bars) for the different minima at both densities. This observation suggest a “putting green like” free energy landscape in which the local minima are like “holes” of varying depth in a nearly flat background. This structure also implies that there is a strong correlation between the depth of a minimum and the height of the barriers that separate it from the other minima: the barriers are higher for deeper minima. C. The Vogel-Fulcher law and entropic effects As first pointed out in Ref. [19], there is a direct connection between our results and the Vogel-Fulcher law [20]. We summarize this connection here. The basic point is that the results for ∆fc as shown in Figs. 5 and 6 are consistent with the form: ∆fc (n∗ , t) = a(t) + b, n∗c − n∗ (8) where a(t) is a weak function of t, b is a constant, and the “critical” density n∗c is found to be independent of t within the accuracy of our results. In Fig.10, we show the data for ∆fc for a L = 15 minimum at times 5000, 10000 and 15000 MCS, and also the best fits of the data to the form of Eq.(8) with b = 0. The parameter values for the best fits are: a = 0.31, n∗c = 1.19 for t = 5000; a = 0.30, n∗c = 1.22 for t = 10000; a = 0.27, n∗c = 1.23 for t = 15000. The form of Eq.(8) leads at once to the Vogel-Fulcher law appropriate for a hard sphere system [21] since the characteristic time should be proportional to the exponential of β∆Fc . The values of n∗c obtained from the fits, particularly at later times, are very close to the random close packing density, n∗rcp ≃ 1.23. This is in agreement with the results of molecular dynamics simulations [21]. The L = 12 data yield similar values of n∗c , but with b ≃ 1.0. The weak dependence of ∆fc on t was also analyzed in detail in Ref. [19] where it was found that this dependence for all values of n∗ and all the minima studied is well-described by the form: ∆fc (n∗ , t) = c(n∗ )t−α + ∆f0 , (9) with α in the range 0.24 − 0.40, and ∆f0 nearly independent of n∗ . The coefficient c(n∗ ) was found to increase with increasing n∗ . Fits of the data to the form of Eq.(9) are shown in Fig. 2 of Ref. [19]. [This form agrees with Eq.(8) if a(t) ∝ t−α and c(n∗ ) ∝ 1/(n∗c − n∗ ); our data are consistent with these conditions.] This result suggests a physical interpretation of the observed Vogel-Fulcher behavior. The quantity ∆f0 (the value of ∆fc in the t → ∞ limit) provides a measure of (β/N times) the height of the lowest free energy barriers that must be crossed in order to reach some of the other local minima of the free energy from the one under consideration. As discussed in detail in Ref. [19], the coefficient c(n∗ ) may be interpreted as a measure of the difficulty of finding low free-energy paths to other minima. The observation that c(n∗ ) increases with n∗ while ∆f0 is nearly independent of n∗ then implies that the increase of the effective barrier height with increasing n∗ is primarily due to “entropic” effects associated with the difficulty of finding low-lying saddle points that connect a minimum with the others. Therefore, a picture emerges from our work as to the origin of the Vogel-Fulcher divergence. As the system evolves over longer and longer times, the probability that it will find paths to other minima involving jumps over lower and lower free energy barriers increases. At early times, the system can only explore nearby paths and must then jump over whatever barrier is available in that region. At longer times, a wider region is explored and the chances of finding a path with a lower barrier increase. What our argument shows is that the Vogel-Fulcher law follows from the fact that the difficulty of finding such low-free-energy paths to other minima increases with increasing n∗ . IV. SUMMARY AND DISCUSSION We have developed and used in this work a numerical method to study the topography of the free energy surface of a dense hard-sphere system characterized by a model free energy functional. At the relatively high densities considered in this study, this system exhibits a complex free energy landscape characterized by the presence of many glassy local minima. The number of accessed glassy local minima is found to depend strongly on the commensurability properties of the discretization scale h and the sample size L used. For fixed values of these parameters, the number of minima 9 is nearly independent of the density in the range studied. In the case where L and h are commensurate, a crystalline minimum is found and the number of glassy minima accessed is large enough to allow for statistical study. The free energy values at its minima are distributed over a broad band whose width increases with increasing density. The phase-space distance between different minima shows a broad distribution with a peak near the high end. However, there are a few pairs of minima which are very close to each other in phase space and are separated by low freeenergy barriers. These, we believe, are examples of “two-level systems” which are expected to be present in all glassy materials. We have found in all cases a strong correlation between the depth of a minimum and the effective height of free energy barriers that separate it from the other minima: deeper minima have higher barriers. The observed density-dependence of the effective barrier height is consistent with the Vogel-Fulcher law. Our results indicate that this Vogel-Fulcher growth is primarily due to an increase in the difficulty of finding low-free-energy paths to other minima as the density is increased. Our results have close connections with those of a number of recent studies of the equilibrium and dynamic properties of dense supercooled liquids. We first discuss the relation of our observations with spin-glass-like theories [7,8,23] of the structural glass transition. These theories are based on the similarity between the phenomenology of the structural glass transition in so-called “fragile” [3] liquids and the behavior found in a class of generalized mean-field spin glass models [6,24] with infinite-range interactions, and also in certain mean-field spin models with complicated multispin interactions but no quenched disorder [9,10]. At high temperatures, the free energy of these mean-field models, expressed as a function of the single-site magnetizations, exhibits only one minimum – the “paramagnetic” one at which all the site magnetizations are zero. As the temperature is lowered, an exponentially large number of non-trivial local minima come into existence at a characteristic temperature Td . A “dynamic transition”, characterized by a breaking of ergodicity, occurs at this temperature. At this “transition”, the system gets trapped in the basin of attraction of one of the newly developed local minima and remains confined in this basin for all subsequent times because the free energy barriers between different local minima diverge in the thermodynamic limit in these models. This “dynamic transition” does not have any signature in the equilibrium behavior of the system. A thermodynamic phase transition occurs at a lower temperature Tc at which the configurational entropy associated with the exponentially large number of free-energy minima becomes non-extensive. In the suggested analogy between these models and the structural glass transition, the paramagnetic minimum of the free energy is identified with the one that represents the uniform liquid, and the role of the non-trivial local minima of the free energy is played by the glassy local minima of the liquid free energy. The “dynamic transition” found at Td in the mean-field spin models is thought to be smeared out in liquids. This is because the free energy barriers between different minima are expected to remain finite in physical systems with finite-range interactions. It has been suggested [7,8,23] that the temperature Td should be identified with the “ideal glass transition” temperature of mode-coupling theories [25] of the dynamics of dense liquids. This temperature is supposed to signal the onset of activated processes in the dynamics. The temperature Tc is interpreted as the “Kauzmann temperature” [26] at which the difference in entropy between the supercooled liquid and the crystalline solid extrapolates to zero. The relaxation time of the supercooled liquid is supposed to diverge at the same temperature. Heuristic arguments that suggest that this divergence is of the Vogel-Fulcher form have been proposed [8,23]. These arguments are based on an entropic mechanism associated with the vanishing of the configurational entropy at Tc . The behavior found in our numerical study is in qualitative agreement with this scenario. We find a characteristic density (we recall once more that the density plays the role of the temperature in the hard-sphere system we consider) at which a large number of glassy minima of the free energy come into existence. We do not yet know whether the number of glassy minima depends exponentially on the sample volume – a study of this question is difficult due to the dependence of the number of minima on the commensurability of h and L. While the number of glassy minima for fixed h and L remains nearly constant as the density is increased, the configurational entropy associated with these minima decreases with increasing density because the width of the band over which the free energy of these minima is distributed increases with density. As discussed above, we have also found evidence for a Vogel-Fulcher-type growth of relaxation times driven by an entropic mechanism. There are, however, a number of differences between the details of our findings and the predictions of spin-glass-like theories. In our earlier work [15,16] on the Langevin dynamics of the model system considered here, we found that the dynamic behavior is governed by activated processes if the dimensionless density n∗ exceeds a crossover value, n∗x , of about 0.95. This value is substantially higher than the value of n∗ (≈ 0.8) at which the glassy minima make their first appearance. These two densities are assumed to be the same in spin-glass-like theories. Another difference lies in the values of the free energy of the glassy minima relative to that of the uniform liquid. We find that the free energy of a glassy minimum becomes lower than that of the uniform liquid minimum as the density is increased above a value that is only slightly higher than the density at which the glassy minimum comes into existence. In particular, the free energies of the glassy minima are substantially lower than that of the uniform liquid one for values of n∗ near n∗x . This is different from the behavior found in the spin glass models. In these systems, the free energies of the non-trivial local minima remain higher than that of the paramagnetic one over the entire temperature range Tc < T < Td . Our 10 results for the distribution of the overlap between different minima are also somewhat different from those for the spin glass models. We cannot rule out that some of these differences arise from finite-size effects which may be significant for the rather small samples considered in our study. Another possibility is that these differences arise in our system from the effects of small fluctuations about a local minimum, which are unimportant in models with infinite-range interactions. A careful investigation of these issues would be very interesting. A number of numerical studies of “aging” phenomena in the non-equilibrium dynamics of simple model liquids have been reported recently [27–29]. In these studies, the system is quenched from a relatively high temperature to a temperature lower than the numerically determined glass transition temperature, and is then allowed to evolve at this low temperature for a certain “waiting time” tw . Then, the two-time correlation function C(tw , tw + t) of an appropriate fluctuating quantity is measured and the dependence of this correlation function on t and tw is analyzed. The simulations show that the decay of C(tw , tw + t) as a function of t becomes slower as tw is increased. Our results about the topography of the free energy landscape provide a qualitative explanation of this observation. When the system is rapidly quenched to a low temperature (or compressed to a high density in our hard-sphere system), it is likely to get trapped in the basin of attraction of one of the glassy local minima that are close in phase space to the initial configuration. Such minima would not, in general, have the lowest free energies. As the system evolves during the waiting time tw , it can be expected to move progressively to the basins of attraction of minima with lower free energies because such minima would have a higher Boltzmann weight. Since the effective barrier height is higher for deeper minima (this follows from the “putting green like” structure of the free energy landscape), the time scale for subsequent relaxation is expected to increase with increasing tw . This is precisely the behavior found in the aging simulations mentioned above. Our study is rather similar in spirit to numerical investigations of the “potential energy landscape” [30–36] of model liquids characterized by simple Hamiltonians. In such studies, a numerical minimization procedure (e.g. the conjugate gradient method) is used to find the local minima of the total potential energy of small samples as a function of the coordinates of the particles. The potential energy function is generally found to exhibit a large number of local minima. These local minima and the potential energy barriers that separate them define a complex “potential energy landscape”. Recently, there have been several attempts [32,34,35] to relate the properties of this landscape to the dynamic behavior of the liquid. While the similarities between these investigations and our work are obvious, there are several important differences between these two approaches, some of which we now discuss. A study of the potential energy landscape is based on a microscopic Hamiltonian defined in terms of the coordinates of the particles, whereas our work involves a model free energy which is a functional of a coarse grained (both in space and time) density field. Information about the microscopic interactions is incorporated in our description through the direct pair correlation function C(r) appearing in Eq.(1). The free energy of a thermal system is, of course, equal to the potential energy at zero temperature. Therefore, the potential energy landscape of such systems becomes identical to the free energy landscape at T = 0. There are some mean-field spin-glass models (e.g. the p-spin spherical spin glass [37]) in which the correspondence between the local minima of the energy and the free energy extends also to non-zero temperatures. Such a correspondence is not likely to be generic, however. A description based on the potential energy landscape is certainly appropriate at low temperatures where entropic effects are relatively unimportant. But it would, in general, be difficult to extend such a description to higher temperatures where entropic effects play a crucial role. In particular, information about the energy landscape alone would not be sufficient to describe the behavior near a phase transition (such as the melting transition of a solid and the order-disorder transition in magnetic systems) driven by a competition between energetic and entropic effects. In contrast, a description based on a model free energy that includes entropic contributions provides a convenient and intuitively appealing starting point for studying the behavior near such phase transitions. For example, the free energy functional used in our work is known [11] to provide a correct description of the crystallization transition of simple liquids. Another well-known example is the Curie-Weiss theory of magnetism (our approach is analogous to an inhomogeneous version of the Curie-Weiss theory). For these reasons, we believe that our free-energy-based approach is more suitable for a description of the behavior of liquids near the glass transition than approaches based on the potential energy landscape. Another important difference between free energy and potential energy landscapes is that the former changes as the appropriate control parameter (density or temperature) is changed, whereas the latter, being determined completely by the Hamiltonian of the system, remains unchanged. Specifically, some of the local minima of the free energy may appear or disappear as the control parameter is varied (for example, the inhomogeneous minima of the free energy used in our study disappear at sufficiently high densities). Also, the heights of free energy barriers between different local minima may change with the control parameter (see, e.g. Fig.10 where the dependence of the height of a typical free energy barrier on the density is shown). In contrast, the potential energy landscape does not show any such variation as the temperature is changed. This difference may be important in understanding some of the results found in recent studies [35,36] based on the energy landscape of simple model liquids. In Ref. [35], an approximate description of the dynamics of a Lennard-Jones system in supercooled and glassy regimes is developed in terms of the numerically determined properties of the local minima of the potential energy and the energy barriers between them. While this 11 description is found to reproduce several interesting features of glassy dynamics, it does not show the expected faster than Arrhenius growth of the viscosity at low temperatures. This may be due to the energy barriers in this description not changing with temperature. It is possible that a free-energy-based description in which the barrier heights change with temperature would lead to a faster than Arrhenius growth of the viscosity. This possibility is clearly illustrated in our study which shows that the dependence of the heights of free energy barriers on the appropriate control parameter leads to Vogel-Fulcher behavior. Ref. [36] describes a numerical study of the local minima and the saddle points of the potential energy surface of small Lennard-Jones clusters. One of the quantities calculated in this paper is an “entropic ratio” R that approximately quantifies the entropic effects on the rate of thermally activated transitions between two local minima of the potential energy function. Values of R > 1 indicate entropic suppression of the transition rate, whereas R < 1 corresponds to an enhancement. The probability of R having values greater than one is found to be small. This result is interpreted as evidence for entropic effects being relatively unimportant. In particular, the authors mention that this observation contradicts our conclusion (described in detail in Ref. [19] and summarized in section III C above) that the growth of relaxation times in simple glassy liquids is primarily entropic in origin. In our opinion, the results reported in Ref. [36] do not necessarily contradict our conclusion. The difference between our conclusion and that of Ref. [36] about the importance of entropic effects is probably just a reflection of the aforementioned fact that the free energy landscape changes with the control parameter, but the potential energy landscape does not. We find in our free-energy-based study that the growth of the height of a typical effective free energy barrier with increasing density is primarily due to entropic effects arising from an increase in the difficulty of finding low-free-energy paths to other minima. This effect is closely related to changes in the topography of the free energy landscape as the density is changed. In contrast, the potential energy landscape studied in Ref. [36] does not depend on the temperature which is the appropriate control parameter for the Lennard-Jones system considered there. In particular, the calculated values of the height of the potential energy barrier between two minima and the entropic factor R do not change as the temperature is changed. Therefore, there is no direct connection between our results (which, as explained above, are about the changes of these quantities as the appropriate control parameter is changed) and those reported in Ref. [36]. We end with a word of caution. Due to the computational complexity of numerical studies of free energy and potential energy landscapes, such studies have been restricted to rather small samples which may exhibit finitesize effects. In our work, we found that certain features of the free energy landscape are quite sensitive to the commensurability properties of the discretization scale and the sample size. Strong dependences on the sample size and the boundary condition have also been found [33,36] in studies of the potential energy landscape. One should, therefore, be careful in extrapolating the results obtained from such studies to the thermodynamic limit. V. ACKNOWLEDGMENTS A part of the numerical work was carried out at the Supercomputer Education and Research Centre of Indian Institute of Science. One of us (CD) acknowledges support from the Theoretical Physics Institute, University of Minnesota, for a visit. ∗ [1] [2] [3] [4] [5] [6] [7] [8] [9] Also at the Condensed Matter Theory Unit, Jawaharlal Nehru Centre for Advanced Scientific Research, Bangalore 560064, India. For a review, see Liquids, Freezing and the Glass Transition, edited by J. P. Hansen, D. Levesque, and J. Zinn-Justin, (Elsevier, New York, 1991). J. Jäckle, Rep. Prog. Phys. 49, 171 (1986). C. A. Angell, J. Phys. Chem. Solids, 49, 863 (1988). P. W. Anderson in Ill Condensed Matter, Lecture Notes of the Les Houches Summer School, 1978, edited by R. Balian, R. Maynard and G. Toulouse (North Holland, Amsterdam, 1979). P. G. Wolynes in Proceedings of International Symposium on Frontiers in Science, (AIP Conf. Proc. No. 180), edited by S. S. Chen and P. G. Debrunner (AIP, New York, 1988). T. R. Kirkpatrick and P. G. Wolynes, Phys. Rev. A 35, 3072 (1987); Phys. Rev B 36, 8552 (1987). T. R. Kirkpatrick and D. Thirumalai, J. Phys. A 22, L149 (1989). T. R. Kirkpatrick, D. Thirumalai and P. G. Wolynes, Phys. Rev. A 40, 1045 (1989). J. P. Bouchaud and M. Mezard, J. Phys. I (France) 4, 1109 (1994). 12 [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] L. F. Cugliandolo, J. Kurchan, G. Parisi and F. Ritort, Phys. Rev. Lett. 74, 1012 (1995). T. V. Ramakrishnan and M. Yussouff, Phys. Rev. B 19, 2775, (1979). C. Dasgupta, Europhys. Lett. 20, 131 (1992). C. Dasgupta and S. Ramaswamy, Physica A 186, 314 (1992). L. M. Lust, O. T. Valls, and C. Dasgupta, Phys. Rev. E 48, 1787 (1993). C. Dasgupta and O. T. Valls, Phys Rev. E 50, 3916 (1994). O. T. Valls and C. Dasgupta, Transport Theory and Stat. Physics 24, 1199 (1995). C. Dasgupta and O. T. Valls, Phys. Rev. E 53, 2603 (1996). P. W. Anderson, B. I. Halperin and C. M. Varma, Philos. Mag. 25, 1 (1972); W. A. Phillips, J. Low. Temp. Phys. 7, 351 (1972). C. Dasgupta and O. T. Valls, Phys. Rev. E 58, 801 (1998). H. Vogel, Z. Phys. 22, 645 (1921); G. S. Fulcher J. Amer. Ceram. Soc. 8, 339 (1925). L. V. Woodcock and C. A. Angell, Phys. Rev. Lett. 47, 1129 (1981). J. P. Hansen and I. R. McDonald, Theory of Simple Liquids (Academic, London, 1986). G. Parisi in Complex Behaviour of Glassy Systems: Proceedings of the XIV Sitges Conference, edited by M. Rubi and C. Perez-Vicente (Springer, Berlin, 1997). T. R. Kirkpatrick and D. Thirumalai, Phys. Rev. B 36, 5388 (1987). See, for example, W. Götze in Ref. [1]. W. Kauzmann, Chem. Rev. 48, 219 (1948). W. Kob and J.-L. Barrat, Phys. Rev. Lett. 78, 4581 (1997). J.-L. Barrat and W. Kob, 1998 preprint (cond-mat/9806027). G. Parisi, J. Phys. A 30, L765 (1997); Phys. Rev. Lett 79, 3660 (1997). F. H. Stillinger and T. A. Weber, Phys. Rev. A 28, 2408 (1983). T. A. Weber and F. H. Stillinger, Phys. Rev. B 32, 5402 (1985). F. H. Stillinger, Science 267, 1935 (1995). A. Heuer, Phys. Rev. Lett. 78, 4051 (1997). S. Sastry, P. G. Debenedetti, and F. H. Stillinger, Nature 393, 554 (1998). L. Angelani, G. Parisi, G. Ruocco, and G. Villani, 1998 preprint (cond-mat/9803165). G. Daldoss, O. Pilla, G. Villani, and G. Ruocco, 1998 preprint (cond-mat/9804113). J. Kurchan, G. Parisi, and M. A. Virasoro, J. Phys. I (France) 3, 1819 (1993). FIG. 1. The density correlation function g(r), as defined in Eq.(3), plotted as a function of distance (in hard sphere units) for two typical initial configurations. Note the glassy character of the correlations for both lattice sizes. FIG. 2. The quantity d2 (t), which characterizes the “phase space distance” between two points, (see Eq.(4)) plotted as a function of Monte Carlo time for three values of the free energy increment. The regions of sharp changes in the curves are discussed in the text. FIG. 3. Example of the determination of the “critical” value ∆fc , defined as the value of the free energy increment ∆f at which the transition probability P is 1/2. The black dots mark the intersections of the plots with the line P = 0.5 (see text for a complete discussion). The data shown are for a sample of size L = 12. FIG. 4. Another example of the determination of ∆fc (see the caption of Fig.3 for the details), but for a sample of size L = 15. FIG. 5. Results for the critical value, ∆fc , of the free energy increment, as a function of the Monte Carlo time t for four different values of the dimensionless density n∗ . The results shown are for a L = 12 minimum. FIG. 6. Results for ∆fc as a function of the Monte Carlo time t at four different densities for a L = 15 minimum. In this figure and in Fig.5, results for other times and densities studied interpolate smoothly with the results shown. FIG. 7. The “density of states” for glassy free energy minima, defined as the probability of finding a glassy minimum with free energy in a given range (see text). Results for L = 12 samples at two densities are shown. 13 FIG. 8. Histogram representing the fraction of pairs of free energy minima found to differ in their real space density configurations by an amount dm as defined in Eq.(5). Results for L = 12 samples are shown at two different densities. FIG. 9. The value Fc of the free energy at which transitions to other minima begin to occur with a high probability, plotted as a function of F0 , the free energy at the starting minimum. Results are shown for L = 12 at two densities. One can see that, at a given n∗ , the dependence of Fc on the starting value is quite weak. FIG. 10. Vogel-Fulcher fits of the data for ∆fc obtained for a L = 15 minimum at three different values (5000, 10000 and 15000 MCS) of t. The continuous lines are the best fits of the data to the form of Eq.(8) with b = 0. The parameter values for the best fits are given in the text. 14