SciPost Phys. 10, 001 (2021)

Spatial structure of unstable normal modes
in a glass-forming liquid
Masanari Shimada1? , Daniele Coslovich2,3 , Hideyuki Mizuno1 and Atsushi Ikeda1
1 Graduate School of Arts and Sciences, The University of Tokyo, Tokyo 153-8902, Japan
2 Laboratoire Charles Coulomb (L2C), Université de Montpellier, CNRS, Montpellier, France
3 Dipartimento di Fisica, Università di Trieste, Strada Costiera 11, 34151, Trieste, Italy? masanari-shimada444@g.ecc.u-tokyo.ac.jp

Abstract
The phenomenology of glass-forming liquids is often described in terms of their underlying, high-dimensional potential energy surface. In particular, the statistics of stationary
points sampled as a function of temperature provides useful insight into the thermodynamics and dynamics of the system. To make contact with the real space physics,
however, analysis of the spatial structure of the normal modes is required. In this work,
we numerically study the potential energy surface of a glass-forming ternary mixture.
Starting from liquid configurations equilibrated over a broad range of temperatures using a swap Monte Carlo method, we locate the nearby stationary points and investigate
the spatial architecture and the energetics of the associated unstable modes. Throughthis spatially-resolved analysis, originally developed to study local minima, we corroborate recent evidence that the nature of the unstable modes changes from delocalized to
localized around the mode-coupling temperature. We find that the displacement amplitudes of the delocalized modes have a slowly decaying far field, whereas the localized
modes consist of a core with large displacements and a rapidly decaying far field. The
fractal dimension of unstable modes around the mobility edge is equal to 1, consistent
with the scaling of the participation ratio. Finally, we find that around and below themode-coupling temperature the unstable modes are localized around structural defects,
characterized by a disordered local structure markedly different from the liquid’s locallyfavored structure. These defects are similar to those associated to quasi-localized vibrations in local minima and are good candidates to predict the emergence of localized
excitations at low temperature.
Copyright M. Shimada et al.
This work is licensed under the Creative Commons
Attribution 4.0 International License.
Published by the SciPost Foundation.

Received 18-09-2020
Accepted 15-12-2020
Check forupdatesPublished 04-01-2021
doi:10.21468/SciPostPhys.10.1.001

Contents
1 Introduction

2

2 Methods
2.1 Sample preparation
2.2 Normal mode analysis

4
4
5
1

SciPost Phys. 10, 001 (2021)

2.3

2.2.1 Participation ratio
2.2.2 Decay profile
2.2.3 Fractal dimension
2.2.4 Energy profile
2.2.5 One-particle dynamical matrix
Local structure
2.3.1 Locally favored structures
2.3.2 Structural order parameter

5
6
6
6
7
7
7
7

3 Results
3.1 Spatial structure of unstable modes
3.2 Local structure of saddles and local minima
3.3 Local structure of unstable and stable cores

8
8
12
14

4 Conclusions

17

A Difference between id and ie

18

B Local minima

19

C Size dependence

19

References

20

1

Introduction

Glass formation is driven by the rapid increase of structural relaxation times that takes placewhen liquids are cooled fast enough to avoid crystallization [1, 2]. From a purely operational
point of view, the glass transition occurs at a temperature, called Tg , at which the structural
relaxation time reaches about 100s. Below Tg , the supercooled liquid behaves, on typical experimental time scales, as an amorphous solid. However, the lack of a sharp divergence indynamic and thermodynamic quantities, as would be observed in conventional phase transitions, makes it difficult to define the glassy state in a clear-cut way, and the mechanisms
responsible for the enormous slowing down of the dynamics are still highly debated [3–5].
A long-standing line of thought describes glass formation in terms of the exploration of ahigh-dimensional potential energy surface (PES) [6, 7]. Stationary points of the PES are configurations such that the gradient of the total potential energy U vanishes, and correspond toeither local minima or saddles. These configurations play an important role in the PES-based
description of glass formation. In his seminal work in 1969, Goldstein [8] argued that the
dynamics of supercooled liquids is dominated by activated barrier crossing between neighboring local minima below some temperature Tx , at which the structural relaxation time is oforder 10−9 s. In this scenario, Tx marks the crossover between the normal liquid dynamics andactivated glassy dynamics over large energy barriers.
This phenomenological description later evolved into a more complex picture, which suggests the existence of an additional dynamic regime [7]. The early stages of the slowing downof several supercooled liquids are indeed reasonably well described by mode-coupling theory (MCT) [9]. MCT provides a semi-quantitative description of several non-trivial dynamic
features above the critical temperature TMCT , at which the theory predicts a power law diver2

SciPost Phys. 10, 001 (2021)

gence of the relaxation times. This singularity is smeared out in actual supercooled liquids by
the presence of thermally activated processes, which are not accounted for by the theory andturn the transition into a crossover. According to a series of theoretical studies [10, 11], however, a sharp change in the PES still underlies the MCT crossover: above TMCT , the dynamics
is governed by the motion along unstable saddle modes, which the system can exploit without resorting to thermal activation; below TMCT , instead, the system mostly vibrates around
local minima, occasionally relaxing via activated barrier crossing, and one can thus identify
TMCT ≈ Tx . This theoretical scenario is inspired by the behavior of a class of mean-field spinmodels, for which the (schematic) MCT equations are exact and the nearest stationary point
to an equilibrium configuration indeed changes from a saddle to a minimum at the MCT transition [12–14].
Finding a clear signature of the MCT crossover in the PES of actual supercooled liquids
is quite challenging and this has led to contradictory claims [15–25]. These difficulties were
due to two main factors. First, conventional simulations cannot easily sample the PES atequilibrium below TMCT , due to exceedingly long equilibration times. Second, all the above
mentioned studies neglected the role of the modes’ spatial localization. In particular, it isclear that thermally activated processes in finite dimensional liquids involve only localized
portions of the system, whereas the unstable modes of saddles in mean-field models are spatially delocalized. A recent numerical study [26], based on an efficient swap Monte Carloalgorithm [27,28], has tackled these two issues at once, establishing that the stationary pointsof several three-dimensional model liquids indeed show a sharp change around TMCT : the fraction of delocalized unstable modes vanishes at the MCT crossover temperature. At any temperature, however, saddles also possess a finite fraction of localized unstable modes, which only
involve a finite number of particles, as originally envisaged by Goldstein. The extent to which
the MCT transition is avoided in actual supercooled liquids thus depends on the concentrationof such localized excitations [26].
Despite these advances, the spatial structure of the unstable modes remains largely unknown. In fact, the simple finite size scaling analysis of Ref. [26] was only based on the bulk
average mode participation and further investigation is needed to characterize the modes’
structure. Moreover, since the saddle modes progressively stabilize at low temperature, it isnatural to inquire the connection between them and the stable modes populating the lowfrequency portion of the vibrational spectrum. The lowest-frequency modes of local minima
display quasi-localized vibrations (QLVs) [29, 30], which are at present much better understood than the unstable saddle modes [29–37]. QLVs consist of an energetically unstable core
and a stable far-field [38] and their eigenenergy is determined by the competition between
these two components. The vibrational amplitudes decay with a distinct power law r −(d−1) ,where r is the distance from the core and d is the spatial dimension, unless hybridization
with acoustic modes occurs [29, 33, 39]. Moreover, the QLVs spatial structure brings useful
insight into several local phenomena in low temperature glasses, such as the response to a
local dipolar force or plastic events associated to shear-transformation zones [40–45].
In this work, we apply the spatially resolved mode analysis developed to study QLVs [29,38] to account for the saddle modes spatial structure. To achieve this, we improved the optimization protocol used in Ref. [26], which suffered from poor convergence when searchingfor stationary points in moderately large system sizes. Through this analysis, we confirm thatthe delocalized unstable modes that populate the saddles’ spectrum at high temperature havean extended character, at least up to the length scale we could probe. We also show that the
localized unstable modes are unique to the saddle structure. Their cores consist of only a few
particles displaying very large displacements, while their far field decays rapidly compared tothat of delocalized modes and of the QLVs. In local minima, truly localized vibrations would
only be found in the high-frequency tail of the spectrum. However, we provide evidence that3

SciPost Phys. 10, 001 (2021)

