Low-frequency excitations and their localization properties in glasses

Besides the dynamical slowing down signaled by an enormous increase of the viscosity approaching the glass transition, structural glasses show interesting anomalous thermodynamic features at low temperatures that hint at peculiar deviations from Debye's law at low enough frequencies. Theory, numerical simulations, and experiments suggest that deviation from Debye's law is due to soft-localized glassy modes that populate the low-frequency spectrum. We study the localization properties of the low-frequency modes in a three-dimensional supercooled liquid model. The density of states $D(\omega)$ is computed considering the inherent structures of configurations well thermalized at parental temperatures close to the dynamical transition $T_\text d$. We observe a crossover in the probability distribution of the inverse of the participation ratio that happens approaching $T_\text d$ from high temperatures. We show that a similar crossover is observed at high parental temperature when the translational invariance of the system is explicitly broken by a random pinning field.


Introduction
Although glasses and amorphous materials are widespread in nature since the ancient times [1,2], a unified and coherent theoretical framework for describing their thermodynamical and dynamical properties remains still a challenge that attracts the attention of a wide scientific community [3,4]. Glasses can be obtained by cooling fast enough a liquid in order to avoid crystallization. In experiment and numerical simulations, the glass transition temperature is defined as the temperature at which the structural relaxation time τ α overcomes some threshold value. A glass can thus be seen as a fluid that does not flow anymore [5]. Under this perspective, static observables that are usually suitable for revealing positional order as the radial distribution function g(r) and its Fourier transform, i.e., the static structure factor S(q), do not indicate remarkable differences between the liquid and the glassy state.
Looking at glasses with the lens of solid state physics, it turns to be natural to study their lowenergy excitations. In particular, on large scales, glasses are continuum media and thus, at small enough frequencies, the density of states D(ω) follows Debye's law [6], i.e., D(ω) ∼ ω d−1 in d spatial dimensions. Debye's law assumes that the only low energy excitations are phonons. Debye's law provides precise theoretical predictions for the thermodynamic quantities such as the specific heat at low temperatures. However, differently from crystalline solids, glasses show anomalies in thermal conductivity and specific heat as temperature decreases towards zero [7]. For instance, the specific heat C v scales with T deviating from Debye's law. Moreover, thermodynamic anomalies are shared by different glassy systems suggesting that universal mechanisms are responsible for that [7].
Theoretical models as the two-level system model or the soft potential model face the problem looking at other excitations mechanisms besides phonons that should be taken into account for correctly describing the low excitations in disordered media. In particular, Gurarie and Chalker in reference [8] pointed out that non-Goldstone, and thus non-phononic, excitations in disordered systems contribute to D(ω) with a low-frequency sector that is universal, i.e., independent of the spatial dimensions, with a scaling D(ω) ∼ ω 4 . Moreover, such glassy modes are spatially localized and not extended like phonons. Since in five or less spatial dimensions, the predicted non-Goldstone contribution is subdominant with respect to Debye spectrum, it is hard to detect in both numerical simulations and experiments. Recently, a few numerical strategies have been developed in numerical simulations for taking access to the non-Goldstone contribution [9][10][11][12][13].
A simple strategy for probing the non-Goldstone sector of the spectrum consists of removing lowfrequency phonons. This can be done introducing an external field that breaks the translational invariance of the system [9]. Random pinning has been widely adopted in numerical simulations, analytical computations, and experiments to gain an insight into glassy transition, as a strategy for reaching the Kauzmann temperature and for measuring static correlations lengths [14][15][16][17][18][19][20][21]. In a previous work, we showed that random pinning can be employed for probing the non-Debye spectrum [12].
In reference [12] we showed that the low-frequency spectrum in a three-dimensional model of glass can be written as D(ω) ∼ ω s(p) , with p being the fraction of frozen particles. The exponent s(p) turns to be bounded by two extreme values, i.e., s(p) = 2 for p → 0 and s(p) = 4 above a threshold value p th that is of the order of 50 % of frozen particles. Such a phenomenology has a simple interpretation: as the number of frozen particles increases, phonons are pushed at higher frequencies and the moving particles rattle in a random environment. In particular, their equilibrium positions result to be randomly displaced with respect a crystalline configuration and thus these vibrations naturally give rise to a Rayleigh-type scattering mechanism [7].
Frozen particles provide an artificial tool for introducing heterogeneous regions with different elastic properties. We showed that a remarkably similar phenomenology emerges approaching the dynamical transition [22]. In particular, as it has been observed in reference [23], the low-frequency spectrum of D(ω) depends on the parental temperature T. We observed that one can write D(ω) ∼ ω s(T ) with s(T) = 2 for parental temperatures T T d , with T d the dynamical temperature, i.e., the temperature where the system undergoes the dynamical arrest. As T → T d we observed an increase in s(T), i.e., s(T) → 4 for T → T d . This happens because of the dynamical heterogeneities [24] that proliferate as temperature decreases towards the dynamical temperature. Since the behavior of s(T) mirrors that of s(p), it has been shown that the growing of spatially heterogeneous regions can be measured comparing the two systems. In this way, one can extract a dynamical correlation length ξ pin ∝ p as a function of the parental temperature T. ξ pin (T) shows a mild divergence at T d which is in agreement with the behavior of the dynamical correlation length ξ dyn computed through multi-point correlation functions.
In this paper, we investigate the localization properties of the low-frequency modes of a threedimensional glass former. Measuring the degree of localization of a mode of frequency ω through the inverse of its participation ratio R(ω), we find that R(ω) depends on the parental temperature T for modes populating the low-frequency spectrum. The low-frequency spectrum is defined as the frequencies below the boson peak that contribute to D(ω). In particular, approaching the dynamical transition, the probability distribution function of R(ω) shifts towards higher values. Defining R a as the first moment of the distribution, it turns out that R a undergoes a smooth crossover mirroring that in D(ω) ∼ ω s(T ) . We then compare the localization properties of low-frequency modes obtained considering a fraction p of particles frozen during the minimization of the mechanical energy. We obtain that also in the case of the pinned system, R a starts growing as p increases. We then compare the two protocols showing that it is possible to extract the behavior of a typical length scale ξ 3 ∝ p through the relation R a (T, p = 0) = R a (T = ∞, p). The mapping confirms that the properties of the pinned system at p → p th , with p th ∼ 0.5 provide complementary information on the same system at p = 0 and T → T d . In particular, the dependency p = p(T) is in agreement with other estimates [22].

