The Glass Transition: A Topological Perspective

1. Introduction

As of today, the glass transition still stands as an open problem of contemporary physics. Indeed, as glasses are amorphous solids, no symmetry breaking is associated with this transition, which is therefore out of the domain of applicability of Landau theory [1]. It has further been argued that glass forming materials (at least strong glass formers in Angell’s classification [2]) do not even exhibit a transition in the conventional sense, as it yields no dynamical singularity; the transition temperature $T_g$ is usually regarded as purely conventional, and defined by the passing of a threshold value of the viscosity or of the relaxation time.

In a number of references [3,4,5], it has been shown that glass transitions most likely correspond to geometric transition. Namely, the critical points of the potential energy were studied, and in particular their instability index (i.e. the number of negative eigenvalues of the Hessian matrix). It was found that the average index density vanishes at the so-called mode-coupling temperature (MCT).

Yet it is known, from Morse theory [1,6], that any change of stability indices of a surface is accompanied by a change of its topology.

Furthermore, the relatively recent topological theory of phase transitions [1,7,8,9], unraveled a deep link between classical phase transitions and changes in the topology of the potential level sets (PLS) ${\Sigma }_{\Phi }$, i.e. the iso-energy hypersurfaces. One advantage of the topological theory of phase transitions is that it applies to small systems (mesoscopic and nanoscopic scales), thus escaping the thermodynamic limit dogma upon which is built the Yang–Lee theory of phase transitions. Furthermore, it applies to phase transitions in the absence of symmetry-breaking (hence in the absence of a well-defined order parameter). That second point is of great interest to us, as the glass transition notoriously falls in the latter category.

These observations compel us to investigate glass-forming systems, resorting to a few elegant theorems linking the topological invariants and geometric quantities such as the mean curvature.

The topological theory of phase transition is usually defined at equilibrium, as it requires the study of the whole PLS, regardless of metastability and dynamical slowing down. This could be considered an important drawback to its application to the glass transition, considered an inherently out-of-equilibrium phenomenon. However, glass formation can be regarded as a transitional, metastable state, induced by the slowing down of the dynamics due to a frustration on the way to crystallization [10,11]. This phenomenon is in turn necessarily rooted into the equilibrium properties of matter. We thus deem in principle relevant such a study of a glass-former’s equilibrium properties.

Moreover, in systems that do not exhibit a proper crystalline phase due to a high degree of frustration, the free energy global minimum (or, in the microcanonical ensemble, the entropy global maximum at a given energy) might very well be in fact a disordered, glassy state. It seems that, in this situation, the glass phase is not considered an equilibrium state for the sole reason that it breaks ergodicity, the system being stuck into potential wells for extensive amounts of times.

While this proposition holds by definition of thermodynamic equilibrium, we argue for its reevaluation under some specific conditions. In this picture of a highly frustrated system, the PLS at low energy is partitioned into a multitude of isolated wells, from which we cannot escape, in realistic conditions. However, ergodicity can be artificially restored, resorting to the variety of numerical tricks that allow the system to “jump out” from well to well. In such a situation, the glass phase can be studied as if it was an equilibrium state, in the sense that measurements can be performed out of a state-space importance sampling; provided that there is no possible crystalline configuration that can “corrupt” this sampling, and that only glass states can be found at the energy of interest, we can then say that we somewhat performed equilibrium measurements of the glass phase.

In Section 3.1, we present the model of glass-former under investigation. Section 2.1 is devoted to the presentation of geometric quantities that are linked to topology by a handful of useful theorems; we show how these quantities can be measured in our simulations. In Section 2.2, we present our numerical results. We found a two-step transition marked by peaks of the specific heat, that we classify as a second order and a weakly first order transition points, using analytic tools developed specifically for the microcanonical ensemble. Furthermore, we show that these transitions correspond to jumps of bond-orientational order parameters, as well as modifications of the spatial distribution profiles; these critical points are thus accompanied by modifications of the short-range structural properties of the system. Finally, we observe singular behaviors of various geometric quantities, in correspondence with the observed transition, conclusively implying an underlying change in the topology of the potential energy level sets.

