# The thermodynamical liquid-glass transition in a Lennard-Jones binary mixture

\[

## Abstract

We use the results derived in the framework of the replica approach to study the liquid-glass thermodynamic transition. The main results are derived without using replicas and applied to the study of the Lennard-Jones binary mixture introduced by Kob and Andersen. We find that there is a phase transition due to the entropy crisis. We compute both analytically and numerically the value of the phase transition point $`T_K`$ and the specific heat in the low temperature phase.

\]

In recent times there have been many progresses in the analytic understanding of the thermodynamics transition of glasses . Beyond technicalities the basic assumption is that a glass is very near to a frozen liquid. This quite old statement can be rephrased as a liquid is very near to a heated glass. In other words the configurations of a glass at low temperature are not far from those of a liquid. Therefore, if we use some smart method to explore the phase space in the liquid phase, we can find the properties of the low temperature glassy phase.

This strategy has been put into action using the replica theory. One finds, as output of an explicit computation, that the glass transition is characterized by the vanishing of the configurational entropy (the so called complexity) and the low temperature phase is described by one step replica symmetry breaking  with a non-vanishing non-ergodic parameter at the phase transition point.

Replica theory is a very powerful tool, but it has the disadvantage that many of the underlying physical hypothesis cannot be seen in a clear way. We will rederive some of the main result of  without using the replica formalism . We suppose that below some temperature (to be identified with the mode coupling transition ), the phase space of the system can be approximately divided into regions (which we will call valleys) which are separated by high barriers . This approximation becomes better and better when we decrease the temperature and it becomes exact below the thermodynamics phase transition temperature $`T_K`$ of the glass , where the viscosity should diverges  (unfortunately no quantitative predictions have been obtained microscopically on the behaviour of the viscosity beyond the Adams Gibbs argument). In a first approximation each valley can be associated to one inherent structure, i.e. one minimum of the potential energy .

Let us consider a generic system with $`N`$ particles with Hamiltonian $`H(C)`$, $`C`$ denoting the generic configuration of the system. The partition function both in the liquid phase at low temperature and in the glass phase can be written as

$`Z(\beta )={\displaystyle \underset{a}{}}\mathrm{exp}(\beta Nf_a(\beta ))=`$ (1)   
$`{\displaystyle 𝑑𝒩(f,\beta )\mathrm{exp}(\beta Nf)},`$ (2)   
where $`f_a(\beta )`$ is the free energy density of the valley labeled by $`a`$ at the temperature $`\beta ^1`$ and $`𝒩(f,\beta )`$ is the number of valleys with free energy density (per particle) less than $`f`$. In the glassy phase this sum is dominated by the valleys with minimal free energy, while in the liquid phase an exponentially large number of valleys contributes to the partition function. It is normally assumed that in glass forming systems, for large $`N`$ and $`f>f_0(\beta )`$, we have that

$$𝒩(f,\beta )=\mathrm{exp}(N\mathrm{\Sigma }(f,\beta )),$$
(3)   
where the configurational entropy, or complexity, $`\mathrm{\Sigma }(f,\beta )`$ is positive in this region and vanishes at $`f=f_0(\beta )`$. The quantity $`f_0(\beta )`$ is the minimum value of the free energy: $`𝒩(f,\beta )`$ is zero for $`f<f_0(\beta )`$.

The key point is the computation of the function $`\mathrm{\Sigma }(f,\beta )`$ for $`f>f_0(\beta )`$. It can be done using liquid theory because in this region we are still in the liquid phase. The location of the zero of $`\mathrm{\Sigma }(f,\beta )`$ will tell us the value of the free energy in the glassy phase.

At this end it is convenient to consider the generalized partition function

$$Z(\gamma ;\beta )\mathrm{exp}(N\gamma \mathrm{\Phi }(\gamma ;\beta ))=\underset{a}{}\mathrm{exp}(\gamma Nf_a(\beta )).$$
(4)   
The physical meaning of $`\gamma `$ is will be clearer later.

It is evident that

