# Energy Landscape and Overlap Distribution of Binary Lennard-Jones Glasses

## 1 Distance Measure in Configuration Space

We study correlations among an ensemble of such metastable states, labelled by $`a`$. A glassy state $`a`$ is characterized by a spontaneous breaking of translational invariance as signalled by a non-zero expectation value of the Fourier-components $`\mathrm{exp}(i\stackrel{}{k}\stackrel{}{r}_i)_a0`$ of the local density. Here $`\stackrel{}{r}_i`$ denotes the position of the $`i`$-th particle, $`i=1,\mathrm{},N`$. The local density is analogous to the local magnetization, which is non-zero in the spin glass phase, $`s_i_a0`$, indicating the spontaneous breaking of the Ising symmetry. An overlap between two glassy configurations can be defined as

$$Q^{a,b}(\stackrel{}{k})=\frac{1}{N}\underset{i=1}{\overset{N}{}}\mathrm{exp}(i\stackrel{}{k}\stackrel{}{r}_i)_a\mathrm{exp}(i\stackrel{}{k}\stackrel{}{r}_i)_b$$
(1)   
in naive analogy to the spin glass order parameter $`Q^{a,b}=\frac{1}{N}_{i=1}^Ns_i_as_i_b`$. For the glassy state, however, $`Q^{a,b}(\stackrel{}{k})`$ is not a good measure, because it depends on the labelling of the particles. In particular two identical configurations with different labels can have a very small overlap according to the above definition. Permutations of identical particles do not give rise to new states in the space of physically distinguishable configurations. Hence, we need a distance measure which is invariant under all permutations $`\pi `$ of identical particles in one state only. This can be achieved by maximizing over all such permutations

$$Q^{a,b}(\stackrel{}{k})=\underset{\pi }{\mathrm{max}}\frac{1}{N}\underset{i=1}{\overset{N}{}}\mathrm{exp}(i\stackrel{}{k}\stackrel{}{r}_i)_a\mathrm{exp}(i\stackrel{}{k}\stackrel{}{r}_{\pi (i)})_b.$$
(2)   
This overlap has the required invariance property, but its computation requires the enumeration of all permutations, which is computationally hard and can at most be done for very small systems.

The computationally hard problem can be avoided by working with densities. We specialize to zero temperature and only consider configurations of locally minimal energy. Such a configuration is uniquely characterized by a point measure $`\rho ^a(\stackrel{}{r})=(1/N)_{i=1}^N\delta (\stackrel{}{r}\stackrel{}{r}_i^a)`$, which by construction is independent of the labelling of the particles. To obtain a distance measure in configuration space, we regularise the point particle densities $`\delta (\stackrel{}{r})`$ to homogeneous $`\eta `$-spheres

$$\delta (\stackrel{}{r})\mathrm{\Delta }_\eta (\stackrel{}{r})=\frac{1}{V_\eta }\theta (\eta |\stackrel{}{r}|),$$
(3)   
Here $`V_\eta =4\pi \eta ^3/3`$ is the volume of the $`\eta `$-sphere and we choose $`2\eta <r_{\mathrm{min}}^a=\mathrm{min}_{i,j}(|\stackrel{}{r}_i^a\stackrel{}{r}_j^a|)`$ to guarantee that the spheres do not overlap and that the positions $`\stackrel{}{r}_i^a`$ can still be uniquely reconstructed from the regularised density. Such a regularisation is necessary because products of $`\delta `$ functions are ill defined and the probability of coincidence of any two 3-dimensional vectors which are distributed smoothly in space is equal to zero.

A natural distance measure between two regularised densities is

$$Q^{a,b}=\frac{V_\eta }{N}\mathrm{d}^3x\underset{i,j=1}{\overset{N}{}}\mathrm{\Delta }_\eta (\stackrel{}{x}\stackrel{}{x}_i^a)\mathrm{\Delta }_\eta (\stackrel{}{x}\stackrel{}{x}_j^b).$$
(4)   
If we restrict $`\eta `$ to even smaller values, $`4\eta <r_{\mathrm{min}}:=\mathrm{min}(r_{\mathrm{min}}^a,r_{\mathrm{min}}^b)`$, then each sphere in state $`a`$ can overlap with at most one sphere in state $`b`$. Consequently the product of two local regularised densities is only non-zero, for the permutation $`j=\pi (i)`$, which identifies particles within the same $`\eta `$-sphere. The overlap is then given by the volume fraction of overlapping spheres in the two configurations, so that $`0Q^{a,b}1`$.