2. Material and Methods

2.1. Geometric Signatures of Topological Changes

Entropy stems as a fundamental building block of thermodynamics. In particular, in the microcanonical ensemble, most macroscopic observables can be retrieved from its derivatives. Furthermore, it has recently been proposed [12] a microcanonical ensemble classification of phase transitions, analogous to the notorious classification of Ehrenfest, heuristically associating first and second order phase transitions to discontinuity of the second and third derivatives, respectively, of the entropy.

In standard Hamiltonian systems as ours, that is where H is a quadratic function of the momenta, the kinetic part of the canonical partition is known to be trivial, as it reduces to a constant factor. In the microcanonical ensemble, the dissociation of the kinetic and configurational parts of the partition function is somewhat less evident, but can nevertheless be performed through Laplace transform techniques [13]; in particular, this separation allows for the practical expression of the microcanonial probability density (A2). It results that the relevant information is entirely contained in the configurational entropy

$$S\left(\phi \right)={\displaystyle \frac{1}{3N}}log\int d\Gamma \Theta \left(\phi -\Phi (\Gamma )\right),$$

(1)

where $\Theta \left(x\right)$ is the Heaviside step function, which vanishes for $x<0$ and equates one for $x\ge 0$. Yet the latter expression can be rewritten in terms of the PLS volumes [1,14]

$$S\left(\phi \right)={\displaystyle \frac{1}{3N}}log{\int }_0^{\phi }d\varphi {\int }_{{\Sigma }_{\varphi }}{\displaystyle \frac{d\sigma }{|\nabla \Phi |}},$$

(2)

∇ being here the gradient operator, $d\sigma$ is the elementary volume induced by the immersion in ${\mathbb{R}}^{3N}$ of the PLS ${\Sigma }_{\varphi }$, hypersurface of dimension $3N-1$ defined as

$${\Sigma }_{\varphi }=\left\{\Gamma \in \Omega |\Phi \left(\Gamma \right)=\varphi \right\},$$

(3)

where $\Omega$ is the full state space, and $\varphi$ is the fixed value defining the PLS. Finally, it has been shown that Eq. (Section 2.1) can be expressed in terms of topological invariants, namely [1]

$$S\left(\phi \right)={\displaystyle \frac{1}{3N}}log\left[Vol\left({\mathbb{S}}_1^{3N-1}\right)\sum _{i=0}^{3N}{b}_i\left({\Sigma }_{\phi }\right)+{\mathcal{R}}_1\left(\phi \right)\right]+{\displaystyle \frac{1}{N}}log{\mathcal{R}}_2\left(\phi \right),$$

(4)

where ${\mathcal{R}}_1,{\mathcal{R}}_2$ are smooth functions of the potential, $Vol\left({\mathbb{S}}_1^{3N-1}\right)$ is the volume of the unit ball of dimension $3N-1$, and $b_i\left({\Sigma }_{\phi }\right)$ is the ith Betti number of the manifold ${\Sigma }_{\phi }$.

Eq. (Section 2.1) highlights the dependence of the configurational entropy on topological aspects of the PLS, encoded in the Betti numbers $b_i$. This observation was at the root of the topological hypothesis, stating that the deep mathematical origin of a phase transition was to be found in a topological change of the PLS. It is worth noting that, whereas any phase transition is rooted in a topological change, not all topological changes entail a phase transition.

In the present work, we thus aim at establishing a correspondence of the glass transition with topological changes of the PLS.

Probing the topology of the high-dimensional manifolds that are the PLS is by no means a simple task. To our best knowledge, there exists no way of fully characterizing it by means of measurable average observables, namely the tools accessible to us. For lack of a complete reconstruction of the topology of the submanifolds of interest, it is however possible to probe topological changes, which are, in the end, our true object of study. There fortunately exist a few theorems of differential topology drawing sufficiently strong links between geometrical and topological quantities, allowing us to observe, when they are present, sharp topological changes.