$$\gamma \mathrm{\Phi }(\gamma ;\beta )=\gamma f\mathrm{\Sigma }(\beta ,f),f=\frac{(\gamma (\mathrm{\Phi }(\gamma ;\beta )))}{\gamma }.$$
(5)   
The complexity can be simply obtained from $`\mathrm{\Phi }(\gamma ;\beta )`$ in the same way as the entropy can be obtained from the usual free energy.

The crucial step consists in writing

$`Z(\gamma ;\beta )=`$ (6)   
$`{\displaystyle 𝑑C\mathrm{exp}\left(\gamma H(C)N\gamma f(\beta ,C)+N\gamma f(\gamma ,C)\right)}=`$ (7)   
$`{\displaystyle 𝑑C\mathrm{exp}\left(\gamma H(C)N\gamma \widehat{f}(\beta ,C)+N\gamma \widehat{f}(\gamma ,C)\right)},`$ (8)   
where $`\widehat{f}(\beta ,C)=f(\beta ,C)f(\mathrm{},C)`$ and $`f(\beta ,C)`$ is a function that is constant in each valley and it is equal to the free energy density of the valley to which the configuration $`C`$ belongs. In deriving eq. (8) we have assumed that the value of $`\gamma `$ is high enough that all the configurations $`C`$, which contribute to the integral for large $`N`$, belong to some valley: We have also assumed that there is a one to one correspondence among the valley at inverse temperatures $`\beta `$ and at $`\gamma `$.

Before entering into the computation of $`\widehat{f}(\beta ,C)`$ it is useful to make the so called quenched approximation, i.e. to make the following approximation inside the previous integral:

$$\mathrm{exp}(A\widehat{f}(\beta ,C))=\mathrm{exp}(A\widehat{f}_\gamma (\beta )),$$
(9)   
where $`\widehat{f}_\gamma (\beta )`$ is the expectation value of $`\widehat{f}(\beta ,C)`$ taken with the probability distribution proportional to $`\mathrm{exp}(\gamma H(C))`$. The quenched approximation would be exact if the temperature dependance of the energy of all the valleys would be the same, apart from an overall shift at zero temperature. In other words we assume that the minima of the free energy have different values of the free energy but similar shapes. More refined computations can be done to compute systematically the corrections to the quenched approximation. This approximation would be certainly bad if we were using the free energy $`\widehat{f}(\beta ,C)`$ at the place of $`f(\beta ,C)`$ because the zero temperature energy strongly varies when we change the minimum.

We finally find

$$\mathrm{\Phi }(\gamma ;\beta )=F_L(\gamma )+\widehat{f}_\gamma (\beta )\widehat{f}_\gamma (\gamma ),$$
(10)   
$`F_L(\gamma )`$ being the free energy of the liquid ($`S_L(\gamma )`$ will be the entropy of the liquid). A simple algebra shows that

$$\mathrm{\Sigma }(\gamma ;\beta )=S_L(\gamma )S_\gamma (\gamma )+\widehat{f}_\gamma ^{}(\gamma )\widehat{f}_\gamma ^{}(\beta ),$$
(11)   
where $`\widehat{f}_\gamma ^{}(\beta )=\widehat{f}_\gamma (\beta )/\gamma `$.

In the liquid phase, we find out that the configurational entropy is given by

$$\mathrm{\Sigma }(\beta )=\mathrm{\Sigma }(\beta ;\beta )=S_L(\beta )S_\beta (\beta ),$$
(12)   
as expected: the entropy of the liquid is the entropy of the typical valley plus the configurational entropy.

The thermodynamic transition is characterized by the condition

$$\mathrm{\Sigma }(\beta _K)=0.$$
(13)   
In the glassy phase the free energy can be found by first computing the value of $`\gamma (\beta )`$ such that

$$\mathrm{\Sigma }(\gamma (\beta );\beta )=0$$
(14)   
and evaluating the corresponding free energy.

The quantity $`\gamma (\beta )`$ is the inverse of the effective temperature of the valley. It is easy to show (following ) that the previous formulae are completely equivalent to the replica approach.

