Spatial structure of unstable normal modes in a glass-forming liquid

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. Through this 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 the mode-coupling temperature the unstable modes are localized around structural defects, characterized by a disordered local structure markedly different from the liquid's locally favored 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.

The phenomenology of glass-forming liquids is often described in terms of their underlying, highdimensional 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. Through this 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 the mode-coupling temperature the unstable modes are localized around structural defects, characterized by a disordered local structure markedly different from the liquid's locally favored 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.

I. INTRODUCTION
Glass formation is driven by the rapid increase of structural relaxation times that takes place when 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 T g , at which the structural relaxation time reaches about 100s. Below T g , the supercooled liquid behaves, on typical experimental time scales, as an amorphous solid. However, the lack of a sharp divergence in dynamic 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 underlying mechanisms responsible for the enormous slowing down of the dynamics are still highly debated [3][4][5].
A long-standing line of thought describes glass formation in terms of the exploration of a high-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 to either 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 T x , at which the structural relaxation time is of order 10 −9 s. In this scenario, T x marks the crossover between the normal liquid * masanari-shimada444@g.ecc.u-tokyo.ac.jp dynamics and activated 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 down of 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 T MCT , at which the theory predicts a power law divergence of the relaxation times. This singularity is expected to be smeared out in actual supercooled liquids, because of the presence of thermally activated processes not accounted for by the theory, which turn 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 T MCT , the dynamics is governed by the motion along unstable saddle modes, which the system can exploit without resorting to thermal activation; below T MCT , instead, the system mostly vibrates around local minima of the PES, occasionally relaxing via activated barrier crossing, and one can thus identify T MCT ≈ T x . This theoretical scenario is largely inspired by the behavior of a class of mean-field spin models, 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][13][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][16][17][18][19][20][21][22][23][24][25]. These difficulties were due to two main factors. First, conventional simulations cannot easily sample the PES at equilibrium below T MCT , due to exceedingly long equilibration times. Second, all the above mentioned studies neglected the role of the modes' spatial localization. In particular, it is clear that thermally activated processes in finite dimensional liquids involve only localized portions of the system, whereas the unstable modes of saddles in meanfield models are spatially delocalized. A recent numerical study [26], based on an efficient swap Monte Carlo algorithm [27,28], has tackled these two issues at once, establishing that the stationary points of several threedimensional model liquids indeed show a sharp change around T MCT : 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 concentration of 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 is natural to inquire the connection between them and the stable modes populating the low-frequency 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. It is now established that the vibrational density of states of the QLVs shows a quartic dependence on the frequency [29][30][31][32][33][34][35][36][37], while they are suppressed when pre-stress is removed [30,38]. As for their spatial structure, QLVs consist of an energetically unstable core and a stable far-field [39] 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 a distance from a core and d is the spatial dimension, unless hybridization with acoustic modes occurs [29,33,40]. 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 sheartransformation zones [41][42][43][44][45][46].
In this work, we apply the spatially resolved mode analysis developed to study QLVs [29,39] 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 searching for stationary points in moderately large system sizes. Through this analysis, we confirm that the delocalized unstable modes that populate the saddles' spectrum at high temperature have an 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 to that 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 that around and below T MCT 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 [47]. 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. II A, we introduce our model and describe how we located saddles and local minima starting from equilibrium configurations. In Sec. II B and II C, we define the quantities used to characterize these stationary points. We provide our main results in Sec. III. We focus on the spatial structure of unstable saddle modes in Sec. III A and the local structure of unstable cores in Sec. III C. We compare saddles and local minima in Sec. III B. Finally, we conclude our work with a summary in Sec. IV.