We now introduce very roughly a few notions of differential extrinsic geometry, that will be useful to the development of our topological probing. For a more extensive development of this framework, we refer the reader to Ref. [1,14,15].

In order to alleviate our notations, we now drop the dependence in $\Gamma$.

We first present Pinkall’s theorem, relating the average dispersion of principal curvatures with the weighted sum of the Betti numbers

$${\displaystyle \frac{{\int }_{{\Sigma }_{\varphi }}{\sigma }_{\kappa }^2(\Gamma )d\Gamma }{{\int }_{{\Sigma }_{\varphi }}d\Gamma }}={\left[Vol\left({\mathbb{S}}_1^D\right)\sum _{i=1}^D{\left({\displaystyle \frac{i}{D-1}}\right)}^{D/2-i}{b}_i\left({\Sigma }_{\varphi }\right)\right]}^{2/D}-r\left({\Sigma }_{\varphi }\right),$$

(5)

where D is the dimension of the manifold, ${\sigma }_{\kappa }^2=\langle {\kappa }_i^2\rangle -{\langle {\kappa }_i\rangle }^2$ is the dispersion of the principal curvatures, and $r\left({\Sigma }_{\varphi }\right)$ is a remainder, which stays small provided that ${\sigma }_{\kappa }^2$ doesn’t exhibit too large variations on the submanifold ${\Sigma }_{\varphi }$.

Another geometric quantity that connects to topological invariants is the range of variability of the sectional curvatures ${\Delta }_{sec}$. Overholt’s theorem indeed states that it provides an upper bound to the sum of Betti numbers

$${\Delta }_{sec}\ge {\left[{\displaystyle \frac{Vol\left({\mathbb{S}}^D\right){\sum }_{i=0}^D{b}_i\left({\Sigma }_{\varphi }\right)}{2Vol\left({\Sigma }_{\varphi }\right)}}\right]}^{2/D}.$$

(6)

In turn, ${\Delta }_{sec}$ is related to the variance of the scalar curvature $R_{\Sigma }$, as the latter is simply defined as the sum of all the sectional curvatures at a given point $\Gamma$

$$R_{\Sigma }=\sum _{i\ne j}{K}_{ij}=\sum _{i\ne j}{\kappa }_i{\kappa }_j,$$

(7)

where $K_{ij}$ is the sectional curvatureof sectional plane $(u_i,u_j)$, and ${\left\{u_j\right\}}_{j=1,\cdots ,D}$ forms an orthonormal basis in the tangent space at this point.

It results

$${\displaystyle \frac{\langle R_{\Sigma }^2\rangle -{\langle R_{\Sigma }\rangle }^2}{N(N-1)}}\approx {\Delta }_{sec},$$

(8)

Up to our best knowledge, the simplest way to compute these quantities in the context of numerical simulations, is by considering the Weingarten operator, also called shape operator. A most useful tool characterizing the extrinsic geometry of hypersurfaces, it is the operator such that, for $\mathit{X}\in T\Sigma$ a vector field in the tangent bundle of $\Sigma$, we have

$${\mathcal{W}}_{\mathit{n}}\left(\mathit{X}\right)=-{\nabla }_{\mathit{X}}\mathit{n},$$

(9)

where ${\nabla }_{\mathit{X}}$ is the Levi-Civita connexion on $\Sigma$, and

$$\mathit{n}={\displaystyle \frac{\nabla \varphi }{|\nabla \varphi |}}\left(\Gamma \right)$$

is the normal to $\Sigma$ at a given point $\Gamma$.