It is obvious that $`\eta `$ cannot be made arbitrarily small, because then the distance measure looses its ability to discriminate between different structures. Very small $`\eta `$ will simply lead to maximal distances in most cases and ultimately reproduce the point measure. We have varied $`\eta `$ over the range $`0<\eta <r_{\mathrm{min}}/4`$ and computed the number of overlapping spheres for two configurations which are random and two configurations which are similar (about 80% overlap). We find that the number of overlapping spheres is hardly sensitive to the choice of $`\eta `$ in the range $`r_{\mathrm{min}}/8<\eta <r_{\mathrm{min}}/4`$. For $`\eta <r_{\mathrm{min}}/8`$ the number of overlapping spheres decreases drastically with $`\eta `$ even for the strongly correlated configurations so that $`Q^{a,b}`$ can no longer discriminate between strongly and poorly correlated states. In the following we fix $`\eta =r_{\mathrm{min}}/5`$.

Other distance measures have been used in the literature , in particular

$$D(a,b)=\underset{\pi }{\mathrm{min}}\underset{i=1}{\overset{N}{}}\left(\stackrel{}{r}_i^a\stackrel{}{r}_{\pi (i)}^b\right)^2,$$
(5)   
which requires the solution of a hard computational problem. It has the further disadvantage to be not directly related to quantities which may be obtained from experiment (like the density) or which appear naturally in any of the existing theories of supercooled liquids and structural glasses.

A system in free space with pairwise interactions possesses translational and rotational symmetries, which may be broken explicitly by boundary conditions. The most convenient boundary conditions for computer simulations, namely a cubic box with periodic boundary conditions, leave a subgroup of the free space symmetries unbroken. This subgroup $`𝒮`$ is generated from the translations and the cubic point group symmetries. As a consequence, there exists a whole orbit of symmetry related, degenerate minima of $`H`$ to every single minimum found in the simulations.

We are interested in the rate of occurrence, i.e. the probability distribution, of overlaps $`P(Q)=_{a,b}P_aP_b\delta (QQ^{a,b})`$. Here the summation includes all symmetry related states, which all have the same rate of occurrence $`P_a`$. The large number of symmetry related states imposes two severe difficulties. First, many translated states (out of a continuum of states) have to be generated to obtain a reliable estimate of $`P(Q)`$ from numerical simulations. Second, the $`Q`$ values generated by all symmetry related states may scatter considerably in the interval $`0Q1`$. This may smear out structures in $`P(Q)`$ which would be clearly marked in an ensemble without residual symmetries. To avoid this problem we have broken the translational symmetry explicitly by fixing the centre-of-mass. Since total momentum is conserved in our simulations, this can be achieved by an appropriate choice of initial conditions (see below). The cubic point group symmetry has been left untouched: for each minimum which was found by our algorithm, we have generated 48 equivalent states and included them in the histogram for $`Q^{a,b}`$. To check, whether this procedure smears out relevant structure, we have also maximized the overlap over all 48 equivalent states.

## 2 Simulations

The system under consideration is a binary mixture of large (L) and small (S) particles with $`80\%`$ large and $`20\%`$ small particles. Small and large particles only differ in diameter, but have the same mass. They interact via a Lennard–Jones potential of the form $`U_{\alpha \beta }(r)=4ϵ_{\alpha \beta }[(\sigma _{\alpha \beta }/r)^{12}(\sigma _{\alpha \beta }/r)^6]`$. All results are given in reduced units, where $`\sigma _{LL}`$ was used as the length unit and $`ϵ_{LL}`$ as the energy unit. The other values of $`ϵ`$ and $`\sigma `$ were chosen as follows: $`ϵ_{LS}=1.5,\sigma _{LS}=0.8,ϵ_{SS}=0.5,\sigma _{SS}=0.88`$. The systems were kept at a fixed density $`\rho 1.2`$. Periodic boundary conditions have been applied and the potential has been truncated appropriately according to the minimum image rule , and shifted to zero at the respective cutoff. For the systems with $`N=30`$ particles the cut-off is $`r_c=1.43`$ and for $`N=60`$ the cut-off is at $`r_c=1.7`$. The choice of the Lennard-Jones parameters follows recent simulations of Lennard-Jones glasses ; it is known to suppress recrystallization of the system on molecular dynamics time scales. The glass transition is supposed to occur at the temperature $`T_g0.45`$ . Throughout this study we present results for systems with $`N=30`$ and $`N=60`$ particles, noticing that most of the results have been verified for $`N=100`$.