A. Sample preparation
We present numerical results for the ternary mixture model introduced by Gutirrez et al. [48]. The model consists of three species of particles, which interact via a repulsive inverse power potential with exponent 12 where α, β = A, B, C are species indices. The additional terms in Eq.(1) ensure continuity of the potential up to its second derivative at the cutoff distance r cut = 1.25σ αβ . See Ref. [49] for the expressions of c 0 , c 2 , and c 4 . The interactions are additive, the size ratio is σ AA σ BB = σ BB σ CC = 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. Results for N = 1000 are shown in Appendix C. The estimated mode-coupling critical temperature is T MCT ≈ 0.29 [50]. Starting from these equilibrium configurations, we performed three kinds of optimizations: (i) potential energy minimizations using the l-BFGS method to locate the local minima of the PES; (ii) total square-force minimizations (or W -minimizations) using the l-BFGS method [51] to locate local minima of where F i is the force acting on particle i; (iii) eigenvectorfollowing (EF) optimizations [52] to locate stationary points of the PES with a prescribed number n u of unstable modes.
Our goal is to locate stationary points, either local minima or saddles, in the neighborhood of 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 order n u , but the value of n u has to be provided as input. To define a mapping between a given configuration and its neighboring stationary point we have followed the approach used in Ref. [26]: for each starting configuration, we first perform a W -minimization. The order n u of the W -minimized configuration is then used as the target order for an EF optimization. All the saddles analyzed in Sec. III 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 system size of interest in this work (N = 3000). To fix this issue, in this work we have 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 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. 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 was larger than a threshold (10 4 ), and decreased all the trust radii until the potential energy of the new configuration dropped below the threshold. Except for these small but important practical details, the algorithm was the same as the one used in Ref. [26].

B. Normal mode analysis
We performed a standard normal mode analysis [53] for local minima and saddles. We diagonalized the dynam-ical matrix, which is the Hessian of the total potential energy U , where is the contribution from the pair ij , u ij (r ij ) is the pair interaction potential, n ij = r ij /r ij is the unit vector along the pair, and I is the d × d identity matrix. We denote its eigenvalues by λ α 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 three eigenmodes corresponding to global translations. Using these eigenmodes we computed the quantities defined below.

Participation ratio
The participation ratio [54][55][56] of a mode α is given by 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.

Decay profile
The decay profile [29,40] d α (r) of a mode α is a function of the distance r from the particle with the largest | e α,i |, which we denote by i d * . The decay profile is defined as where r i is the distance of a particle i from the particle i d . In the case of saddles, we show the averaged decay profile for a certain eigenvalue λ defined as * The distance r used here differs slightly from the one used in the paper by Gartner and Lerner [40]. They defined the center of an eigenmode using w particles with the largest | e α,i | 2 and measured the distance r from this center. They mainly used w = 4 while our definition corresponds to w = 1. This difference does not affect the our results qualitatively.

Fractal dimension
To further characterize the spatial structure of the eigenmodes, we introduce another function of the norms | e α,i |. We identify the particles contributing to a mode α as those with the P (λ α ) largest norms and then determine the number N α (r) of such contributing particles up to a distance r from the particle i d . We define the averaged function N (r) by the same procedure as in Eq. (7). From this function, we can estimate the fractal dimension D of mode of eigenvalue λ from N (r) ∼ r D .

Energy profile
We define the local vibrational energy of a particle i on mode α as where e α,ij = e α,i − e α,j is the relative displacement between the pair and ( e ⊥ α,ij ) 2 = ( e α,ij ) 2 − ( n ij · e α,ij ) 2 is the transverse relative displacement. The energy profile [39] Λ α (r) of a mode α is a function of the distance r from the particle with the smallest δE α,i , which we denote by i e . Λ α (r) is defined as where r i is the distance of a particle i from the particle i e . In the following, we show the averaged normalized energy profile defined as for saddles, while we used only the lowest frequency QLVs for local minima. Since i δE α,i = λ α ,Λ(r) → −1 as r → ∞ when λ < 0 whileΛ(r) → 1 when λ > 0. Finally, we note that the difference between i d and i e does not matter in practice, because there is a negative correlation between | e α,i | and δE α,i (see Appendix A).

One-particle dynamical matrix
Finally, using M ij in Eq. (4), we define the one-particle dynamical matrix as Clearly, the full dynamical matrix reduces to H ij when only a tagged particle is allowed to move and the others are kept fixed. We denote its eigenvalues by µ β and the eigenvectors by f β , where β = 1, 2, . . . , N d.

C. Locally favored structure
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 [57]. The radical Voronoi tessellation requires as input the typical diameters r α of particles of each type α, which we estimated from the positions of the first peaks of the partial radial distribution functions g αα (r) measured for minimized configurations at T = 0.28, namely r 1 = 1.34, r 2 = 1.11, r 3 = 0.93. Following a well-established approach [58], we characterized the shape the Voronoi polyhedra using the signature (n 3 , n 4 , n 5 , . . . ), where n k is the number of faces of the polyhedron with k faces. We found that the (0,0,12) signature, corresponding to an icosahedral local arrangement, becomes the most frequent at low temperature, in either minimized or instantaneous configurations. Following Ref. [24], we therefore identify the icosahedron as the locally favored structure of the model.