The trace of the shape operator and of its square can be expressed in terms of mere derivatives of $\Phi$, namely

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill Tr\left[{\mathcal{W}}_{\mathit{n}}\right]& ={\displaystyle \frac{\Delta \Phi }{|\nabla \Phi |}}-{\displaystyle \frac{\nabla {\Phi }^{*}·Hess\left(\Phi \right)·\nabla \Phi }{{|\nabla \Phi |}^3}}\hfill \\ \hfill Tr\left[{\mathcal{W}}_{\mathit{n}}^2\right]& ={\displaystyle \frac{Tr\left[Hess{\left(\Phi \right)}^2\right]}{{|\nabla \Phi |}^2}}+{\displaystyle \frac{{|\nabla {\Phi }^{*}·Hess\left(\Phi \right)·\nabla \Phi |}^2}{{|\nabla \Phi |}^6}}-2{\displaystyle \frac{{|Hess\left(\Phi \right)·\nabla \Phi |}^2}{{|\nabla \Phi |}^4}},\hfill \end{array}\end{array}$$

(10)

where $\Delta \Phi$ and $Hess\left(\Phi \right)$ denotes respectively the Laplacian and the Hessian of $\varphi$, and “·” the scalar product.

The eigenvalues of ${\mathcal{W}}_{\mathit{n}}$ are the D principal curvatures ${\kappa }_i$. It results that the above-mentioned geometric quantities can all be expressed with combinations of $Tr\left[{\mathcal{W}}_{\mathit{n}}\right]$ and $Tr\left[{\mathcal{W}}_{\mathit{n}}^2\right]$, namely

$$\begin{array}{cc}\hfill M_{\Sigma }& ={\displaystyle \frac{Tr\left[{\mathcal{W}}_{\mathit{n}}\right]}{D}}\hfill \\ \hfill {\sigma }_{\kappa }^2& ={\displaystyle \frac{Tr\left[{\mathcal{W}}_{\mathit{n}}^2\right]}{D}}-{\displaystyle \frac{Tr{\left[{\mathcal{W}}_{\mathit{n}}\right]}^2}{D^2}}\hfill \\ \hfill R_{\Sigma }& =Tr{\left[{\mathcal{W}}_{\mathit{n}}\right]}^2-Tr\left[{\mathcal{W}}_{\mathit{n}}^2\right].\hfill \end{array}$$

(11)

The combination of formulae (10) and (11) clearly provide a straightforward way of obtaining the quantities of interest in the context of numerical simulations, simply by computing and combining the gradient and Hessian of the potential function $\Phi$ at each measurement step.

The latter geometric quantities pertain to the geometrical characteristics of the PLS, while our simulations are performed at constant total energy E. However, at large N, the fluctuations of $\Phi$ and K tend to vanish, and surfaces of constant $K\left({\left\{{\mathbf{p}}_i\right\}}_{i=1,\cdots ,N}\right)={\sum }_i{p}_i^2/2$ are diffeomorphic to $3N$-hyperspheres; the energy level sets can then be seen as product manifolds ${\Sigma }_E\sim {\mathbb{S}}_K^{3N}\times {\Sigma }_{\varphi }$. We thus consider $\varphi \left(E\right)$ stable enough for the corresponding PLSs to be diffeomorphic to one another, and for the general behavior of the above-defined geometric quantities to be trusted.

2.2. Numerical Methods

Using a microcanonical ensemble Monte Carlo scheme described in the Appendix A, we explored the behavior of the model defined by (12).

We simulated a system of size $N=216$ particles. Periodic boundary conditions and smooth cutoffs have been implemented, described in Appendix B.

It is worth noting that this simulations was, as is often the case in so-called glassy systems, very time consuming and hard to equilibrate.

An exact estimation of the total computation time is in practice difficult to assert, partly due to the fact that the set of energies we considered was changed multiple times during this extensive work. To provide a rough idea of the involved time scales, the equilibration phase, added with the earlier simulations aimed at code-testing and optimization, took more than 400 hours in CPU time per replica, with 50 replicas of the system. The efficient computation time, over which we performed retrieved the equilibrium data presented in this work, was 175 hours per replica, with 100 replicas.

Throughout the following of this article, we present data acquired at equilibrium, over $5.7·{10}^6$Monte Carlo sweeps (that is, a Monte Carlo step per particle), with a sampling rate of $1/1000$. Replica exchanges were attempted every 2000 sweeps.