A strong simplification happens if we assume that the entropy of the valley can be evaluated in the harmonic approximation where we only keep the vibrational contributions. In this case we obtain for a system with $`M`$ degrees of freedom that the harmonic entropy of the valley near to a configuration $`C`$ is given by

$$S(\beta (C))=\frac{M}{2}\mathrm{ln}\left(\frac{2\pi e}{\beta }\right)\frac{1}{2}\text{Tr}\left(\mathrm{ln}((C))\right),$$
(15)   
where $`(C)`$ is an $`M\times M`$ Hessian matrix (e.g. if $`H`$ depends on the coordinates $`x_i`$ we have that $`_{i,k}=^2H/x_ix_k`$).

If our approximation were fully consistent, we should find that all the eigenvalues of $``$ (the so called INN, Instantaneous Normal Modes ) are positive. This is not the case, however the number of negative eigenvalues becomes very small at low temperature, still in the liquid phase, signaling that valleys can be approximately defined in this region.

How can we evaluate the harmonic entropy in this case? Two rather similar possibilities are: (a) We compute $`\text{Tr}\left(|\mathrm{ln}((C))|\right)`$ instead of $`\text{Tr}\left(\mathrm{ln}((C))\right)`$. (b) We find the minimum of the Hamiltonian which is the nearest to $`C`$ (i.e. the corresponding inherent structure) and we use the spectrum of the Hessian at this point. We have checked numerically that the two methods give rather similar results  and we will follow the second one in the numerical simulations presented in the paper.

The framework has been set. We now segue into an explicit computation. Two possibilities are open: (a) We do an analytic computation of all the quantities which appear in the previous equations; (b) We extract them from numerical simulations. The first possibility is open only in relatively simple systems and further approximations are needed, the second one is viable for all systems.

In these letter we shall explore both possibilities in the case of a realistic model for glasses, the binary mixture of particles (80% large particles, 20 % smaller particles) interacting via a Lennard-Jones potential, introduced by Kob and Andersen . This Hamiltonian should mimick the behaviour of some metallic glasses and it is one of the best studied and simplest Hamiltonian which do not lead to crystalization at low temperature.

The numerical procedure does not present any serious difficulty. We have studied via Monte Carlo simulations a system of $`N=260`$ particles, in a cubic box with periodic boundary conditions at density $`\rho =1.2`$, starting from $`\beta ^{1/4}=0.02`$ and increasing $`\beta `$ by steps of $`\mathrm{\Delta }\beta ^{1/4}=0.02`$; we performed up to $`\mathrm{4\; 10}^6`$ MC steps at each $`\beta `$. We need to simulate the system at high temperatures in order to compute the entropy of the liquid (we use as reference entropy that of a perfect gas at $`\beta =0`$). The entropy is obtained using the formula $`S(\beta )=S(0)+_0^\beta 𝑑\beta ^{}(E(\beta )E(\beta ^{}))`$. Given an equilibrium configuration the near minimum of the potential is found by steepest descendent and the computation of the 780 eigenvalues of $`(C)`$ does not present any particular difficulty.

The results for the total entropy of the liquid and for the harmonic part are shown in fig. (1) as function of $`T^{.4}`$ (a more detailed description of the simulations can be found in ref. ). The entropy of the liquid is remarkable linear when plotted versus $`T^{.4}`$, as it happens in many cases .

In fig (2) we show the configurational entropy as function of $`T^{.4}`$. We fit it with a polynomial of second degree in $`T^{.4}`$. The extrapolated configurational entropy becomes zero at a temperature $`T_K=.31\pm .04`$, where the error contains systematic effects due to the extrapolation (similar conclusions have been reached in ref. ).

There are many methods to compute analytically the free energy in the liquid phase which lead to integral equations for the correlation functions. Here we follow a simple procedure proposed in ref. . It mixes the HNC (hypernetted chain) and MSA (mean spherical approximation) closures by means of a single parameter $`\alpha (T)`$ that is chosen in order to reduce thermodynamic inconsistencies, minimizing the difference between two different ways of computing the compressibility. This technique allow us to compute the internal energy in the liquid with a reasonable approximation. The resulting integral equations for the correlations function of the two kinds of particles are transformed in a set on non linear differential equations by discretizing them with a space resolution of $`2^5`$ and a large distance cutoff of $`2^4`$ (we have varied these parameters and checked that this choice gives a reasonable accuracy). The resulting equations in $`3512`$ unknown have been solved using a package from the IMSL library.