III. RESULTS
A. Spatial structure of unstable modes 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 the participation ratio in Fig. 1 In Ref. [26], it was established that the average participation ratio of the modes satisfies P (λ < λ e )N −1/3 → 0 as N → ∞ below the mobility edge λ e , while P (λ > λ e )N −1/3 → ∞ above it. 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] are shown as dotted lines in Fig. 1; note that λ e vanishes around T MCT ≈ 0.29. Table I shows the specific values of λ and λ e we used to perform the average in Eq. (7). We used a bin 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 and 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 P (λ α ) particles having the largest norms, where [. . . ] denotes the integer part. Different colors indicate different particle types. The particle with the largest norm is placed at the center of the box and the displacements are shown by black arrows.    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 the cores of the localized modes modes in Fig. 2(a) and (b), while it is difficult to identify similar cores in the 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 this question, we have analyzed the network of the [P ] particles which, for a given mode, have the largest 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 to a future investigation.
Let us start our analysis by looking at the decay profiles. We remind that the decay profiles of 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 their eigenvalues as described in Sec. II B 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. I: 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 our data for the QLVs, we cannot observe a clear asymptotic behavior d(r) ∝ r −2 because our 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 that the 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 achieve this, 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 This figure shows that the modes with λ ∼ λ e (shown in gray) follow N (r) ∼ r 1 . This is consistent with the definition of the mobility edge, at which P grows linearly with L. We note, however, that N (r) is an averaged quantity and there are strong mode-to-mode fluctuations (not shown). Thus, we should 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 of the localized modes rapidly converge to their final values while the delocalized modes and the QLVs grow faster than linear with r. Nevertheless, we can also see a difference between the delocalized unstable modes and the QLVs, particularly at T ≥ 0.35: the latter saturate at large distances. This means that the particles contributing to the QLVs are denser and more compact than those contributing to the delocalized modes. This observation 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 vertical lines are placed at half of the box length L/2. We find that localized and delocalized unstable modes display qualitatively different energy profiles. The former have pronounced dips at r ∼ 2, which are the cores of the modes [32], and for r 2 they rapidly converge to −1 from below. Thus, the localized modes have weak farfield 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 unstable core 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, the delocalized modes have the large farfield components and display a pronounced temperature dependence. The energy profiles of the delocalized modes also show a small dip (or kink) at r ∼ 2, which indicates that even these modes may possess a core. However, at T ≥ 0.32 we observeΛ(r ∼ 2) > −1. Therefore, the energy profiles converge to −1 from above. In other words, both the core and the far field are unstable. This is significantly different from the localized modes and also from the QLVs.
For T ≤ 0.30, i.e. below the localization transition, the energy profiles converge to −1 from below and the functional form is similar to that of the QLVs. This indicates that the softest localized unstable modes that survive below the mode-coupling crossover [26] have spatial properties similar to the soft stable modes, see also Sec. III C. We emphasize, however, that the energy profiles discussed here characterize only the average spatial features of the vibrational modes. Dips or kinks at r ∼ 2 in the energy profiles indicate the presence of regions that are particularly unstable on average. Our analysis does not tell, however, how many 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 the mode is strongly localized, as observed in Fig. 2(a) and in Fig. 6(a), the particle i e 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" [41], although the two definitions differ in practice.
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 in Fig. 2 and the particles with the most negative δE α,i are placed at the center of the squares. Particles having more negative energy are colored with darker colors, while those with positive energy are colored white, regardless of the value of δE α,i . We can see a few dark particles at the centers of the boxes in all cases; they correspond to the small dips at r ∼ 2 observed in Fig. 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 unstable backgrounds, i.e., there are scattered yellow regions. This is consistent with the differences between localized and delocalized modes mentioned above.   8. PDF of the eigenvalues of the one-particle dynamical matrix. All parameters are the same as in Fig. 7.