Initially $`N`$ atoms are placed randomly inside a cubic simulation box of side-length $`L`$. The following steps are performed repeatedly:

* Heat up the system to a temperature $`T_{\mathrm{run}}`$.
* Let the system evolve for a time $`\tau _{\mathrm{run}}`$ using molecular dynamics.
* Locate the nearest local minimum by quenching down the system to $`T=0K`$ using the steepest descent path.

Explicit breaking of translational invariance can be achieved by fixing the centre-of-mass of the system. A unique definition of the centre-of-mass $`\stackrel{}{R}_s`$ is not obvious for periodic boundary conditions. We start from the extended zone scheme representation of the simulation box. The infinite periodic density pattern is cut off at finite volume $`V`$, so that the cubic point group symmetry stays intact. The centre-of-mass is the same for all such volumes, including the smallest one, the simulation box itself. The regularised $`\stackrel{}{R}_s`$ remains unchanged in the limit $`V\mathrm{}`$.

Each run is initialized with the same positions and different velocities of the particles. The latter are drawn from a Maxwell-Boltzmann distribution with temperature $`T_{\mathrm{run}}`$, subject to the constraint that the total momentum vanishes. The total momentum is conserved in the molecular dynamics simulation as well as in the steepest descent algorithm. Thereby we generate an ensemble, such that in all states the centre-of-mass stays fixed. In the following we shall also consider states which are related by cubic point group symmetries. To guarantee that the centre of mass remains unaffected by point group operations we choose $`\stackrel{}{R}_s=0`$. We use standard molecular dynamics with the velocity form of the Verlet algorithm . Steepest descent is achieved with the conjugate gradient algorithm .

## 3 Results

Simulating small systems requires careful checks to avoid sampling crystalline configurations. First, it can be quite helpful to look at the structures. Second, we have applied a common neighbor analysis  to detect the amount of polytetrahedral order (amorphous) and closed-packed structures (crystalline). Occasionally the structural stability of low energy states has been verified for temperatures below the glass transition.

In fig. 2 we show a histogram of the energy per particle of locally minimal states for $`N=30`$ and $`N=60`$. For each system size $`1000`$ minima have been generated, starting from $`T_{\mathrm{run}}=0.5`$. For the small system, crystalline states have the dominating basins of attraction, whereas for the larger system we found very few crystalline states (3 out of 1000).

Energy Landscape. The set of sampled local minima depends on $`T_{\mathrm{run}}`$. We have generated $`200`$ configurations for $`N=60`$ particles and fixed $`T_{\mathrm{run}}`$ chosen in the interval $`0.3T_{\mathrm{run}}T=1.8`$. For each $`T_{\mathrm{run}}`$ we have computed the mean energy, the lowest energy found and the variance of the energy. The results are shown in fig. 2. For each configuration the length of the molecular dynamics simulation has been chosen as $`\tau _{\mathrm{run}}=0.510^6`$. This is sufficient to allow for equilibration at $`T_{\mathrm{run}}>T_g`$.

