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 also interesting anomalous thermodynamic features at low temperatures that hint 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 is computed considering the inherent structures of configurations well thermalized at parental temperatures close to the dynamical transition . We observe a crossover in the probability distribution of the inverse of the participation ratio that happens approaching 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.
Although glasses and amorphous materials are widespread in nature since the ancient time Binder and Kob (2011); Leuzzi and Nieuwenhuizen (2007), 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 Berthier and Biroli (2011); Debenedetti and Stillinger (2001). Glasses can be obtained 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 see as a fluid the does not flow anymore Cavagna (2009). Under this perspective, static observables that are usually suitable for revealing positional order as the radial distribution function and its Fourier transform, i. e., the static structure factor , 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 low-energy excitations. In particular, on large scales, glasses are continuum media and thus, at small enough frequencies, the density of states follows Debye’s law Kittel (2005), i. e., in spatial dimensions. Deby’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 Zeller and Pohl (1971). For instance, the specific heat scales with deviating from Debye’s law. Moreover, thermodynamic anomalies are shared by different glassy systems suggesting that universal mechanisms are responsible for that Zeller and Pohl (1971).
Theoretical models as the two-level system model or the soft potential model face the problem looking at other excitations mechanisms besides phonons that have to be taken into account for describing correctly the low excitations in disordered media. In particular, Gurarie and Chalker in Ref. Gurarie and Chalker (2003) pointed out that non-Goldstone, and thus non-phononic, excitations in disordered systems contribute to with a low-frequency sector that is universal, i. e., independent of the spatial dimensions, with a scaling . 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 Baity-Jesi et al. (2015); Lerner et al. (2016); Mizuno et al. (2017); Angelani et al. (2018); Wang et al. (2019).
A simple strategy for probing the non-Goldstone sector of the spectrum consists of removing low-frequency phonons. This can be done introducing an external field that breaks the translational invariance of the system Baity-Jesi et al. (2015). Random pinning has been widely adopted in numerical simulations, analytical computations, and experiments for gaining insight into glassy transition, as a strategy for reaching the Kauzmann temperature and for measuring static correlations lengths Cammarota and Biroli (2013, 2012); Szamel and Flenner (2013); Kob et al. (2012); Gokhale et al. (2014); Karmakar and Parisi (2013); Kob and Berthier (2013); Ozawa et al. (2015). In a previous work, we showed that random pinning can be employed for probing the non-Debye spectrum Angelani et al. (2018).
In Ref. Angelani et al. (2018) we showed that the low-frequency spectrum in a three-dimensional model of glass can be written as , with the fraction of frozen particles. The exponent turns to be bounded by two extreme values, i. e., for and above a threshold value that is of the order of 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 Zeller and Pohl (1971).
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 Paoluzzi et al. (2019). In particular, as it has been observed in Ref. Lerner and Bouchbinder (2017), the low-frequency spectrum of depends on the parental temperature . We then observed that one can write with for parental temperatures , with the dynamical temperature, i. e., the temperature where the system undergoes the dynamical arrest. As we observed an increase in , i. e., for . This happens because of the dynamical heterogeneities Kob et al. (1997) that proliferate as temperature decreases towards the dynamical temperature. Since the behavior of mirrors that of , 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 as a function of the parental temperature . shows a mild divergence at that is in agreement with the behavior of the dynamical correlation length computed through multi-point correlation functions.
In this paper we investigate the localization properties of the low-frequency modes of a three-dimensional glass former. Measuring the degree of localization of a mode of frequency through the inverse of the participation ratio , we find that depends on the parental temperature for modes populating the low-frequency spectrum. The low-frequency spectrum is defined as the frequencies below the Boson peak that contribute to . In particular, approaching the dynamical transition, the probability distribution function of shifts towards higher values. Defining as the first moment of the distribution, it turns that undergoes a smooth crossover that follows the crossover observed in . We then compare the localization properties of low-frequency modes obtained considering a fraction of particles frozen during the minimization of the mechanical energy. We obtain that also in the case of the pinned system starts growing as increases. We then compare the two protocol showing that it is possible to extract the behavior of a typical length scale through the relation . The mapping confirms that the properties of the pinned system at , with , provides complementary information on the same system at and . In particular, the dependency is in agreement with other estimates Paoluzzi et al. (2019).
We consider a three-dimensional system composed of a binary mixture of soft spheres confined in a cubic box of side whit periodic boundary conditions and interacting through a pure repulsive pairwise potential Bernu et al. (1987); Grigera and Parisi (2001). We label large particles with and the small ones with . The total number of particles reads and the corresponding density is . The radii are and with and Grigera and Parisi (2001). The side of the box is such that . Indicating with the position of the particle , with , two particles interact via the potential where . We impose a cutoff to the potential at in a way that for . The coefficient and guaranty continuity to up to the first derivative at .
ii.1 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 Montecarlo moves Grigera and Parisi (2001). In particular, in order to generate thermalized configurations, we propose an update of the system according to swap moves every time steps. We consider system sizes and averaging over independent configurations. In the following we report all quantities in reduced units considering , where is the mobility of the Brownian particles.
Fig. (1) reports the behavior of the internal energy , where the angular brackets indicate averages over trajectories, i.e., , with a generic observable, is chosen such that the starting configuration is equilibrated at the temperature , and . Blue symbols refers 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 Rosenfeld and Tarazona (1998) 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 has been computed fitting the structural relaxation time with a power law . is defined as . is the self-overlap function between two configurations of the system, the first one taken at and the second one at LaÄeviÄ et al. (2003). Dynamical quantities have been computed considering the Brownian evolution of configurations that were previously thermalized throw Brownian/Swap dynamics.
ii.2 Inherent Structures and Density of States
After thermalization, we compute the corresponding inherent structures through the Limited-memory Broyden-Fletcher-Goldfarb-Shanno algorithm Bonnans et al. (2006). Let be a configuration of the system, i. e., . The mechanical energy of the configuration is . We indicate with the configuration that minimizes , We define the Inherent Structure energy , with indicating the average over independent configurations. The temperature dependence of is shown in Fig. (2).
The spectrum of the harmonic oscillations around is then obtained considering a perturbed configuration . The mechanical energy now reads with with the elements Hessian matrix , where latin indices indicate the particles and greek symbols the cartesian coordinates. For estimating the correlation length , we also considered configurations where a finite number of particles , with the particle fraction, are maintained frozen during the minimization of . The details about the minimization of with pinned particles can be found in Ref. Angelani et al. (2018). We have computed all the eigenvalues , with , using gsl-GNU libraries for sizes up to . The corresponding eigenvalues of are .
To identify the low-frequency spectrum, we focus our attention on the cumulative =, where = is the density of states. is the number of non-zero modes that is for translational invariant systems.
The localization properties of the normal-modes have been investigated through the inverse participation ratio defined as where is the eigenvector of the mode Bell and Dean (1970). For a mode completely localised on a single particle, one has , while a mode extended over all the particles corresponds to . We also compute that is the probability distribution of the inverse participation ratio for the modes with frequency smaller than a threshold frequency . is a normalization constant. After computing the distribution , we measure .
As it has been shown in Ref. Angelani et al. (2018), the non-Goldstone sector becomes clearly visible in , i. e., their weight in the density of states overcomes phonons, in systems with a finite fraction of frozen particles. In particular, non-Goldstone modes result to be soft, i.e., with and localized, i .e., does not scale as . These fact are in agreement with the soft potential model that predicts the scaling for non-Goldstone excitations around the absolute minima in a one-dimensional random energy landscape Gurarie and Chalker (2003). Moreover, at low frequencies changes its features when it is computed considering configurations thermalized with different protocols Lerner and Bouchbinder (2017). The inherent structures obtained thermalizing the system closer and closer the dynamical transition show the same type of crossover that can be described through s scaling Paoluzzi et al. (2019). The same crossover can be documented through different observables, i. e., the effective exponent , 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 of the low-frequency modes and their probability distribution function . We thus consider the inherent structures obtained starting from configurations thermalized at parental temperature that approaches the dynamical temperature .
iii.1 Localization of the glassy modes
In Fig. (3) we show the distribution at temperature (green) and (yellow). The distribution has been computed consider only low-frequency modes. For doing that, we impose a cutoff frequency that is chosen below the Boson peak in the region where the power law scaling holds Paoluzzi et al. (2019). As one can appreciate, the distribution changes shape and shifts towards higher values as temperature decreases indicating that extended modes have been progressively suppressed.
To gain insight into the role played by the cutoff frequency , we have computed for different values of . The cutoff is taken below the frequency of the Boson peak that, in our models, is around . The location of the Boson peak in three spatial dimensions can be obtained looking at the maximum of as shown, for instance, in Ref. Grigera et al. (2003). We thus choose . The results are shown in Fig. (4a) in the case of . The presence of extended modes at frequencies larger than dramatically changes the shape of . This is made evident when one looks at as a function of temperature as it is shown in Fig. (4b). Choosing in the low-frequency region, i.e., , undergoes a smooth crossover as temperature decreases towards . In particular, data for (blue symbols), i.e., a frequency that is below the boson peak, show that starts to increase for .
This finding is confirmed in Fig. (5) where we plot as a function of for different values of (blue, green, and red symbols, respectively). We observe a region that is independent of both, and parental temperature . Below , increases as decreases but also as 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 as a function of mirrors that of as shown in Fig. (6b) where green symbols are obtained from the cumulative distribution Paoluzzi et al. (2019).
iii.2 Comparison between and protocol
Now we are going to compare the localization properties of the eigemodes obtained at high parental temperatures but considering a finite fraction of frozen particle during the minimization of the mechanical energy with those obtained in the non-pinned system as a function of . It has been shown that, as increases, with undergoing a smooth crossover from to above a threshold value Angelani et al. (2018). We refer to this protocol as , since one consider configuration thermalized at high temperatures, i. .e, away from . In the previous section, we have considered a protocol where the inherent structures were computed considering configutations thermalized at parental temperatures close to . We refer to this protocol as since energy minimization has been computed without any artificially frozen particle.
Also in in the case of pinned system the modes modes responsible for are localized. This is because they live in between the frozen regions induced by the random pinning protocol. Since as the number of frozen particle increases, the lowest frequency of the spectrum naturally shifts towards higher values and the threshold value as well. Random pinning explicitly destroys the translational invariance removing the corresponding three zero modes from the spectrum and attenuating any extended mode. This have a strong effect on that grows towards for , as it is shown in Fig. (6a), blue symbols. Also in this case, similarly whit what we observed in the previous protocol, i.e., decreasing the parental temperature without any artificially frozen particle, undergoes a crossover mirroring that in (black symbols in the same panel). It is worth noting that in both protocols when the density of states undergoes a crossover , .
Since frozen particles occupy a finite volume fraction , we can associate to the volume fraction a typical length scale . In Ref. Paoluzzi et al. (2019) it has been shown that, looking at the solution of , one can extract the behavior of as a function of the parental temperature . Here we can extract looking at the degree of localization of the low-frequency modes. For doing that, we solve and invert numerically the equation . In particular, we compared the growing of with The result is shown in Fig. (7), red symbols. Green symbols refer to the solution of . As one can see, the two data sets are in a good agreement.
In this paper, we have investigated the localization properties of the low-frequency modes in a three-dimensional 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 to be more localized than those populating the density of states at higher values. This is consistent with results presented in Paoluzzi et al. (2019); Lerner and Bouchbinder (2017) and also with simulations on well-equilibrated polydisperse glass forming Wang et al. (2019). In particular, the increasing in localization takes place near the dynamical transition temperature that is where the exponent of the power law approaches . This finding confirms the interplay between dynamical and zero-temperature structural properties in glasses Paoluzzi et al. (2019). At lower temperatures, we observed also a dependency of on the cutoff frequency . This shift reminds the effect of random pinning on the density of states Angelani et al. (2018).
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 emerging of soft localized excitations is due to the breaking of translational invariance in the system. With increasing , moving particles rattle into small islands that are surrounded by the frozen ones. At higher values, i. e., for , the phononic spectrum is totally destroyed giving rise to extremely localized modes, i. e., . 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 with those obtained breaking explicitly the spatial translational invariance, i. e., crossover in 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 , in the protocol, with , in the protocol. In this paper, we provide evidences for a crossover from extended to localized modes at low-frequencies as temperature decreases towards that is in agreement un recent works Wang et al. (2019); Coslovich et al. (2018). We also showed that, in analogy with Ref. Paoluzzi et al. (2019), can be employed for mapping structural into dynamical properties. In particular, the degree of localization measured through is regulated by the proliferation of dynamical heterogeneous regions in the protocol. The protocol allows to define a length scale that is an external and tunable parameter. We can thus study just looking at the solution of , with a generic observable. As it has been shown in Ref. Paoluzzi et al. (2019), mirrors the behavior of the dynamical length scale Biroli et al. (2006). Here we showed that from we can extract the behavior of that is in agreement with those observed through other observables Paoluzzi et al. (2019).
As a future direction, it would be interesting to study the density of states in systems thermalized at parental temperature close to with a fraction of pinned particles. In this way, in the spirit of early works on random pinning Cammarota and Biroli (2013, 2012); Szamel and Flenner (2013); Kob et al. (2012); Gokhale et al. (2014); Karmakar and Parisi (2013); Kob and Berthier (2013); Ozawa et al. (2015), it would be possible to take access to the properties of glassy modes towards the Kauzmann temperature.
MP acknowledges the financial support of the Joint Laboratory on “Advanced and Innovative Materials”, ADINMAT, WIS-Sapienza .
- Binder and Kob (2011) K. Binder and W. Kob, Glassy materials and disordered solids: An introduction to their statistical mechanics (World Scientific, 2011).
- Leuzzi and Nieuwenhuizen (2007) L. Leuzzi and T. M. Nieuwenhuizen, Thermodynamics of the glassy state (CRC Press, 2007).
- Berthier and Biroli (2011) L. Berthier and G. Biroli, Rev. Mod. Phys. 83, 587 (2011).
- Debenedetti and Stillinger (2001) P. G. Debenedetti and F. H. Stillinger, Nature 410, 259 (2001).
- Cavagna (2009) A. Cavagna, Physics Reports 476, 51 (2009).
- Kittel (2005) C. Kittel, Introduction to solid state physics (Wiley, 2005).
- Zeller and Pohl (1971) R. C. Zeller and R. O. Pohl, Phys. Rev. B 4, 2029 (1971).
- Gurarie and Chalker (2003) V. Gurarie and J. T. Chalker, Phys. Rev. B 68, 134207 (2003).
- Baity-Jesi et al. (2015) M. Baity-Jesi, V. Martín-Mayor, G. Parisi, and S. Perez-Gaviro, Phys. Rev. Lett. 115, 267205 (2015).
- Lerner et al. (2016) E. Lerner, G. Düring, and E. Bouchbinder, Phys. Rev. Lett. 117, 035501 (2016).
- Mizuno et al. (2017) H. Mizuno, H. Shiba, and A. Ikeda, Proceedings of the National Academy of Sciences 114, E9767 (2017).
- Angelani et al. (2018) L. Angelani, M. Paoluzzi, G. Parisi, and G. Ruocco, Proceedings of the National Academy of Sciences 115, 8700 (2018).
- Wang et al. (2019) L. Wang, A. Ninarello, P. Guan, L. Berthier, G. Szamel, and E. Flenner, Nature communications 10, 26 (2019).
- Cammarota and Biroli (2013) C. Cammarota and G. Biroli, The Journal of chemical physics 138, 12A547 (2013).
- Cammarota and Biroli (2012) C. Cammarota and G. Biroli, Proceedings of the National Academy of Sciences 109, 8850 (2012).
- Szamel and Flenner (2013) G. Szamel and E. Flenner, EPL (Europhysics Letters) 101, 66005 (2013).
- Kob et al. (2012) W. Kob, S. Roldán-Vargas, and L. Berthier, Nature Physics 8, 164 (2012).
- Gokhale et al. (2014) S. Gokhale, K. H. Nagamanasa, R. Ganapathy, and A. Sood, Nature communications 5, 4685 (2014).
- Karmakar and Parisi (2013) S. Karmakar and G. Parisi, Proceedings of the National Academy of Sciences 110, 2752 (2013).
- Kob and Berthier (2013) W. Kob and L. Berthier, Phys. Rev. Lett. 110, 245702 (2013).
- Ozawa et al. (2015) M. Ozawa, W. Kob, A. Ikeda, and K. Miyazaki, Proceedings of the National Academy of Sciences 112, 6914 (2015).
- Paoluzzi et al. (2019) M. Paoluzzi, L. Angelani, G. Parisi, and G. Ruocco, arXiv (2019).
- Lerner and Bouchbinder (2017) E. Lerner and E. Bouchbinder, Phys. Rev. E 96, 020104 (2017).
- Kob et al. (1997) W. Kob, C. Donati, S. J. Plimpton, P. H. Poole, and S. C. Glotzer, Phys. Rev. Lett. 79, 2827 (1997).
- Bernu et al. (1987) B. Bernu, J. P. Hansen, Y. Hiwatari, and G. Pastore, Phys. Rev. A 36, 4891 (1987).
- Grigera and Parisi (2001) T. S. Grigera and G. Parisi, Phys. Rev. E 63, 045102 (2001).
- Rosenfeld and Tarazona (1998) Y. Rosenfeld and P. Tarazona, Molecular Physics 95, 141 (1998).
- LaÄeviÄ et al. (2003) N. LaÄeviÄ, F. W. Starr, T. B. SchrÃ¸der, and S. C. Glotzer, The Journal of Chemical Physics 119, 7372 (2003).
- Bonnans et al. (2006) J.-F. Bonnans, J. C. Gilbert, C. Lemaréchal, and C. A. Sagastizábal, Numerical optimization: theoretical and practical aspects (Springer Science & Business Media, 2006).
- Bell and Dean (1970) R. Bell and P. Dean, Discussions of the Faraday society 50, 55 (1970).
- Grigera et al. (2003) T. Grigera, V. Martin-Mayor, G. Parisi, and P. Verrocchio, Nature 422, 289 (2003).
- Coslovich et al. (2018) D. Coslovich, A. Ninarello, and L. Berthier, arXiv preprint arXiv:1811.03171 (2018).
- Biroli et al. (2006) G. Biroli, J.-P. Bouchaud, K. Miyazaki, and D. R. Reichman, Phys. Rev. Lett. 97, 195701 (2006).