B. 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 expect that this reflects a structural difference between saddles and local minima. Otherwise, it may 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 to each 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 [39]. 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 and saddles using a Voronoi tessellation, see Sec. II C. 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 Voronoi signatures in saddles and minima have markedly 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 in minima and saddles becomes very similar sampled at low temperature (T = 0.28), see Tab. II, 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 that these structural differences fade away close to T MCT . The strong structural similarity between saddles and minima at low temperature suggests that the localized unstable modes find their origin in very subtle features of the fluid structure.

C. 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 [41,47]. To address this question, 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 analysing 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. II B 5 for details). Here we show that its eigenmodes give us an insight into the cores of the localized unstable modes of 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 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 oneparticle 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 as where 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 with smaller P (λ α ) have larger overlap with the eigenmodes of H (1) . Therefore, the motions of the strongly 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-body correlation functions, as follows: we restrict the calculation of the radial distribution function g(r) to central particles of species 3 that form the cores of the unstable modes, i.e., the particles whose index is i e † . Note that any other neighbors at a distance r is used in this calculation, irrespective of its species. Similarly, we compute the restricted bond-angle distribution b(θ) obtained from the angles formed between a central particle of species 3 belonging to the cores of the unstable modes and two of its nearest neighbors, irrespective or 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 the local structure around core particles themselves.
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, the g(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 core particles are unlikely to align linearly with their neighbors [62], which is consistent with the absence of a peak at 180 degrees in the bond-angle distribution. Overall, the structure around core particles is almost featureless and resembles the one of the fluid at higher temperature.
Finally, we also analyzed the signatures of the Voronoi cells surrounding the core particles and found 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 II  and III. 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 structural features: 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 unstable modes and QLVs differ. Overall, the results indicate that, around and below the MCT crossover, localized un-stable modes and soft stable modes have a similar structural origin.
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 in granular packings [41] and in a model metallic glass [47]. We note, however, that while the locally favored structure is the same in our model and in the one of Ref. [47], the Voronoi signatures of the structural defects have a very broad, model-specific distribution. Thus, this kind of Voronoi-based analysis does not clearly tell us what structural defects are -at best what they are not. Devising a 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 [61,64,65].

IV. SUMMARY AND 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 kinds of unstable modes: (i) delocalized modes, which are characterized by spatially extended displacement fields and which disappear below the mode-coupling crossover temperature T MCT and (ii) localized modes, whose properties are largely independent of temperature. We have also performed a detailed comparison to the features of the softest stable modes observed in local 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 system sizes 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 meanfield glass models. By contrast, localized unstable modes find their origin in structural defects of the liquid, and are 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 useful framework 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 [66]. See also Ref. [67] 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 → T MCT . Our independent analysis confirms the presence of a qualitative change of the spatial localization features across the mobility 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 a very 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-mode level.
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 the local structure of local minima and saddles differ significantly, and the decay profiles and fractal dimensions of QLVs do not clearly match the ones of either delocalized or localized unstable modes. However, as the temperature drops around and below T MCT , 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 and 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. The origin of these subtle differences is probably encoded in a complex way in the energy function and cannot be fully revealed yet by our methodology. Although the computational costs 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 the saddle 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 and energy 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 localized around similar structural defects. Our study connects previous investigations of the PES, which have 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. The origins of the functions d(r) and Λ(r) are different. However, it was shown in Ref. [39] 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 11 shows the data at T = 0.32. Different colors in (a) indicate different eigenvalues shown in Tab. I. To average | e α,i | and δ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. To compare with the results of saddles, we show the results of local minima. The decay profile, N (r)/N , and the energy profile are shown in Fig. 12(a), (b), and (c), respectively. Different symbols and colors indicate different temperatures. Unlike the case of saddles, we averaged them over the lowest-frequency QLVs in each sample. We cannot see any strong dependence on the temperature although the modes are slightly more localized when the temperature decreases.  To study the size dependence, we show the data of N = 1000. Figure 13 shows the corresponding decay profile. We show the data of N = 1000 by lines and those of N = 3000 by symbols. The specific eigenvalues used to compute the average is shown in Tab. IV. We can hardly see any size dependence. In contrast, we can see the size dependence of the energy profile in Fig. 14. In particular, the delocalized unstable modes show a strong size dependence at high temperatures while the localized modes do not depend on N . This is natural because the 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 depend on N in any case.