The mean energies per particle of the sampled minima are approximately constant ($`ee_h=5.2`$) at high temperatures and decrease significantly within a transition region. The lowest energy minima $`e_l5.45`$ are found for $`T_{\mathrm{run}}T_g`$, whereas for $`T_{\mathrm{run}}<T_g`$ sampling is restricted to one quasi-ergodic component. Thus extensive sampling of locally stable states is most effective for $`T_{\mathrm{run}}T_g`$. At high temperatures, most of the weight is found in a broad band of locally stable states with high energies around $`e_h`$ . The low energy structures are not sampled, although these minima are favoured by a relative Boltzmann factor $`\mathrm{exp}(N(e_he_l)/k_BT_{\mathrm{run}})`$. This implies that the configuration space of the high energy basins is much larger than that of the low energy basins, so that it overcompensates the relative Boltzmann factor. When $`T_{\mathrm{run}}`$ is lowered to temperatures around $`T_g`$, the Boltzmann factor becomes more important and eventually overruns the configuration space volumes of the $`e_h`$ states. These results are in agreement with recent work of Sastry et al. , who showed that the transition is accompanied by the onset of typical, glass like relaxation phenomena.

To assess the performance of our algorithm in sampling low lying energy states quantitatively, we compare the lowest energies sampled with the results of a recent study on the optimization properties of traditional simulation algorithms applied to the Lennard-Jones glasses . The authors reported on a lowest energy configuration found for $`N=60`$ at $`e=5.38`$ by cooling down slowly from the liquid phase using standard molecular dynamics. The lowest state shown in fig. 2 is at $`e=5.45`$. The difference in energy $`\delta e=0.07`$ corresponds to a factor of $`100`$ in computer time according to the estimate given in refs. . For $`T_{\mathrm{run}}0.5`$ around $`15\%`$ of all states sampled using our method are lower in energy then $`e=5.38`$.

Distribution of overlaps. To investigate the distribution of overlaps we have generated $`2000`$ amorphous states for the $`N=60`$ particles system keeping the centre-of-mass fixed. Each molecular dynamics trajectory has been initialized at $`T_{\mathrm{run}}=0.5`$ and has been run for $`\tau _{\mathrm{run}}=10^6`$ time steps. We compute the overlap of the large particles for all pairs $`(a,b)`$. This implies that the weights $`P_a`$ are approximated by the rate of occurrence of state $`a`$ in our simulation. Several arguments support our choice of ensemble: 1) The generated states are as low in energy as the best ones from other optimization procedures. 2) In agreement with recent simulations of Sastry et al. , we find a qualitative change in the properties of sampled states at a well defined temperature $`T_g0.5`$, such that above $`T_g`$ a broad range of states are sampled, whereas below $`T_g`$ the system is landscape dominated, i.e. trapped in a metastable state. 3) The Stillinger-Weber method is controlled and reproducible and sufficiently simple to allow for analytical calculations.

A histogram of overlaps is shown in fig. 4. We observe a most probable value around $`Q^{a,b}0.2`$ and no significant structure in the distribution. It is instructive to also compute the distribution of overlaps of liquid configurations. This distribution is compared to the histogram for glassy states in fig. 4. The distribution of overlaps for glassy states is $`35\%`$ broader and centered at a $`20\%`$ higher value. Otherwise no significant changes occur. Thus we conclude that the glassy minima are as different as they can be. This conclusion also holds, if we restrict our ensemble of glassy states to the lowest energy states.

The residual cubic point group symmetry does not affect the overlap distribution significantly. We have calculated overlaps by either maximizing over the $`48`$ symmetry related states or by including all symmetry related states with the same weight. The resulting distributions are shown in Fig. 4. As one would expect the distribution obtained by maximizing over all rotations is peaked at a higher value of $`Q`$, but does not reveal any additional structure.

Recently Coluzzi and Parisi  computed the distribution of distances $`P(D)`$, using the distance measure of eq. (5) and simulating small systems ($`28`$ to $`36`$ particles) confined by soft walls. They find highly non-trivial distributions $`P(D)`$ in the glassy phase, which strongly depend on particle number. We have also applied our method to generate an ensemble of low energy states by rapid cooling for the system of ref. . A straightforward evaluation of $`P(Q)`$ yields the bimodal distribution, shown in fig. 5. A closer inspection, however, reveals that the structure is due to many imperfect crystalline states. Soft walls increase the probability for crystalline states, as compared to periodic boundary conditions, and for small systems ($`N30`$ particles) most of the particles sit on the surface of the sample.