Model
We consider a three-dimensional system composed of a 50 : 50 binary mixture of N soft spheres confined in a cubic box of side L with periodic boundary conditions and interacting through a pure repulsive pairwise potential [25,26]. We label large particles with A and the small ones with B. The total number of particles reads N = N A + N B and the corresponding density is ρ = N/L 3 . The radii are σ A and σ B with σ A /σ B = 1.2 and σ A + σ B ≡ σ = 1 [26]. The side of the box is L == N 1/3 such that ρ = 1. Indicating with r i the position of the particle i, with i = 1, . . . , N, two particles i, j interact via the potential φ( We impose a cutoff to the potential at r c = √ 3σ in a way that φ(r) = 0 for r > r c . The coefficients k 0 and k 2 ensure a continuity to φ(r) up to the first derivative at r = r c .

Equilibrium dynamics
For the dynamics, we have considered hybrid Brownian/Swap Monte Carlo simulations obtained combining the numerical integration of the equations of motion with Swap Monte Carlo moves [26]. In particular, in order to generate thermalized configurations, we propose an update of the system according to swap moves every 2 × 10 3 time steps. We consider system sizes N = 10 3 , 12 3 and averaging over 400 independent configurations. In what follows we report all quantities in reduced units considering σ = = µ = 1, where µ is the mobility of the Brownian particles. Figure 1 reports the behavior of the internal energy observable, t 0 is chosen such that the starting configuration is equilibrated at the temperature T, and t fin t 0 . Blue symbols refer to purely Brownian simulations, red symbols are hybrid Brownian/Swap simulations. As one can see, data obtained through Brownia/Swap simulations are well fitted by Rosenfeld and Tarazona (RT) formula [27] indicating that they are well thermalized. On the contrary, blue symbols deviate from RT meaning that the corresponding configurations did not reach thermal equilibrium. The dynamical temperature of the model T d has been computed fitting the structural relaxation time τ α with a power law τ α ∼ (T − T d ) −δ . τ α is defined as Q(τ α ) = e −1 . Q(t) is the self-overlap function between two configurations of the system, the first one taken at t = 0 and the second one at t [28]. Dynamical quantities have been computed considering the Brownian evolution of configurations that were previously thermalized through Brownian/Swap dynamics.

Inherent Structures and Density of States
After thermalization, we compute the corresponding inherent structures through the Limited-memory Broyden-Fletcher-Goldfarb-Shanno algorithm [29]. Let r be a configuration of the system, i.e., r ≡ (r 1 , . . . , r N ). The mechanical energy of the configuration r is E . The details about the minimization of E[r] with pinned particles can be found in reference [12]. We have computed all the 3N eigenvalues λ κ , with κ = 1, . . . , 3N, using gsl-GNU libraries for sizes up to N = 12 3 . The corresponding eigenvalues of M are ω 2 κ = λ κ . To identify the low-frequency spectrum, we focus our attention on the cumulative F(ω)= is the density of states. N is the number of non-zero modes that is 3N − 3 for translational invariant systems.
The localization properties of the normal-modes have been investigated through the inverse participation ratio R(ω) defined as R(ω) ≡ i |e i (ω)| 4 / i |e i (ω)| 2 2 where e i (ω) is the eigenvector of the mode ω [30]. For a mode ω completely localized on a single particle, one has R(ω) = 1, while a mode extended over all the particles corresponds to R(ω) ∼ N −1 . We also compute P(R) ≡ Q ω λ :ω λ <ω th δ [R(ω λ ) − R] that is the probability distribution of the inverse participation ratio for the modes with frequency ω smaller than a threshold frequency ω th . Q is a normalization constant. After computing the distribution P(R), we measure R a ≡ log R .

Results
As it has been shown in reference [12], the non-Goldstone sector becomes clearly visible in D(ω), i.e., its weight in the density of states overcomes phonons, in systems with a finite fraction p of frozen particles. In particular, non-Goldstone modes result to be soft, i.e., D(ω) ∼ ω s(p) with s(p) > 2 and localized, i.e., R(ω) does not scale as 1/N. These facts are in agreement with the soft potential model that predicts the scaling D(ω) ∼ ω 4 for non-Goldstone excitations around the absolute minima in a one-dimensional random energy landscape [8]. Moreover, D(ω) at low frequencies changes its features  when it is computed considering configurations thermalized with different protocols [23]. Thermalizing the system closer and closer to the dynamical transition, the corresponding inherent structures show the same type of crossover that can be described through a scaling D(ω) ∼ ω s(T ) [22]. Moreover, the same crossover can be documented through different observables, i.e., the effective exponent s(T), the distribution of displacements travelled by particles for reaching the inherent structures, the mean-distance travelled by a particle for reaching the optimal configuration. Here, we focus our attention on the R of the low-frequency modes and their probability distribution function P(R). We thus consider the inherent structures obtained starting from configurations thermalized at parental temperature T that approaches the dynamical temperature T d .

Localization of the glassy modes
In figure 3 we show the distribution P(R) at temperature T/T d = 1.42 (blue) and T/T d = 1.03 (red). The distribution that has been computed considers only low-frequency modes. To do that, we impose a cutoff frequency ω th that is chosen below the boson peak in the region where the power law scaling D(ω) ∼ ω s(T ) holds [22]. As one can appreciate, the distribution changes the shape and shifts towards higher R values as temperature decreases indicating that extended modes have been progressively suppressed.
To gain an insight into the role played by the cutoff frequency ω th , we have computed P(R) for different values of ω th . The cutoff ω th is taken below the frequency of the boson peak ω BP which, in our models, is around ω BP ∼ 0.1. The location of the boson peak in three spatial dimensions can be obtained looking at the maximum of D(ω)/ω 2 as shown, for instance, in reference [31]. We thus choose ω th ∈ [0.02, 0.1]. The results are shown in figure 4 (a) in the case of T/T d = 1.16. The presence of extended modes at frequencies larger than 0.04 dramatically changes the shape of P(R). This is made evident when one looks at R a as a function of temperature as it is shown in figure 4 (b). Choosing ω th in the low-frequency region, i.e., ω th < 0.08, R a undergoes a smooth crossover as temperature decreases towards T d . In particular, data for ω th = 0.02, 0.03 (blue symbols), i.e., a frequency that is below the boson peak, show that R a starts to increase for T/T d − 1 < 0.5.
This finding is confirmed in figure 5 where we plot R a as a function of ω th for different values of T/T d = 1.03, 1.29, 1.93 (blue, green, and red symbols, respectively). We observe a region that is independent of both, ω th and parental temperature T. Below ω th ∼ 0.05, R a increases as ω th decreases but also as T decreases. This is a clear signal that in the low-frequency region extended modes become attenuated in configurations thermalized at lower parental temperatures. It is worth noting that the behavior of R a as a function of T mirrors the behavior of s(T) as shown in figure6 (b) where green symbols are s(T) + 1 obtained from the cumulative distribution F(ω) [22].

Comparison between (T , p = 0) and (T = ∞, p) protocol
Now we are going to compare the localization properties of the eigemodes obtained at high parental temperatures but considering a finite fraction of a frozen particle during the minimization of the mechanical energy with those obtained in the non-pinned system as a function of T. It has been shown that, as p increases, D(ω) ∼ ω s(p) with s(p) undergoing a smooth crossover from s = 2 to s = 4 above a threshold value p th ∼ 0.5 [12]. We refer to this protocol as (T = ∞, p), since one considers configuration thermalized at high temperatures, i.e., away from T d . In the previous section, we have considered a protocol where the inherent structures were computed considering configutations thermalized at parental temperatures T close to T d . We refer to this protocol as (T, p = 0) since energy minimization has been computed without any artificially frozen particle. Also in the case of the pinned system, modes responsible for ω 4 are localized. This is because they live in between the frozen regions induced by the random pinning protocol. Since the number of frozen particle increases, the lowest frequency of the spectrum naturally shifts towards higher values and the threshold value ω th as well. Random pinning explicitly destroys the translational invariance removing the corresponding three zero modes from the spectrum and attenuating any extended mode. This has a strong effect on R a that grows towards R a ∼ 0 for p → 1, as it is shown in figure 6 (a), blue symbols. Also in this case, similarly to what we observed in the previous protocol, i.e., decreasing the parental temperature without any artificially frozen particle, R a undergoes a crossover mirroring that in s(p) (black symbols in the same panel). It is worth noting that in both protocols when the density of states undergoes a crossover s → 4, R a ∼ −1.5.
Since frozen particles occupy a finite volume fraction pV, we can associate to the volume fraction a typical length scale ξ≡ (pV) 1/3 . In reference [22] it has been shown that, looking at the solution of s(p, T = ∞) = s(T, p = 0), one can extract the behavior of ξ 3 as a function of the parental temperature T. Here, we can extract ξ 3 looking at the degree of localization of the low-frequency modes. To do that, we solve and invert numerically R a (T, The result is shown in figure 7, red symbols. Green symbols refer to the solution of s(p, T = ∞) = s(T, p = 0). As one can see, the two data sets are in a good agreement.

Conclusions
In this paper, we have investigated the localization properties of low-frequency modes in a threedimensional model of supercooled liquid. In particular, we have focused our attention on the role played by the parental temperature on the localization of the soft glassy modes. Our findings show that low- frequency vibrational modes at lower parental temperature turn out to be more localized than those populating the density of states at higher T values. This is consistent with the results presented in [22,23] and also with simulations on well-equilibrated polydisperse glass former [13]. In particular, the increasing in localization takes place near the dynamical transition temperature T d that is where the exponent s(T) of the power law D(ω) ∼ ω s(T ) approaches s(T d ) → 4. This finding confirms the interplay between dynamical and zero-temperature structural properties in glasses [22]. At lower temperatures, we also observed a dependency of R a on the cutoff frequency ω th . This shift reminds the effect of random pinning on the density of states [12].
We have thus investigated how the localization of the lowest eigenmodes takes place in the same system with random pinning. In particular, we observed that the same scenario of progressive localization of glassy modes takes place as the number of frozen particles increases. In the pinning protocol, the emergence of soft localized excitations is due to the breaking of translational invariance in the system. With increasing p, moving particles rattle into small islands that are surrounded by the frozen ones. At higher values, i.e., for p > p th , the phononic spectrum is totally destroyed giving rise to extremely localized modes, i.e., R → 1. This marks a difference with a system thermalized at lower parental temperatures where the translational invariance is preserved. Nevertheless, configurations thermalized at lower temperatures show a spectrum of harmonic vibrations whose properties are remarkably similar to those obtained breaking explicitly the spatial translational invariance, i. e., crossover in D(ω) from Debye to non-Debye, localization of low-frequency modes, caging effects during minimization. We can thus extract useful and complementary information comparing the region p < p th , in the (T = ∞, p) protocol, with T > T d , in the (T, p = 0) protocol. In this paper, we provide evidences for a crossover from extended to localized modes at low-frequencies as temperature decreases towards T d which is in agreement with recent works [13,32]. We also showed that, in analogy with reference [22], R a can be employed for mapping structural properties into dynamical ones. In particular, the degree of localization measured through R a is regulated by the proliferation of dynamical heterogeneous regions in the (T, p = 0) protocol. The (T = ∞, p) protocol allows one to define a length scale ξ ∼ p 1/3 that is an external and tunable parameter. We can thus study p = p(T) just looking at the solution of O(T, p = 0) = O(T = ∞, p), with O being a generic observable. As it has been shown in reference [22], p(T) mirrors the behavior of the dynamical length scale ξ dyn [33]. Here, we showed that from R a we can extract the behavior of p(T) which is in agreement with those observed through other observables [22].
As a future direction, it would be interesting to study the density of states in systems thermalized at parental temperature T close to T d with a fraction p of pinned particles. In this way, in the spirit of early works on random pinning [14][15][16][17][18][19][20][21], it would be possible to take access to the properties of glassy modes towards the Kauzmann temperature.