3. Results

3.1. Model

The system we chose here to consider consists in a binary Lennard-Jones mixture, first introduced in [16], of Hamiltonian $H\left({\left\{{\mathit{q}}_i\right\}}_{i=1,\cdots ,N}\right)=K+\Phi \left({\left\{{\mathit{q}}_i\right\}}_{i=1,\cdots ,N}\right)$, where K is the total kinetic energy, the ${\mathit{q}}_i$ are the N 3-dimensional position variables, and

$$\begin{array}{cccccccccccccc}\hfill 37& \Phi (\Gamma )=\hfill & \hfill & {\Phi }_{11}(\Gamma )\hfill & \hfill & +\hfill & \hfill & {\Phi }_{22}(\Gamma )\hfill & \hfill & +\hfill & & {\Phi }_{12}(\Gamma )& & \\ \hfill & =\hfill & \hfill & \sum _{i,j\in {\Lambda }_1}4{ϵ}_{11}{\left({\displaystyle \frac{{\sigma }_{11}}{r_{ij}}}\right)}^{12}\hfill & \hfill & +\hfill & \hfill & \sum _{i,j\in {\Lambda }_2}4{ϵ}_{22}{\left({\displaystyle \frac{{\sigma }_{22}}{r_{ij}}}\right)}^{12}\hfill & \hfill & +\hfill & & \sum _{\begin{array}{c}i\in {\Lambda }_1,\\ j\in {\Lambda }_2\end{array}}4{ϵ}_{12}\left[{\left({\displaystyle \frac{{\sigma }_{12}}{r_{ij}}}\right)}^{12}-{\left({\displaystyle \frac{{\sigma }_{12}}{r_{ij}}}\right)}^6\right]& & ,\end{array}$$

(12)

where we introduced the shorthand notation $\Gamma ={\left\{{\mathit{q}}_i\right\}}_{i=1,\cdots ,N}$ for the instantaneous configuration of the system, ${\Lambda }_1,{\Lambda }_2$ are the set of particles belonging to species 1 and 2 respectively, and $r_{ij}=|{\mathit{q}}_i-{\mathit{q}}_j|$. The interaction parameters are set as

$$\begin{array}{cccc}\hfill {\sigma }_{22}/{\sigma }_{11}& =0.85\hfill & \hfill {\sigma }_{12}/{\sigma }_{11}& =0.49\hfill \\ \hfill {ϵ}_{12}/{ϵ}_{11}& =6\hfill & \hfill {ϵ}_{22}/{ϵ}_{11}& =1,\hfill \end{array}$$

and ${ϵ}_{11}=1$, ${\sigma }_{11}=1$. The density $\rho =1.6$ and the respective concentrations of the two species are $c_1\approx 0.33$ and $c_2\approx 0.67$. Precisely, we considered in our numerical study a sample of $N_1=75$ particles of species 1, and $N_2=145$ particles of species 2, for a total of $N=216$ particles. As is it customary, we set the Boltzmann constant $k_B=1$.

In this numerical study, we do not take into account the microscopic details of the kinetic energy K; in the Monte Carlo scheme we employed, described in Appendix A, K is in fact employed as a “demon”, allowing us to keep the total energy E constant.

3.2. Characterization of the Phase Transition

Preliminary to our analysis of the geometry and topology of the PLS, we show here the that a phase transition is indeed occurring, and try to determine its precise nature.

To this end, we examined quantities that are usually expected to exhibit singular behavior at the transition.

3.2.1. Specific Heat, Caloric Curve and Entropy Derivatives

The specific heat $c_v$ typically displays these critical behaviors in most phase transitions; this can be due to the presence of latent heat, in the case of first order phase transition, or to critical fluctuations, in the case of continuous phase transitions.

In the microcanonical ensemble, the specific heat can be computed according to

$$c_v={\left({\displaystyle \frac{dT}{dE}}\right)}^{-1},$$

(13)