around and below TMCT localized unstable modes originate from structural defects similar to
those associated to the cores of the QLVs, or “soft spots”, which have been identified in previous
studies on metallic glasses [46]. Overall, our study casts a bridge between so far disconnected
investigations of the PES and provides a spatially resolved picture of the progressive stabilization of saddle modes as a supercooled liquid turns into a glass.
This paper is organized as follows. In Sec. 2.1, we introduce our model and describe how
we located saddles and local minima starting from equilibrium configurations. In Sec. 2.2
and 2.3.1, we define the quantities used to characterize these stationary points. We provide
our main results in Sec. 3. We focus on the spatial structure of unstable saddle modes inSec. 3.1 and the local structure of unstable cores in Sec. 3.3. We compare saddles and local
minima in Sec. 3.2. Finally, we conclude our work with a summary in Sec. 4.

2

Methods

2.1

Sample preparation

We present numerical results for the ternary mixture model introduced by Gutiérrez et al. [47].
The model consists of three species of particles, which interact via a repulsive inverse power
potential with exponent 12
uαβ (r) = ε

σ

αβ

r

12

+ c4

σ

αβ

r

−4

+ c2

σ

αβ

r

−2

+ c0 ,(1)

where α, β = A, B, C are species indices. The additional terms in Eq.(1) ensure continuity of thepotential up to its second derivative at the cutoff distance rcut = 1.25σαβ . See Ref. [48] for theσAAσexpressions of c0 , c2 , and c4 . The interactions are additive, the size ratio is σBB= σCBBC = 1.25
and the chemical compositions are x A = 0.55, x B = 0.30, and x C = 0.15. In the following, we
will express energies and distances in units of ε and σAA, respectively.
The equilibrium configurations we used for this study were obtained using swap Monte
Carlo simulations in Ref. [26]. In the following, we will focus on systems composed of
N = 3000 particles at a number density ρ = 1.1 and temperatures ranging from 0.45 to
0.28. Note that a few samples crystallized at T = 0.28 and were therefore removed from theanalysis. Results for N = 1000 are shown in Appendix C. The mode-coupling critical temperature, estimated in Ref. [49] by fitting the structural relaxation times to a power law, is
TMCT ≈ 0.29. Starting from these equilibrium configurations, we performed three kinds of
optimizations: (i) potential energy minimizations using the l-BFGS method [50] to locate thelocal minima of the PES; (ii) total square-force minimizations (or W -minimizations) using thel-BFGS method to locate local minima of
N

W=

1X
| F~i |2 ,
N i=1

(2)

where F~i is the force acting on particle i; (iii) eigenvector-following (EF) optimizations [51]to locate stationary points of the PES with a prescribed number nu of unstable modes.
Our goal is to locate stationary points, either local minima or saddles, in the neighborhoodof a given equilibrium configuration. In both cases, stationary points are identified as points
for which W is smaller than 10−10 . As is well known, W -minimizations mostly locate so-called
quasi-saddles, i.e., local minima of W with a finite value of the total force and precisely one
inflection mode, while only a small fraction of W -minimizations reach true stationary points.
By contrast, the EF method is guaranteed to converge to a stationary point of prescribed ordernu , but the value of nu has to be provided as input. To define a mapping between a given
4

SciPost Phys. 10, 001 (2021)

configuration and its neighboring stationary point we have followed the approach used inRef. [26]: for each starting configuration, we first perform a W -minimization. The order nu of
the W -minimized configuration is then used as the target order for an EF optimization. All thesaddles analyzed in Sec. 3 were obtained using this protocol and converged to the required
tolerance, W < 10−10 .
The EF optimizations carried out in Ref. [26] displayed poor convergence for the systemsize of interest in this work (N = 3000). One possible reason is that the algorithm gets stuck insome restricted portion of configuration space, without being able to reach the neighborhoodof the stationary point, where convergence is fast. To help the algorithm escape faster fromthose regions, in this work we have allowed for larger steps in configuration space. In the
EF optimization method, the size of the steps is limited by a set of trust radii [51], which
ensure a local harmonic approximation at each step. To mitigate the convergence issues, wehave therefore used a larger tolerance on the deviation from harmonic approximation, namely
we increased (decreased) the trust radii by 20% if the relative error on the corresponding
eigenvalue was smaller (larger) than 200%. The threshold was 100% in Ref. [26]. Also, we
set the initial trust radius to 0.1, instead of 0.2. Empirically, we found that using a larger
tolerance substantially improved the convergence rate for large system sizes. The fraction of
optimizations that converged to a stationary point after 4000 iterations ranged from 100%
at T = 0.28 to 56% at T = 0.45. The fractions obtained using a threshold of 100% ranged
from 48% at T = 0.28 to 3% at T = 0.35. Of course, using a large tolerance sometimes led to
substantial overlaps between particles and thus large values of the potential energy. To prevent
numerical issues, we avoided steps such that the potential energy of the new configuration waslarger than a threshold (104 ), and decreased all the trust radii until the potential energy of the
new configuration dropped below the threshold. Except for these small but important practicaldetails, the algorithm was the same as the one used in Ref. [26].

2.2

Normal mode analysis

We performed a standard normal mode analysis [52] for local minima and saddles. We diagonalized the dynamical matrix, which is the Hessian of the total potential energy U,¨−M i j(i 6= j)
∂ 2U
Hi j == PN,(3)∂ ~ri ∂ ~r j
k=1 Mik (i = j)
where

Mi j = u00i j (ri j )~ni j n~ iTj +

u0i j (ri j ) 
ri j


I −n
~i j n~ iTj ,

(4)is the contribution from the pair 〈i j〉, ui j (ri j ) is the pair interaction potential, n~ i j = ~ri j /ri j is
the unit vector along the pair, and I is the d × d identity matrix. We denote its eigenvaluesby λα and the corresponding eigenvectors by ~eα = (~eα,1 · · · ~eα,N ), where α = 1, 2, . . . , 3N − 3.
We sort them in ascending order as λ1 < λ2 < . . .. Note that we always excluded the threeeigenmodes corresponding to global translations. Using these eigenmodes we computed thequantities defined below.
2.2.1

Participation ratio

The participation ratio [53–55] of a mode α is given by
P(λα ) = PN

1eα,i |4
i=1 |~.

(5)

It quantifies the degree of localization of a given mode. When all particles have equal displacements, P(λα ) = N , while P(λα ) = 1 when the mode is localized on a single particle.
5

SciPost Phys. 10, 001 (2021)

2.2.2

Decay profile

The decay profile [29, 39] dα (r) of a mode α is a function of the distance r from the particle
with the largest |~eα,i |, which we denote by id . The decay profile is defined asdα (r) =medianri ∈[r−∆r/2,r+∆r/2]|~eα,i | ,(6)where ri is the distance of a particle i from the particle id . In the case of saddles, we show theaveraged decay profile for a certain eigenvalue λ defined asd(r) = 〈dα (r)〉λ ,(7)

where 〈•〉λ denotes the average over all the eigenmodes with eigenvalues
λα ∈ [λ − ∆λ/2, λ + ∆λ/2]. For local minima, we averaged the data over the lowest frequency mode of each sample, which are almost always QLVs for N ® 104 [29].
2.2.3

Fractal dimension

To further characterize the spatial structure of the eigenmodes, we introduce another functionof the norms |~eα,i |. We identify the particles contributing to a mode α as those with the dP(λα )elargest norms, where dxe denotes the ceiling function, which returns the least integer greaterthan or equal to x, and then determine the number Nα (r) of such contributing particles up to a
distance r from the particle id . We define the averaged function N (r) by the same procedure asin Eq. (7). From this function, we can estimate the fractal dimension D of mode of eigenvalue
λ from N (r) ∼ r D .
2.2.4

Energy profile

We define the local vibrational energy of a particle i on mode α as
N
T

1X
~eα,i − ~eα, j Mi j ~eα,i − ~eα, j .
2 j=1