The computation of the spectrum is much more involved. We follow a simple approach which becomes exact if we assume that the range of the Hessian is much larger than the typical interatomic distance and we use the superposition approximation

$$g(x_1x_2x_3)=g(x_1x_2)g(x_2x_3)g(x_3x_1)$$
(16)   
in the cases where it is needed . In order to compute the eigenvalues of the Hessian we compute analytically the moments of the eigenvalue distribution

$$M_k=\text{Tr}(^k).$$
(17)   
Each moment can be written as the appropriate integral over the correlations functions. In this long range approximation some terms are dominant over the others. We keep only those terms and we introduce the superposition approximation into the appropriate places. We use the large range expansion to select the diagrams .

After some computations (which we do not report here ) we find a spectral density $`\rho (e)`$ that has support only in the region of positive eigenvalues $`e`$ (real frequencies $`\omega =\sqrt{e}`$) and goes to zero as $`\omega `$ as expected.

We can now put everything together: the final analytic predictions for the liquid and harmonic entropy are shown in fig. (1). The liquid entropy turns out to be very good, while there is a minor discrepancy for the harmonic entropy, likely due to the rather strong approximations we have done in the analytic computation. The analytic configurational entropy is shown in fig. (2). It becomes zero at $`T_K=.32`$, which is our analytic prediction for the thermodynamic transition (as a check we have fitted the analytic configurational entropy using the same procedure as for the numerical one and we have found $`T_K=.34`$).

We have in our hands all the tools to compute analytically the free energy in the low temperature case. In fig. (3) we show the specific heat coming from the $`x`$ dependent part of the Hamiltonian as a function of the temperature (we must add the momentum contribution $`3/2`$ in order to get the total specific heat). Very similar results are obtained by using the extrapolated numerical entropies at the place of the analytic ones. We notice that the Dulong Petit law is extremely well satisfied in the low temperature region. Indeed in the harmonic approximation the Dulong Petit law would be exact if we neglect the $`\gamma `$ dependence of $`S_\gamma (\beta )`$. The value of $`\gamma (\beta )`$ weakly depends on $`\beta `$: its value in the limit $`\beta \mathrm{}`$ is only about 10% higher that its value (i.e. $`\beta _K`$) at the transition temperature.

The results presented here are an explicit realization of a transition driven by an entropy crisis, which has been firstly implemented microscopically in the random energy model . The specific heat has a jump downward, when we decrease the temperature, which is the opposite of the typical behaviour in transitions characterized by the onset of a conventional order parameter (in mean field theory ferromagnets, superconductors, have a jump upward when we decrease the temperature).

We stress that it is possible to remove the approximation of using the harmonic entropy. There are no serious difficulties in computing numerically the true entropy of the valleys; this has been done for a binary mixture of soft sphere and the results are very near to the harmonic ones . Analytically we can also expand in the anharmonicity parameters. It is also possible to take care of the fluctuations of the entropy from valley to valley and go beyond the quenched approximation. It is quite likely that these effects will not strongly change the results. The most important step would be to reach a better theoretical control on the spectrum of the INN in the liquid phase, maybe combining the approximations used here with the low density expansion of .

Summarizing we have found a method which is able to use liquid theory method in the glasses phase putting in practice the old adage a glass is a frozen liquid. We are able to compute with a reasonable approximation the thermodynamics and with a little more effort we can compute the static and the dynamic structure functions.

We owe a lot to M.Mézard and we are very happy to thanks him. We are also very grateful to W.Kob and F.Sciortino for their suggestions and to A.Cavagna and I.Giardina for useful discussions. B.C. would like to thank the Physics Department of Rome University ‘La Sapienza‘ where this work was partially developed during her PhD.