where $T={\displaystyle \frac{2\langle K\rangle }{3N}}$ is the kinetic temperature. Alternatively, we also use the results of [13], which used the Laplace-transform techniques to propose a variety of alternative definitions for usual thermodynamical observables. Among three different formulas for $c_v$, we only display one here, as they all yield the same results, up to the accessible precision.

$$c_v={\displaystyle \frac{3}{2}}{\left[1-{\displaystyle \frac{3N}{2}}\left({\displaystyle \frac{\langle K^2\rangle }{{\langle K\rangle }^2}}-1\right)\right]}^{-1}.$$

(14)

The comparison of the curves obtained with both Eq. (13) and (14) is commonly employed as an equilibration test (see for instance [17]).

Inspection of Figure 1 shows two clear peaks of the specific heat. Except for a few data points, in particular around the second transition point, the two sets of points coincide within their error margin. Given the difficulty of equilibration in glass-formers, especially in the presence of critical fluctuations, we deem this result quite satisfying.

In the following graphs, we flag the estimated positions of the two peaks of Figure 1 with two vertical dotted lines. We denote ${ϵ}_1$ and ${ϵ}_2$ the critical energy density of the first and the second peak, respectively.

Figure 2 makes it clear that, to ${ϵ}_2$ corresponds a plateau of the caloric curve, indicating the occurrence of latent heat in correspondence with an internal arrangement of the Lennard-Jones mixture; the higher energy transition is thus a first order transition. On the other hand, to ${ϵ}_1$ corresponds only a slight inflection of the curve, suggesting a second order phase transition.

Ehrenfest classification of order transitions relies on the loss of analyticity of Helmoltz free energy. However, the relevant thermodynamic potential in the microcanonical ensemble is the entropy, which is widely regarded as a quantity of deeper physical and mathematical meaning. Yet, after (13), the microcanonical specific heat can be rewritten as

$$c_v\left(ϵ\right)=-{\left({\displaystyle \frac{\partial S}{\partial E}}\right)}^2{\left({\displaystyle \frac{{\partial }^{2}S}{\partial E^2}}\right)}^{-1},$$

(15)

emphasizing that the observed singular behavior of $c_v$ can, in principle, find its origin in a divergence of the first order derivative of the entropy, or in the vanishing of its second order derivative. Furthermore, while, in the canonical ensemble, the average specific energy $\langle ϵ\rangle \left(T\right)$ usually displays clear critical behaviors at the transition temperature, the microcanonical ensemble inverse temperature $\beta \left(ϵ\right)={\displaystyle \frac{1}{T\left(ϵ\right)}}$ is often much less sensitive.

Motivated by these observations, in Refs. [12,18,19,20,21], novel methods of classification of phase transitions in the microcanonical ensemble were proposed, relying on the analysis of inflection points of the derivatives of the entropy. In fact, in the absence of a phase transition, all derivatives of S of even order are strictly concave, and those of odd order strictly convex.

In Figure 3, we show the two lowest order derivatives of S with respect to $ϵ$, namely

$$\begin{array}{cc}\hfill \beta \left(ϵ\right)& ={\displaystyle \frac{\partial S}{\partial E}}\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill \gamma \left(ϵ\right)& ={\displaystyle \frac{{\partial }^{2}S}{\partial E^2}}={\displaystyle \frac{\partial \beta }{\partial E}}\hfill \end{array}$$

(17)

To retrieve these quantities, we departed from the kinetic temperature T straightforwardly obtained from simulation, and differentiated $\beta =1/T$ with respect to the energy.

$\beta$ remains convex on the whole range of energies considered, but exhibits visible backbendings at ${ϵ}_1$ and ${ϵ}_2$. $\gamma \left(ϵ\right)$ shows local maxima at both transition point. Around ${ϵ}_2$, it is very close to zero, even reaching the positive region if the error is taken into account. This suggests again the occurrence of a first order phase transition at this critical energy, according to the classification of [20]. The critical behavior of $c_v=-{\beta }^2/\gamma$ thus evidently originates from $\gamma$ approaching zero.