N
u0i j (ri j )1 X 00⊥2
ni j · ~eα,i j )2 +,=u (ri j )(~(~eα,ij)
2 j=1 i jri j

δEα,i =

(8)

where ~eα,i j = ~eα,i − ~eα, j is the relative displacement between the pair and⊥
2(~eα,ieα,i j )2 − (~ni j · ~eα,i j )2 is the squared transverse relative displacement. The energyj ) = (~
profile [38] Λα (r) of a mode α is a function of the distance r from the particle with the smallest δEα,i , which we denote by ie . Λα (r) is defined asXΛα (r) =δEα,i ,
(9)ri <rwhere ri is the distance of a particle i from the particle ie . This measures the vibrational energyof the system when a given mode is excited, and we do not include the total potential energy ofsaddles or minima. In the following, we show the averaged normalized energy profile definedas〈Λα (r)〉λΛ̃(r) =,
(10)〈λα 〉λThe distance r used here differs slightly from the one used in the paper by Gartner and Lerner [39]. They
defined the center of an eigenmode using w particles with the largest |~eα,i |2 and measured the distance r from thiscenter. They mainly used w = 4 while our definition corresponds to w = 1. This difference does not affect the ourresults qualitatively.

6

SciPost Phys. 10, 001 (2021)

for
P saddles, while we used only the lowest frequency QLVs for local minima. Sincei δEα,i = λα , Λ̃(r) → −1 as r → ∞ when λ < 0 while Λ̃(r) → 1 when λ > 0. Finally,
we note that the difference between id and ie does not matter in practice, because there is anegative correlation between |~eα,i | and δEα,i (see Appendix A).
2.2.5

One-particle dynamical matrix

Finally, using Mi j in Eq. (4), we define the one-particle dynamical matrix as(1)
Hi j = δi j

N
X
k=1

Mik .

(11)

(1)

Clearly, the full dynamical matrix reduces to Hi j when only a tagged particle is allowed tomove and the others are kept fixed. We denote its eigenvalues by µβ and the eigenvectors by
f~β , where β = 1, 2, . . . , N d.

2.3

Local structure

2.3.1

Locally favored structures

To provide further insight into the localization features of the modes and their link to the
local structure, we performed a radical Voronoi tessellation of the minimized configurations
using the Voro++ software [56]. The radical Voronoi tessellation requires as input the typical
diameters rα of the particles of type α, which we estimated from the positions of the first peaksof the partial radial distribution functions gαα (r) measured for minimized configurations at
T = 0.28, namely r1 = 1.34, r2 = 1.11, r3 = 0.93. Following a well-established approach [57],we characterized the shape the Voronoi polyhedra using the signature (n3 , n4 , n5 , . . . ), where
nk is the number of faces of the polyhedron having k vertices. We found that the (0,0,12)
signature, corresponding to an icosahedral local arrangement, becomes the most frequent at
low temperature, in both minimized and instantaneous configurations. Following Ref. [24],we therefore identify the icosahedron as the locally favored structure of the model.
2.3.2

Structural order parameter

To provide further insight into the local structure of our ternary mixture model, we calculated the Tong-Tanaka structural order parameter Θ [58]. This order parameter measures theaverage local deviation from close packing of neighboring particles, and has been shown tocorrelate quite strongly with the local dynamics in some models of supercooled liquids [58].
For the calculation of Θ we followed the method described in Ref. [58], except for thefollowing minor modification. Since the contact radius between two neighboring particles isdefined in terms of the nominal radii (σii +σ j j )/2, in our soft sphere model negative deviations
from ideal packing can occur. Negative deviations are harmless, as they only imply a linear
shift of the order parameter, provided the absolute value in Eq. 2 of Ref. [58] is removed. We
checked that this minor modification, which we used in this work, does not affect our calculation. In particular, at T = 0.28, reducing the contact radii by about 10% with respect to theirnominal values effectively removes all negative contributions to Θ and lead to qualitatively
similar results. Finally, we note that as in Ref. [58], we identified the network of neighbors
using a radical Voronoi tessellation, see Sec. 2.3.1.

7

SciPost Phys. 10, 001 (2021)

3

Results

3.1

Spatial structure of unstable modes
Table 1: Range of eigenvalues used to compute the average in Eq. (7), mobility
edge λe , and fraction of the delocalized unstable modes nd /(3N ) at all investigatedtemperatures. Note that a few isolated eigenvalues are included in the averages at
T = 0.28, 0.39, 0.32.

P (λα )

103

T

−λ

−λe

nd /3N

0.45

1.25–23.75

10.9

0.015

0.35

1.25–21.25

4.92

0.005

0.32

1.25–21.25, 36.25

2.71

0.00180.30

1.25–18.75, 26.25

1.40

0.000 48

0.28

1.25–26.25, 31.25––(a) T = 0.45λe(b) T = 0.35λe(c) T = 0.32(d) T = 0.30λe(e) T = 0.28
λe

102
101
100
−30

−20

−10
λα

0

−20

−10
λα

0

−30 −20 −10λα0−30 −20 −10λα0−30 −20 −10λα0Figure 1: Scatter plots of the participation ratio for saddles at (a) T = 0.45, (b)T = 0.35, (c) T = 0.32, (d) T = 0.30, and (e) T = 0.28. In each panel, the solid anddotted lines correspond to λ = 0 and to the mobility edge λe , respectively.
In this section, we quantitatively characterize the spatial structure of the unstable saddle
modes across the MCT crossover temperature. We show the corresponding results for the QLVs
of the local minima in Appendix B, but some of the data of the QLVs are included in this section
for comparison.
To briefly summarize the results of Ref. [26] using our data, we show the scatter plots of theparticipation ratio in Fig. 1. We show the results at (a) T = 0.45, (b) T = 0.35, (c) T = 0.32,(d) T = 0.30, and (e) T = 0.28. From a finite-size scaling analysis of the participation ratio,
we can define the mobility edge λe that separates localized and delocalized modes. It is the
fixed point in P(λ, N )N −1/3 , where P(λ, N ) is the average participation ratio of modes witheigenvalue λ in a system of size N . The modes below the mobility edge are localized because
their participation ratio is generally independent of N , whereas the modes above the mobility
edge are delocalized because their participation grows, on average, at least linearly with thelinear system size N 1/3 . This definition of the mobility edge was originally given in the context of the Anderson localization [59] and later adapted to the study of instantaneous normal
modes by Clapa et al. [60]. The mobility edges determined in Ref. [26] using the above procedure are shown as dotted lines in Fig. 1; note that λe vanishes around TMCT ≈ 0.29. Table 4shows the specific values of λ and λe we used to perform the average in Eq. (7). We used abin width ∆λ = 2.5 and ignored the bins containing only one mode.
Before proceeding to a quantitative analysis of our data, we visualize the real space structure of some selected unstable modes, to grasp the main difference between localized and8

SciPost Phys. 10, 001 (2021)

(a)

(b)

(c)

(f )
(d)

(e)

Figure 2: The unstable modes of a saddle at T = 0.45. We show dP(λα )e particles
with the largest norms shown by black arrows. The particle with the largest norm
is placed at the center in each box. Different colors indicate different particle types.
We scale the norms by factors 20 in (a), (b), (c), and (d), and 40 in (e) and (f) for
visualization purposes. The parameters of each mode are as follows: (a) α = 1,λα = −28.7, and dP(λα )e = 2, (b) α = 11, λα = −13.3, and dP(λα )e = 15, (c)α = 21, λα = −10.6, and dP(λα )e = 45, (d) α = 41, λα = −7.72, and dP(λα )e = 49,(e) α = 81, λα = −4.79, and dP(λα )e = 275, and (f) α = 151, λα = −0.722, anddP(λα )e = 815.
delocalized modes. In Fig. 2 we show six modes of a typical saddle configuration at T = 0.45.
In each box, we only show the dP(λα )e particles having the largest norms. Different colors
indicate different particle types. The particle with the largest norm is placed at the center ofthe box and the displacements are shown by black arrows. We scale the norms by a factor 20
in (a), (b), (c), and (d), and 40 in (e) and (f) for visualization purposes. Details of the parameters used for each mode are given in the caption of Fig. 2. Here, we note that the modes
shown in (a) and (b) are localized and that the others are delocalized. We can clearly see thecores of the localized modes in Fig. 2(a) and (b), while it is difficult to identify similar cores inthe delocalized modes, at least by visual inspection. The snapshots of Fig. 2 are also suggestive of the presence of a percolation transition as P increases. In a first attempt to address thisquestion, we have analyzed the network of the [P] particles which, for a given mode, have thelargest participation ratios. Our preliminary results are qualitatively compatible with a continuum percolation at a threshold participation ratio in the range 300 − 400, but the present
system size is obviously too small to draw firm conclusions. We leave this interesting point toa future investigation.
Let us start our analysis by looking at the decay profiles. We remind that the decay profilesof QLVs in local minima display a r −2 scaling away from the cores [29], which is the same as
the response of the elastic body to a local force dipole. We group the modes according to theireigenvalues as described in Sec. 2.2.2 and analyze the same set of temperatures as in Fig. 1.
The results are shown in Fig. 3. The continuous change in color from blue to red corresponds
to the change of λ shown in Tab. 4: the data above the mobility edge in red and those below
the mobility edge in blue. For comparison, we also plot the data for the QLVs in black. In ourdata for the QLVs, we cannot observe a clear asymptotic behavior d(r) ∝ r −2 because our

9SciPost Phys. 10, 001 (2021)

10−1
∝ r−2

∝ r−2

∝ r−2

∝ r−2

∝ r−2

d(r)

10−2
10−3
10−4

1

λ→0

λ→0

λ→0

λ→0

λ→0

(a) T = 0.45(b) T = 0.35(c) T = 0.32(d) T = 0.30(e) T = 0.28

2

3

4 5 612r

3

4 5 612r

3

4 5 612r

3

4 5 61

2

r

3

4 5 6

r

Figure 3: Decay profiles for saddles at (a) T = 0.45, (b) T = 0.35, (c) T = 0.32, (d)T = 0.30, and (e) T = 0.28. The solid lines are proportional to r −2 . The continuouschange in color from blue to red corresponds to the change of λ shown in Tab. 4. We
show the data above the mobility edge in red, those below the mobility edge in blue.
For comparison, we show the data the QLVs in black.

N (r)/N10−110(a) T = 0.45(b) T = 0.35(c) T = 0.32(d) T = 0.30(e) T = 0.28

−2

10−3
10−4

1

2

3
r

4 5 6 12

3

4 5 6 12r

3

4 5 6 12r

3r

4 5 6 1

2

3

4 5 6

r

Figure 4: The fraction of particles N (r)/N up to a distance r for unstable modes ofsaddles. All parameters are the same as in Fig. 3, but we do not show the data with
λ < −21.25 at T = 0.45, λ < −16.25 at T = 0.35, and λ < −13.75 at T ≤ 0.32 forvisualization purposes. The dashed and solid lines show power laws with the fractal
dimensions D = 1 and 2, respectively.
systems are too small for this purpose † . Nonetheless, all the panels show the same tendency:d(r) decreases more rapidly than r −2 at large r for λ  λe , while d(r) tends to saturate at
finite values at large r for λ  λe . The former implies that the unstable modes much below
the mobility edge are definitely more localized than the QLVs of minima, while the latter thatthe unstable modes much above the mobility edge are definitely delocalized, which cannot
be explained as the response of the elastic body to a local disturbance. Close to the mobility
edge (λ ∼ λe ), unstable modes display d(r) ∼ r −2 in a wide region of r. However, the current
system size is too small to clearly identify this scaling as the asymptotic behavior. To achievethis, we would need to study much bigger systems, which are beyond the reach of this work.
The delocalized modes can be better characterized by N (r)/N , which allows us to identify
the average fractal exponent of the modes. We show these results in Fig. 4. All parameters are
the same as in Fig. 3, but for visualization purposes we do not show the data with λ < −21.25
at T = 0.45, λ < −16.25 at T = 0.35, and λ < −13.75 at T ≤ 0.32. The dashed and solid lines
show power laws with the fractal dimensions D = 1 and 2, respectively. This figure shows thatthe modes with λ ∼ λe (shown in gray) follow N (r) ∼ r 1 . This is consistent with the definitionof the mobility edge, at which P grows linearly with L. We note, however, that N (r) is anaveraged quantity and there are strong mode-to-mode fluctuations (not shown). Thus, weshould not expect that each mode at the mobility edge has a string-like structure, as one may
naively guess. This is also clear from Fig. 2(c), which shows a mode at the mobility edge having
a complex, though somewhat ramified, spatial structure. As for the other modes, the curves ofthe localized modes rapidly converge to their final values while the delocalized modes and the†

In Ref. [29], systems of 106 particles were needed to establish this scaling.

10

SciPost Phys. 10, 001 (2021)

0

0

0

(a) T = 0.45
λ→0

Λ̃(r)

−0.5

0

(b) T = 0.35
−0.5(d) T = 0.30

L/2
0

4

8
r

−1.5
12
0

(e) T = 0.28

−1

−1

−1

−1

0

(c) T = 0.32

−2

−2

4

8

12

−2

0

4

r

8
r

12

−3

−4
0

4

8
r

12

0

4

8

12

r

Figure 5: Energy profiles for saddles. All parameters are the same as in Fig. 3. Weshow the data of the QLVs only in (e). The dashed vertical lines show half of the box
length L/2.
QLVs grow faster than linear with r. Nevertheless, we can also see a small difference between
the delocalized unstable modes and the QLVs, particularly at T ≥ 0.35: the latter saturate atlarge distances. Although subtle, this difference shows that the particles contributing to the
QLVs are denser and more compact than those contributing to the delocalized modes. Thisobservation is consistent with the results of the decay profiles.
To further characterize the mode structure, we analyze the energy profiles, see Fig. 5.
The parameters are the same as in Fig. 3, but we show the data of the QLVs only at T = 0.28
because they display very little temperature dependence (see Appendix B). The dashed verticallines are placed at half of the box length L/2. We find that localized and delocalized unstablemodes display qualitatively different energy profiles. The former have pronounced dips atr ∼ 2, which are the cores of the modes [32], and for r ¦ 2 they rapidly converge to −1 frombelow. Thus, the localized modes have weak far-field components and their energy is mostly
determined by the cores. This is consistent with the results of the decay profile and differs
markedly from the QLVs, whose energy is determined by the competition between the unstablecore and the stable far-field. We also note that the energy profiles of the localized modes
do not depend on the temperature above the MCT crossover temperature. By contrast, thedelocalized modes have the large far-field components and display a pronounced temperaturedependence. The energy profiles of the delocalized modes also show a small dip (or kink) atr ∼ 2, which indicates that even these modes may possess a core. However, at T ≥ 0.32 weobserve Λ̃(r ∼ 2) > −1. Therefore, the energy profiles converge to −1 from above. In otherwords, both the core and the far field are unstable. This is significantly different from thelocalized modes and also from the QLVs.
For T ≤ 0.30, i.e. below the localization transition, the energy profiles converge to −1from below and the functional form is similar to that of the QLVs. This indicates that thesoftest localized unstable modes that survive below the mode-coupling crossover [26] havespatial properties similar to the soft stable modes, see also Sec. 3.3. We emphasize, however,that the energy profiles discussed here characterize only the average spatial features of thevibrational modes. Dips or kinks at r ∼ 2 in the energy profiles indicate the presence ofregions that are particularly unstable on average. Our analysis does not tell, however, howmany cores a given mode possesses. Actually, since delocalized modes have unstable far fields
at high temperatures, one may speculate that they even have multiple cores. Only when themode is strongly localized, as observed in Fig. 2(a) and in Fig. 6(a), the particle ie corresponds
to the unique core of the mode. Finally, we expect that gathering the cores of the QLVs should
provide similar information as the “soft spots” [40], although the two definitions differ inpractice.
To better understand the results of the energy profiles, we visualize in Fig. 6 the spatial
distribution of the local energy δEα,i of the unstable modes. We used the same modes as inFig. 2 and the particles with the most negative δEα,i are placed at the center of the squares.

11

SciPost Phys. 10, 001 (2021)

(a)

(b)

(c)

(d)

(e)

(f )

Figure 6: Spatial distribution of the local energy δEα,i for the same modes as in Fig. 2.
Each snapshot represents a slab of width 1.5. The particles with the most negativeδEα,i are placed at the center of the squares. The particles with darker colors havemore negative energy, while particles with positive energy are colored white.
Particles having more negative energy are colored with darker colors, while those with positiveenergy are colored white, regardless of the value of δEα,i . We can see a few dark particles atthe centers of the boxes in all cases; they correspond to the small dips at r ∼ 2 observed inFig. 5 (a). The localized modes with small P(λα ) have energetically stable backgrounds, i.e.,the background is uniform white, while the delocalized ones with large P(λα ) have unstablebackgrounds, i.e., there are scattered yellow regions. This is consistent with the differences
between localized and delocalized modes mentioned above.
3
Saddles
Minima

2.5

g(r)

2
1.5
1
0.5
0

(a) T = 0.45

11.5(b) T = 0.28

2

1

r

1.5

2

r

Figure 7: Radial distribution function for saddles and local minima of N = 3000.
The temperatures are (a) T = 0.45 and (b) T = 0.28.

3.2

Local structure of saddles and local minima

In the previous section, we showed that the localized modes of saddles are generally more
localized than the QLVs. Since the interaction potential is the same, it is plausible to expectthat this reflects a structural difference between saddles and local minima. Otherwise, it may

12

SciPost Phys. 10, 001 (2021)

10−2

PDF

10−3
10−4
Saddles
Minima

10−5
10−6
10−7

(a) T = 0.45
0

150

300
µ

(b) T = 0.28

450

600 0

150

300
µ

450

600

Figure 8: PDF of the eigenvalues of the one-particle dynamical matrix. All parameters
are the same as in Fig. 7.
100

100

(a) T = 0.45(b) T = 0.28
200

10−3

200

10−2

µmax,αµmax,α10−2maxβ (~eα · f~β )2

10−1

maxβ (~eα · f~β )2

10−1

10−3

0
10−4
100

101

102

0
10−4
100

103

P (λα )

101

102

103

P (λα )

Figure 9: Left axis: the largest overlap between the eigenvectors of the full dynamical
matrix and the one-particle dynamical matrix maxβ (~eα , f~β )2 . Right axis: the corresponding eigenvalue µmax,α . The abscissa is the participation ratio. We show the data
for saddles at (a) T = 0.45 and (b) T = 0.28.
also be that unstable and stable modes are intrinsically different. To investigate this point, we
computed the radial distribution function g(r) for for saddles and local minima, see Fig. 7.
From inspection of these functions, we cannot see any difference between saddles and local
minima at low temperature, T = 0.28. By contrast, the first peaks of local minima shift outwards compared to that of saddles at T = 0.45. Thus, some particles of saddles are closer toeach other than particles of local minima. The very existence of unstable modes in the saddles
can be partly attributed to these structural differences. This is because the closer two particles,
the larger the force, and as is evident from the second term in Eq. (4), a large repulsive force
significantly decreases the eigenenergy of the modes. We note that such destabilization plays
a dominant role also for QLVs [38]. In addition, the dip at r ∼ 1.4 for local minima shifts
slightly downwards at T = 0.45.
To further address this point, we characterize the local structure of the local minima andsaddles using a Voronoi tessellation, see Sec. 2.3.1. As explained in that section, we characterized the local arrangements around each particle using the signature of corresponding Voronoi
polyhedron, as is routinely done to identify the locally favored structures of supercooled liquids [24]. We found that the Voronoi signatures in saddles and minima have different statistics
at T = 0.45. Local minima having a sensibly larger fraction of icosahedral structures, which
we identified as the locally favored structure of the model. By contrast, the statistics of signatures become very similar for minima and saddles sampled at low temperature (T = 0.28),
see Tab. 2, and the fraction of the icosahedral structures is only marginally larger in minima.
This confirms that overall minima are locally more structured than saddles at high T , but thatthese structural differences fade away close to TMCT . The overall structural similarity between
saddles and minima at low temperature suggests that the localized unstable modes find their
origin in subtle features of the fluid structure.

13

SciPost Phys. 10, 001 (2021)

Table 2: Percentage of most frequent Voronoi signatures in saddles and local minima
sampled at T = 0.28. The ? symbol denotes the signature with the largest occurrence
for core particles.
T = 0.28
Saddles

Local minima(0,0,12)17.9(0,0,12)18.9(0,2,8,2)

7.2(0,2,8,2)

7.4(0,1,10,2)

6.1(0,1,10,2)

6.3

?

3.8

(0,2,8,1)

3.8(0,3,6,4)?

(0,2,8,1)

3.8(0,3,6,4)

3.8(0,2,8,4)

3.0(0,2,8,4)

3.0(0,1,10,4)

2.3(0,1,10,4)

2.5(0,3,6,3)

2.1(0,1,10,3)

2.1(0,2,8,5)

2.0(0,3,6,3)

2.1(0,1,10,3)

2.0(0,2,8,5)

2.1

Table 3: Percentage of most frequent Voronoi signatures in saddles and local minima
sampled at T = 0.45.
T = 0.45
Saddles

3.3

Local minima(0,2,8,2)

4.7(0,2,8,2)

7.2(0,3,6,4)

4.4(0,3,6,4)

5.3(0,2,8,1)

3.3(0,2,8,1)

4.7(0,3,6,3)

3.2(0,0,12)4.6(0,1,10,2)

2.5(0,1,10,2)

4.2(0,2,8,4)

2.4(0,3,6,3)

3.5(0,0,12)1.9(0,2,8,4)

3.2(0,3,6,5)

1.6(0,2,8,5)

1.9(1,2,6,3,1)

1.6(0,3,6,5)

1.8(1,2,5,4)

1.5(0,2,8,3)

1.7

Local structure of unstable and stable cores

One may expect that the unstable modes are localized around some sort of structural defects
in the supercooled liquid, in analogy to what found in local minima [40, 46]. To address thisquestion, we characterize the structure of the unstable modes’ cores in three different ways:
by analyzing the one-particle dynamical matrix, by computing restricted few-body correlation
functions, and by analyzing the Voronoi cells surrounding the core particles.
The one-particle dynamical matrix is the dynamical matrix obtained when only a tagged
particle is allowed to move and the others are kept fixed (see Sec. 2.2.5 for details). Here
we show that its eigenmodes give us an insight into the cores of the localized unstable modesof saddles. In Fig. 8 we show the probability distribution functions (PDF) of its eigenvalues.
All the parameters are the same as in Fig. 7. We can see a clear difference between saddles
14

SciPost Phys. 10, 001 (2021)

(a)

(b)

0.04

3

b(θ)

g(r)

4

2

0.02

Saddles
Unstable cores
QLV cores

1
0

1

2

3

0

0

60

r

120

180

θ

Figure 10: (a) Radial distribution function g(r) and (b) bond-angle distribution functions b(θ ) for small particles (black), for small particles at the core of unstable modes(orange) for saddles sampled at T = 0.28, and for small particles at the core of the
QLVs (symbol) for local minima at the same temperature.
and local minima at T = 0.45, as the PDF of saddles has broader tails than that of local
minima. Namely, some particles of saddles feel much steeper potentials in certain directions
than particles of local minima. Thus, the structure of saddles is more inhomogeneous.
The localized unstable modes with small P(λα ) have significant correlations with the eigenmodes of the one-particle dynamical matrix. Here, we consider the overlap between the eigenvectors of H and H(1) , i.e., (~eα · f~β )2 , to investigate whether the motions of the localized unstable modes can be explained at the one-particle level. Figure 9 shows the maximum overlap
maxβ (~eα · f~β )2 (left axis) and the corresponding eigenvalue µmax,α (right axis) defined aswhereµmax,α = µβmax,α ,(12)βmax,α = arg max(~eα · f~β )2 .

(13)

β

The abscissa is the participation ratio P(λα ) of the corresponding eigenmode of H. We show
the data of saddles at (a) T = 0.45 and (d) T = 0.28. This figure shows that the modes withsmaller P(λα ) have larger overlap with the eigenmodes of H(1) . Therefore, the motions of thestrongly localized modes can actually be explained at the one-particle level. The correlation
between µmax,α and P(λα ) is not strong, but µmax,α tends to be negative when P(λα ) < 2.
To gain more insight into the structural features of the cores, we analyze selected few-bodycorrelation functions. In this analysis, we focus on the local environment around particles of
species C. This is motivated by the fact that the vast majority (≈ 90%) of the core particlesof the localized unstable modes belong to this species. Therefore, we restrict the calculationof the radial distribution function g(r) to central particles of species C that form the coresof the unstable modes, i.e., the particles whose index is ie . Note that any other neighbors
at a distance r is used in this calculation, irrespective of its species. Similarly, we computethe restricted bond-angle distribution b(θ ) obtained from the angles formed between a central particle of species C belonging to the cores and two of its nearest neighbors, irrespectiveor their species [61]. Note that the bond-angle distributions are normalized such that the
distribution is flat when angles are drawn randomly on a sphere. We emphasize that these
correlation functions do not entail information of the correlations between the cores but on thelocal structure around the core particles.
We perform our analysis at T = 0.28, where only localized unstable modes are present
and the identification of the cores is therefore unambiguous. Figure 10 reveals that the local
structure around the core particles is markedly different from the average. In particular, theg(r) of the core particles displays a very flat first minimum, suggestive of a much less structured
first coordination shell. The absence of the splitting of the second peak of g(r) means that core15

SciPost Phys. 10, 001 (2021)

20

Saddles
Unstable cores

16

PDF

12
8
4
0

0

0.1

0.2

0.3

Θ

Figure 11: Probability density functions of the Θ order parameter for all particles inthe saddles sampled at T = 0.28 (black line) and for particles belonging to the coresof the unstable modes (orange line).
particles are unlikely to align linearly with their neighbors [62], which is consistent with theabsence of a peak at 180 degrees in the bond-angle distribution. Overall, the structure aroundcore particles is almost featureless and qualitatively resembles the one of the fluid at a higher
temperature. In particular, only at temperatures much higher than T = 0.45 the equilibriumg(r) would display a first minimum as flat as the one found for core particles.
We also analyzed the signatures of the Voronoi cells surrounding the core particles andfound that they never correspond to the locally favored structure of the model, i.e., the icosahedron, which accounts for about 18% of the Voronoi signatures at T = 0.28, see Tables 2
and 3. The Voronoi cells around core particles are characterized instead by a broad distribution of signatures, not corresponding to any symmetric or clearly identified structural arrangement. The most frequent signature around core particles of unstable modes is (0,3,6,4)(4.6%), followed by (0,3,6,6) (2.6%), (0,4,4,4) (2.1%) and several others with similarly small
concentrations. Interestingly, the cores of the QLVs in local minima display similar structuralfeatures: from Fig. 10 we see that their g(r) and b(θ ) resemble the ones of the unstable cores
discussed above. Moreover, we found that the most frequent signature around core particles
is the same for unstable modes and QLVs in local minima, albeit with a higher concentration(7.9%) in the latter. The remaining most frequent signatures around core particles in unstablemodes and QLVs differ. Overall, the results indicate that, around and below the MCT crossover,
localized unstable modes and soft stable modes have a similar structural origin.
To better identify the structural features of the unstable modes’ cores, we calculated the Θ
structural order parameter introduced in Ref [58], which measures the deviations from ideal
close packings of spheres (see Sec. 2.3.2). We carried out this analysis at T = 0.28. We foundthat Θ is significantly larger for particles belonging to the unstable modes’ cores than for thebulk, as demonstrated by the shift in the corresponding probability density p(θ ), see Fig. 11.
Therefore, the cores are associated to less compact local environments. We also calculated theprecision with which Θ can identify the cores, as follows. Say we have nc core particles in agiven saddle configuration. We then select the nc particles having the largest Θ values. The
fraction of particles these subsets have in common defines the precision with which the coresare predicted by Θ. We found that the precision is higher than would be obtained by picking
particles at random, but still relatively low (< 20%). These preliminary results highlight thedifficulty in pinpointing defects in glassy configurations. This problem may be partly attributed
to the presence of spurious second shell neighbors in the radical Voronoi tessellation.
From this analysis we conclude that it is indeed appropriate to identify the core particles
as “structural defects” in the liquid, in agreement with related observations about soft spots
16

SciPost Phys. 10, 001 (2021)

in granular packings [40] and in a model metallic glass [46]. We note, however, that while
the locally favored structure is the same in our model and in the one of Ref. [46], the Voronoi
signatures of the structural defects have a very broad, model-specific distribution. Thus, thiskind of Voronoi-based analysis does not clearly tell us what structural defects are – at best what
they are not. Measurements of the Tong-Tanaka structural order parameter indicate that, in
our model, unstable cores are localized around steric defects associated to loose local packings.
However, it remains difficult to predict which particles will be associated to the cores. Devisinga more general, unsupervised approach to identify structural defects is crucial to predict plastic
events in glasses under shear [63] and also dynamic heterogeneities in supercooled liquids [58,61, 64].

4

Conclusions

In this paper we have provided a spatially-resolved analysis of the structure of the unstable saddle modes of a ternary glass-former across the mode-coupling crossover. Our analysis provides
an independent confirmation of the findings of Ref. [26] concerning the presence of two kindsof unstable modes: (i) delocalized modes, which are characterized by spatially extended displacement fields and which disappear below the mode-coupling crossover temperature TMCTand (ii) localized modes, whose properties are largely independent of temperature. We havealso performed a detailed comparison to the features of the softest stable modes observed inlocal minima. This extended analysis was possible thanks to some simple but crucial tweaks
to the parameters of the EF optimization algorithm, which enabled us to study larger systemsizes than those of Ref. [26].
Our results confirm that the disappearance of delocalized unstable modes below the MCT
crossover is the only trace of the sharp change of the PES observed in mean-field glass models.
By contrast, localized unstable modes find their origin in structural defects of the liquid, andare therefore a genuine finite-dimensional feature. The fact that these two kinds of modes
coexist in a reasonably realistic model glass-former supports the view that the PES is a usefulframework to reconcile mean-field and real space pictures of glasses. In particular, it would be
worth investigating the connection between localized unstable modes and the localized excitations invoked in dynamic facilitated models and observed in some model glass-formers [65].
See also Ref. [66] for a recent attempt along similar lines using instantaneous normal modes.
These two kinds of modes are separated by the mobility edge λe < 0; the localized (delocalized) modes are found below (above) it, and λe → 0 as T → TMCT . Our independent analysis
confirms the presence of a qualitative change of the spatial localization features across themobility edge. We found that modes at the edge have an average fractal dimension of one,consistent with the finite size scaling analysis of Ref. [26]. However, we warn that one should
not interpret this result in a naif geometric sense: unstable modes at the mobility edge have avery open structure but are not clearly string-like and the N (r) ∼ r scaling only holds on average. Our results call for a more in-depth characterization of the mode structure at a per-modelevel.
We further compared the unstable modes of saddles to the QLVs that populate the soft
portion of the vibrational spectrum of local minima. We found that at high temperature thelocal structure of local minima and saddles differ significantly, and the decay profiles andfractal dimensions of QLVs do not clearly match the ones of either delocalized or localized
unstable modes. However, as the temperature drops around and below TMCT , local minima
and saddles become structurally very similar and QLVs share some features with the softest
unstable modes, i.e., those with lowest absolute frequency. Our analysis of the energy profiles
and of the structure of the modes’ cores indicates that at low temperature unstable modes and17

SciPost Phys. 10, 001 (2021)

QLVs are localized around similar structural defects, but in the former the elastic response of
the medium surrounding the core is not strong enough to stabilize the system. In the studied
model, the unstable modes’ cores are associated to steric defects, corresponding to loose local
packings, localized outside the regions of local icosahedral order. Although the computationalcosts are far beyond our reach for the moment, it would also be useful to use much larger
systems and investigate the precise functional forms of the decay and energy profiles of thesaddle modes.
In conclusion, the PES of a simple but realistic model glass-former always possesses purely
localized saddle modes, whose displacements and energies are concentrated on cores composed of few particles, i.e., they are defects unique to the saddle structure. By contrast, delocalized unstable modes, characterized by a slower decay of the displacement amplitudes andenergy profiles, are only accessible at temperatures above the MCT crossover. At low temperature, soft unstable modes and QLVs share similar spatial features and they tend to be localizedaround similar structural defects. Our study connects previous investigations of the PES, whichhave focused either on high-order stationary points or on local minima, and provides a real
space picture of how an equlibrium liquid characterized by several unstable modes converges
to a mechanically stable glass.

Acknowledgments
This work was supported by JSPS KAKENHI Grants No. 18H05225, 19J20036, 19H01812,
19K14670, 20H01868, 20H00128 and partially supported by Initiative on Promotion of Supercomputing for Young or Women Researchers, Supercomputing Division, Information Technology Center, The University of Tokyo. DC acknowledges support as a JSPS International
Research Fellow.

A

Difference between id and ie
1

δEα,i

0
−2

0

−4
−1

−6
−8

(a) Saddles(b) Minima
−2

10−1

100
|~eα,i |10−1100|~eα,i |Figure 12: Scatter plots of |~eα,i | vs. δEα,i at T = 0.32 for (a) saddles and (b) local
minima .
The origins of the functions d(r) and Λ(r) are different. However, it was shown in Ref. [38]using the QLVs of nearly jammed packings that there is a negative correlation between |~eα,i | andδEα,i . We confirm this correlation in the present case. Figure 12 shows the data at T = 0.32.
Different colors in (a) indicate different eigenvalues shown in Tab. 4. To average |~eα,i | and18

SciPost Phys. 10, 001 (2021)

δEα,i , we sorted |~eα,i | in descending order |~eα,1 | > |~eα,2 | > · · · > |~eα,N | and averaged |~eα,i | andδEα,i with fixed i.

Local minima
10−1
(a)

d(r)

∝ r−2

(b)

0
−2

10−2

T = 0.45
10−2 T = 0.35
T = 0.32
T = 0.30
T = 0.28
1

(c)

Λ̃(r)

10−1

N (r)/N

B

−4

L/2

−6

2

3

4 5 6

10−3

1

2

r

3

4 5 6

−8
0

4

8

r

12

r

Figure 13: (a) Decay profile, (b) N (r)/N , and (c) energy profile for the QLVs of IS.
The solid line in (a) is proportional to r −2 . The dashed and solid lines in (b) show
the power laws with the fractal dimensions D = 1 and 2, respectively. The dashedvertical line in (c) shows half of the box length L/2.
To compare with the results of saddles, we show the results of local minima. The decay
profiles, N (r)/N , and the energy profiles are shown in Fig. 13(a), (b), and (c), respectively.
Different symbols and colors indicate different temperatures. Unlike the case of saddles, weaveraged them over the lowest-frequency QLVs in each sample. We cannot see any strong temperature dependence, however the modes are slightly more localized when the temperature
decreases.

C

Size dependence
Table 4: Eigenvalues used to compute the average in Eq. (7), mobility edge λe , andfraction of the delocalized unstable modes nd /(3N ) at all investigated temperatures
for N = 1000 particles.
T

−λ

−λe

nd /3N

0.45

1.25, 6.25, 11.25, 16.25, 21.25

10.9

0.015

0.35

1.25, 6.25, 11.25, 16.25, 21.25

4.92

0.005

0.32

1.25, 6.25, 11.25, 16.25

2.71

0.00180.30

1.25, 6.25, 11.25, 16.25

1.40

0.000 480.28

1.25, 6.25, 11.25, 16.25––

To study the system size dependence, we show the data for N = 1000. Figure 14 shows thecorresponding decay profile. We show the data for N = 1000 by lines and those for N = 1000
by symbols. The specific eigenvalues used to compute the average are shown in Tab. 4. We canhardly see any size dependence. By contrast, the energy profiles do depend on system size,
see Fig. 15. In particular, the delocalized unstable modes show a strong size dependence athigh temperatures while the localized modes do not depend on N . This is natural because the19

SciPost Phys. 10, 001 (2021)

10−1

d(r)

10−2
10−3
10−4

λ→0
(a) T = 0.45
123(b) T = 0.35
4 5 612r(c) T = 0.32

3

4 5 612

3r(d) T = 0.30
4 5 612(e) T = 0.28

3r

4 5 61

2

r

3

4 5 6

r

Figure 14: Finite size effects on the decay profile for saddles at (a) T = 0.45, (b)T = 0.35, (c) T = 0.32, (d) T = 0.30, and (e) T = 0.28. We show the data of
N = 1000 by lines and those of N = 3000 by symbols. The continuous change incolor from blue to red corresponds to the change of λ shown in Tab. 4. We show thedata above the mobility edge in red, those below the mobility edge in blue.
000(a) T = 0.45

Λ̃(r)(d) T = 0.30

0

4

8
r

−1.5

12

−1

4

8

12

−2

−2

−2
0

(e) T = 0.28

−1

−1
−1

0

(c) T = 0.32

−0.5

−0.5

0

(b) T = 0.35

0

4

r

8
r

12

−3

0

4

8
r

12

−4

0

4

8

12

r

Figure 15: Finite size effects on the energy profile for saddles. All parameters are thesame as in Fig. 14.
energy of the delocalized modes is determined by the entire system while the core dominates
the energy of the localized modes. Note that, however, the qualitative behaviors do not dependon N in any case.

References[1] C. A. Angell, Formation of glasses from liquids and biopolymers, Science 267, 1924 (1995),
doi:10.1126/science.267.5206.1924.
[2] L. Berthier and M. D. Ediger, Facets of glass physics, Phys. Today 69, 40 (2016), doi:
10.1063/pt.3.3052.
[3] L. Berthier and G. Biroli, Theoretical perspective on the glass transition and amorphousmaterials, Rev. Mod. Phys. 83, 587 (2011), doi:10.1103/RevModPhys.83.587.
[4] M. Wyart and M. E. Cates, Does a growing static length scale control the glass transition?,
Phys. Rev. Lett. 119, 195501 (2017), doi:10.1103/PhysRevLett.119.195501.
[5] L. Berthier, G. Biroli, J.-P. Bouchaud, and G. Tarjus, Can the glass transition be explained without a growing static length scale?, J. Chem. Phys. 150, 094501 (2019), doi:
10.1063/1.5086509.
[6] A. Heuer, Exploring the potential energy landscape of glass-forming systems: From inherent
structures via metabasins to macroscopic transport, J. Phys.: Condens. Matter 20, 373101(2008), doi:10.1088/0953-8984/20/37/373101.

20

SciPost Phys. 10, 001 (2021)[7] A. Cavagna, Supercooled liquids for pedestrians, Phys. Rep. 476, 51 (2009), doi:
10.1016/j.physrep.2009.03.003.
[8] M. Goldstein, Viscous liquids and the glass transition: A potential energy barrier picture, J.
Chem. Phys. 51, 3728 (1969), doi:10.1063/1.1672587.
[9] W. Götze, Complex dynamics of glass-forming liquids: A mode-coupling theory, Oxford
University Press (2009).
[10] A. Cavagna, Fragile vs. strong liquids: A saddles-ruled scenario, Europhys. Lett. 53, 490(2001), doi:10.1209/epl/i2001-00179-4.
[11] T.-S. Grigera, Geometrical properties of the potential energy of the soft-sphere binary mixture, J. Chem. Phys. 124, 064502 (2006), doi:10.1063/1.2151899.
[12] J. Kurchan and L. Laloux, Phase space geometry and slow dynamics, J. Phys. A: Math. Gen.
29, 1929 (1996), doi:10.1088/0305-4470/29/9/009.
[13] A. Cavagna, I. Giardina, and G. Parisi, Stationary points of the Thouless-Anderson-Palmer
free energy, Phys. Rev. B 57, 11251 (1998), doi:10.1103/PhysRevB.57.11251.
[14] A. Cavagna, I. Giardina and G. Parisi, Role of saddles in mean-field dynamics above the glasstransition, J. Phys. A: Math. Gen. 34, 5317 (2001), doi:10.1088/0305-4470/34/26/302.
[15] K. Broderix, K. Bhattacharya, A. Cavagna, A. Zippelius, and I. Giardina, Energy landscapeof a Lennard-Jones liquid: Statistics of stationary points, Phys. Rev. Lett. 85, 5360 (2000),
doi:10.1103/PhysRevLett.85.5360.
[16] L. Angelani, R. Di Leonardo, G. Ruocco, A. Scala and F. Sciortino, Saddles in the energy landscape probed by supercooled liquids, Phys. Rev. Lett. 85, 5356 (2000), doi:
10.1103/PhysRevLett.85.5356.
[17] T. Grigera, A. Cavagna, I. Giardina, and G. Parisi, Geometric approach to the dynamic glasstransition, Phys. Rev. Lett. 88, 055502 (2002), doi:10.1103/PhysRevLett.88.055502.
[18] L. Angelani, R. Di Leonardo, G. Ruocco, A. Scala and F. Sciortino, Quasisaddles as relevant
points of the potential energy surface in the dynamics of supercooled liquids, J. Chem. Phys.
116, 10297 (2002), doi:10.1063/1.1475764.
[19] L. Berthier and J.-P. Garrahan, Real space origin of temperature crossovers in supercooled
liquids, Phys. Rev. E 68, 041201 (2003), doi:10.1103/PhysRevE.68.041201.
[20] J. P. K. Doye and D. J. Wales, Saddle points and dynamics of Lennard-Jones clusters, solids,
and supercooled liquids, J. Chem. Phys. 116, 3777 (2002), doi:10.1063/1.1436470.
[21] J. P. K. Doye and D. J. Wales, Comment on “Quasisaddles as relevant points of the potential
energy surface in the dynamics of supercooled liquids”, J. Chem. Phys. 118, 5263 (2003),
doi:10.1063/1.1553754.
[22] D. J. Wales and J. P. K. Doye, Stationary points and dynamics in high-dimensional systems,
J. Chem. Phys. 119, 12409 (2003), doi:10.1063/1.1625644.
[23] B. Doliwa and A. Heuer, Energy barriers and activated dynamics in a supercooled LennardJones liquid, Phys. Rev. E 67, 031506 (2003), doi:10.1103/PhysRevE.67.031506.

21

SciPost Phys. 10, 001 (2021)[24] D. Coslovich and G. Pastore, Understanding fragility in supercooled Lennard-Jones
mixtures. II. Potential energy surface, J. Chem. Phys. 127, 124505 (2007), doi:
10.1063/1.2773720.
[25] M. Ozawa, W. Kob, A. Ikeda, and K. Miyazaki, Equilibrium phase diagram of arandomly pinned glass-former, Proc. Natl. Acad. Sci. USA 112, 6914 (2015), doi:
10.1073/pnas.1500730112.
[26] D. Coslovich, A. Ninarello and L. Berthier, A localization transition underlies the mode-coupling crossover of glasses, SciPost Phys. 7, 077 (2019), doi:
10.21468/SciPostPhys.7.6.077.
[27] D. Gazzillo and G. Pastore, Equation of state for symmetric non-additive hard-sphere fluids:
An approximate analytic expression and new Monte Carlo results, Chem. Phys. Lett. 159,
388 (1989), doi:10.1016/0009-2614(89)87505-0.
[28] T. S. Grigera and G. Parisi, Fast Monte Carlo algorithm for supercooled soft spheres, Phys.
Rev. E 63, 045102 (2001), doi:10.1103/PhysRevE.63.045102.
[29] E. Lerner, G. Düring, and E. Bouchbinder, Statistics and properties of low-frequency
vibrational modes in structural glasses, Phys. Rev. Lett. 117, 035501 (2016), doi:
10.1103/PhysRevLett.117.035501.
[30] H. Mizuno, H. Shiba and A. Ikeda, Continuum limit of the vibrational properties of amorphous solids, Proc. Natl. Acad. Sci. USA 114, E9767 (2017), doi:
10.1073/pnas.1709015114.
[31] E. Lerner and E. Bouchbinder, Effect of instantaneous and continuous quenches on thedensity of vibrational modes in model glasses, Phys. Rev. E 96, 020104 (2017), doi:
10.1103/PhysRevE.96.020104.
[32] M. Shimada, H. Mizuno and A. Ikeda, Anomalous vibrational properties in the continuum
limit of glasses, Phys. Rev. E 97, 022609 (2018), doi:10.1103/PhysRevE.97.022609.
[33] G. Kapteijns, E. Bouchbinder and E. Lerner, Universal nonphononic density ofstates in 2D, 3D, and 4D glasses, Phys. Rev. Lett. 121, 055501 (2018), doi:
10.1103/PhysRevLett.121.055501.
[34] L. Wang, A. Ninarello, P. Guan, L. Berthier, G. Szamel and E. Flenner, Low-frequency
vibrational modes of stable glasses, Nat. Commun. 10, 26 (2019), doi:10.1038/s41467018-07978-1.
[35] E. Lerner, Finite-size effects in the nonphononic density of states in computer glasses, Phys.
Rev. E 101, 032120 (2020), doi:10.1103/PhysRevE.101.032120.
[36] P. Das, H. George E. Hentschel, E. Lerner and I. Procaccia, Robustness of densityof low-frequency states in amorphous solids, Phys. Rev. B 102, 014202 (2020), doi:
10.1103/PhysRevB.102.014202.
[37] D. Richard, K. González-López, G. Kapteijns, R. Pater, T. Vaknin, E. Bouchbinder and E. Lerner, Universality of the nonphononic vibrational spectrum acrossdifferent classes of computer glasses, Phys. Rev. Lett. 125, 085502 (2020), doi:
10.1103/PhysRevLett.125.085502.

22

SciPost Phys. 10, 001 (2021)[38] M. Shimada, H. Mizuno, M. Wyart and A. Ikeda, Spatial structure of quasilocalized
vibrations in nearly jammed amorphous solids, Phys. Rev. E 98, 060901 (2018), doi:
10.1103/PhysRevE.98.060901.
[39] L. Gartner and E. Lerner, Nonlinear plastic modes in disordered solids, Phys. Rev. E 93,
011001 (2016), doi:10.1103/PhysRevE.93.011001.
[40] M. L. Manning and A. J. Liu, Vibrational modes identify soft spots in a sheared disorderedpacking, Phys. Rev. Lett. 107, 108302 (2011), doi:10.1103/PhysRevLett.107.108302.
[41] L. Yan, E. DeGiuli and M. Wyart, On variational arguments for vibrational modes near
jamming, Europhys. Lett. 114, 26003 (2016), doi:10.1209/0295-5075/114/26003.
[42] E. Lerner and E. Bouchbinder, A characteristic energy scale in glasses, J. Chem. Phys. 148,
214502 (2018), doi:10.1063/1.5024776.
[43] M. Shimada, H. Mizuno and A. Ikeda, Vibrational spectrum derived from local mechanicalresponse in disordered solids, Soft Matter 16, 7279 (2020), doi:10.1039/D0SM00376J.
[44] C. Rainone, E. Bouchbinder and E. Lerner, Pinching a glass reveals key properties of its soft
spots, Proc. Natl. Acad. Sci. USA 117, 5228 (2020), doi:10.1073/pnas.1919958117.
[45] C. Rainone, E. Bouchbinder and E. Lerner, Statistical mechanics of local force
dipole responses in computer glasses, J. Chem. Phys. 152, 194503 (2020), doi:
10.1063/5.0005655.
[46] J. Ding, S. Patinet, M. L. Falk, Y. Cheng and E. Ma, Soft spots and their structural
signature in a metallic glass, Proc. Natl. Acad. Sci. USA 111, 14052 (2014), doi:
10.1073/pnas.1412095111.
[47] R. Gutiérrez, S. Karmakar, Y. G. Pollack and I. Procaccia, The static lengthscale characterizing the glass transition at lower temperatures, Europhys. Lett. 111, 56009 (2015),
doi:10.1209/0295-5075/111/56009.
[48] S. Karmakar, E. Lerner, I. Procaccia and J. Zylberg, Effect of the interparticle potential on the yield stress of amorphous solids, Phys. Rev. E 83, 046106 (2011), doi:
10.1103/PhysRevE.83.046106.
[49] A. Ninarello, L. Berthier, and D. Coslovich, Models and algorithms for the next generation ofglass transition studies, Phys. Rev. X 7, 021039 (2017), doi:10.1103/PhysRevX.7.021039.
[50] D.-C. Liu and J. Nocedal, On the limited memory BFGS method for large scale optimization,
Math. Program. 45, 503 (1989), doi:10.1007/BF01589116.
[51] D. J. Wales, Rearrangements of 55-atom Lennard-Jones and (C60 )55 clusters, J. Chem. Phys.
101, 3750 (1994), doi:10.1063/1.467559.
[52] C. Kittel, Introduction to solid state physics, Wiley (2005).
[53] V. Mazzacurati, G. Ruocco and M. Sampoli, Low-frequency atomic motion in a model glass,
Europhys. Lett. 34, 681 (1996), doi:10.1209/epl/i1996-00515-8.
[54] H. R. Schober and G. Ruocco, Size effects and quasilocalized vibrations, Philos. Mag. 84,
1361 (2004), doi:10.1080/14786430310001644107.
[55] S. N. Taraskin and S. R. Elliott, Anharmonicity and localization of atomic vibrations invitreous silica, Phys. Rev. B 59, 8572 (1999), doi:10.1103/PhysRevB.59.8572.
23

SciPost Phys. 10, 001 (2021)[56] C. H. Rycroft, VORO++: A three-dimensional Voronoi cell library in C++, Chaos 19,
041111 (2009), doi:10.1063/1.3215722.
[57] M. Tanemura, Y. Hiwatari, H. Matsuda, T. Ogawa, N. Ogita and A. Ueda, Geometricalanalysis of crystallization of the soft-core model, Progr. Theor. Phys. 58, 1079 (1977), doi:
10.1143/PTP.58.1079.
[58] H. Tong and H. Tanaka, Revealing hidden structural order controlling both fast andslow glassy dynamics in supercooled liquids, Phys. Rev. X 8, 011041 (2018), doi:
10.1103/PhysRevX.8.011041.
[59] B. I. Shklovskii, B. Shapiro, B. R. Sears, P. Lambrianides and H. B. Shore, Statistics ofspectra of disordered systems near the metal-insulator transition, Phys. Rev. B 47, 11487(1993), doi:10.1103/PhysRevB.47.11487.
[60] V. Iulian Clapa, T. Kottos and F. W. Starr, Localization transition of instantaneous normal
modes and liquid diffusion, J. Chem. Phys. 136, 144504 (2012), doi:10.1063/1.3701564.
[61] J. Paret, R. L. Jack and D. Coslovich, Assessing the structural heterogeneity of supercooled liquids through community inference, J. Chem. Phys. 152, 144502 (2020), doi:
10.1063/5.0004732.
[62] A. S. Clarke and H. Jónsson, Structural changes accompanying densification of random
hard-sphere packings, Phys. Rev. E 47, 3975 (1993), doi:10.1103/PhysRevE.47.3975.
[63] D. Richard et al., Predicting plasticity in disordered solids from structural indicators (2020),
arXiv:2003.11629.
[64] E. Boattini, S. Marín-Aguilar, S. Mitra, G. Foffi, F. Smallenburg and L. Filion, Autonomously revealing hidden local structures in supercooled liquids, Phys. Rev. Mater. 4,113609 (2020), doi:10.1103/PhysRevMaterials.4.113609.
[65] A. S. Keys, L. O. Hedges, J. P. Garrahan, S. C. Glotzer and D. Chandler, Excitations are
localized and relaxation is hierarchical in glass-forming liquids, Phys. Rev. X 1, 021013(2011), doi:10.1103/PhysRevX.1.021013.
[66] W. Zhang, J. F. Douglas, and F. W. Starr, What does the instantaneous normal mode spectrum
tell us about dynamical heterogeneity in glass-forming fluids?, J. Chem. Phys. 151, 184904(2019), doi:10.1063/1.5127821.

24