These results are remarkably similar to those displayed in [19,21] for the solid-solid and solid-liquid transition of an elastic flexible polymer.

3.2.2. Translational and Orientational Order

Local ordering emerges at low energies, as exemplifies the projective views of Figure 4, where structures appear at low ($ϵ<{ϵ}_1$) and intermediate (${ϵ}_1<ϵ<{ϵ}_2$) energies. Depending on the projection axis, aperiodic pentagonal arrangement or square-shaped seemingly periodic structures are visible, in both these regimes. However, such rough observations are very unreliable, because the visibility of ordered structures is highly dependent on the chosen projection axis. Moreover, despite this apparent order, our samples are far from crystallization and still exhibit a high degree of disorder that makes it very inconvenient for the human eye to distinguish particular geometries.

The observation of a local order is confirmed by the profiles of the radial partial pair distribution functions, displayed in Figure 5. At high energy, these functions are in good agreement with the results obtained in [16] for $T=0.39$, corresponding to an energy just above ${ϵ}_2$. At low energy, the first peak of the three functions becomes sharper, and a few other bumps appear at larger distances in $g_{11}\left(r\right)$ and $g_{22}\left(r\right)$. we observe the appearance of several new peaks of the density profile. This is a clear signature of the emergence of translational order. Ill-formed periodic profiles emerge at low energy, in particular in $g_{1}1\left(r\right)$, indicating partial crystallization.

To further characterize the nature of these configurational changes, we inspect the bond-orientational order parameters $Q_l$, first defined in [22] to characterize crystalline order in Lennard-Jones liquids. For a given $l\in \mathbb{N}$, it writes

$$Q_l=\sqrt{{\displaystyle \frac{4\pi }{2l+1}}\sum _{m=-l}^l{\left|Q_{lm}\right|}^2},$$

(18)

with

$$Q_{lm}={\displaystyle \frac{1}{n_B}}\sum _{(i,j)\in B}{Y}_{lm}\left(\theta \left({\mathbf{r}}_{ij}\right),\phi \left({\mathbf{r}}_{ij}\right)\right),$$

(19)

where B is the considered set of bonds, $n_B$ its cardinality, $\theta \left({\mathbf{r}}_{ij}\right)$ and $\phi \left({\mathbf{r}}_{ij}\right)$ are respectively the azymutal and polar angles of the bond vector ${\mathbf{r}}_{ij}$ in a fixed reference frame, and the $Y_{lm}$ are spherical harmonics.

Two particles $i,j$ are considered bonded if $r_{ij}<c_b$, where $c_b$ is an arbitrary cutoff. As is often prescribed [23,24,25], we set $c_b$ to be the approximate position of the second minimum of the radial distribution functions right after the first peak (displayed in Figure 5).

The authors of Ref. [16] showed that this model exhibits a short to medium-range order, namely a local tetrahedral ordering, coined as a tetrahedral network. The results displayed in Figure 6 seem to corroborate this observation, as the order parameters defined in Eq. (18) exhibit clear steps at the transition energies ${ϵ}_1$, ${ϵ}_2$, for certain values of l.

Interestingly, at energy density ${ϵ}_1$, most of the $Q_l$ exhibit singular behaviors, whether it be a positive of a negative peak, suggesting a temporary rearrangement of the particles upon cooling.

The profiles of the functions $Q_l\left(ϵ\right)$ clearly indicate the emergence, at low energy, of some kind of orientational order, that is the repetition, throughout the whole sample, of some preferred angles between neighboring atoms. However, we found no correspondence between the combination of values we found for the $Q_l$, and that of a well-defined known crystalline structure, as described e.g. in [22].

We thus conclude that our samples underwent frustrated crystallization, with different types of arrangement competing in the route to minimize the potential, most likely, considering the mapping given in [22], a mixing between fcc, hcp and icosahedral orders. This kind of behavior has been proposed to be a mechanism at the origin of the glass transition [11]. This competition could indeed explain the dramatic dynamical slowing down associated to this class of transitions, and that we observed in this study, in the form of a tremendously long equilibration time.

3.3. Topological Changes

Finally, in correspondence with these two transitions, we also found important changes of the geometrical properties of the PLS.

Figure 7 shows, in correspondence with the two $c_v$-peaks, inflection points of the total mean curvature of the submanifold ${\Sigma }_{\varphi }$, indication of a change of the landscape of this hypersurface. These inflection points are clear discontinuities, a sharp corner in ${ϵ}_1$, and a step in ${ϵ}_2$.

The average variance of the principal curvatures, shown in Figure 8, exhibit clear changes of slope at these transition points, providing a strong indication of a change of the topology of ${\Sigma }_{\varphi }$, according to Pinkall’s theorem (5). Namely, such a change is necessarily due to a change in the values of the Betti numbers, hence of the topological properties of ${\Sigma }_{\varphi }$; though the precise nature of these changes is not accessible to our analysis, we can expect that, from the high energy chaotic phase to the low energy glass phase, ${\Sigma }_{\varphi }$ loses connectivity and the system is more easily confined to restricted regions of state-space.

Finally, the variance of the scalar curvature, shown in Figure 9, jumps in ${ϵ}_1$, exhibiting a wide peak in the intermediary region, and a second, smaller peak in ${ϵ}_2$. As per Overholt’s theorem (6), this quantity is an upper bound to the alternate sum of the Betti numbers. These sharp changes are thus indications of major topological changes of the PLS.

All of these observations suggest that there are indeed important changes in the topology of ${\Sigma }_{\varphi }$ at play in correspondence with these two critical points of $c_v$. In this intermediary region, the overall shape of the manifold ${\Sigma }_{\varphi }$ most probably undergo dramatic changes. These changes correspond to the rearrangement of particle configurations in pseudo-crystalline orders upon cooling.

4. Discussion

In this work, we addressed the phenomenon of glass transition using a novel approach, considering its relation with topological changes of the PLS. Using a Monte Carlo algorithm, we studied on the equilibrium properties of a glass-forming system.

Although the glass transition is often considered a fundamentally dynamical phenomenon, we gathered several indications that the system, despite having reached equilibrium, did not fully crystallize and yet exhibits clear signatures of phase transitions. This is due to the frustrated nature of the model, a Lennard-Jones binary mixture with competing interactions, that lacks a well-defined crystalline phase at the studied density.

To overcome ergodicity breaking, we employed advanced numerical techniques that allowed the system to “jump” between confining energy minima. This, combined with the absence of a crystalline phase in this frustrated system, enabled us to perform equilibrium measurements of the glass phase, based on importance sampling of the state space.

Through microcanonical analysis of the entropy derivatives, we identified two distinct transitions: a low-energy, second-order transition and a higher-energy, first-order transition. We conjecture that the first order transition is in fact a glass-liquid transition, while the second order one is a reconfiguration of the glass into different local orderings. This is supported by our finding of important changes of the radial pair distribution function and of the bond-orientational order parameters at the transitions points. These results do not correspond to what could be expected in a crystal, and are instead in agreement with recent studies [10,11] suggesting that the glass transition can be understood in terms of frustration on the way to crystallization.

Finally, in agreement with the topological theory of phase transition, we observed that several quantities that are deeply linked with the topology of the PLS exhibit inflection points and discontinuities at the transition energies. This study introduces a new approach in our understanding of the phenomenon of vitrification.

Future perspectives include the repetition of our experiment with different system sizes, to be able to perform finite-sized scaling, confirming our findings and calculating critical exponents. It would also be of interest to identify crystalline seeds in the material, and perform localized measurements of the same quantities investigated in the present work, to compute spatial correlation functions. Doing so, spatial correlation functions could be computed, which we expect to present divergent behavior around the transition energies (see [11]. Furthermore, it would reveal to which extent the topological analysis in this context is sensible to the scaling. In other words, whether the well-known degeneracy of the energy minima in the glass phase is an inherently global phenomenon, and if the system finds locally well-defined energy minima.