Relation between heterogeneous frozen regions in supercooled liquids and non-Debye spectrum in the corresponding glasses

TL;DR

This study links the low-frequency vibrational spectrum in glasses to heterogeneous regions, showing how parental temperature influences non-Goldstone modes and their localization, with implications for understanding glass dynamics.

Contribution

It reveals how parental temperature affects the ratio of Goldstone to non-Goldstone modes and introduces a method to estimate heterogeneity size from vibrational spectra.

Findings

Parental temperature modifies the GM/NGM ratio in vibrational spectra.

The low-frequency spectrum follows a power law with an exponent dependent on temperature.

Heterogeneity size diverges mildly as temperature approaches the dynamical crossover.

Abstract

Recent numerical studies on glassy systems provide evidences for a population of non-Goldstone modes (NGMs) in the low-frequency spectrum of the vibrational density of states $D(\omega)$. Similarly to Goldstone modes (GMs), i. e., phonons in solids, NGMs are soft low-energy excitations. However, differently from GMs, NGMs are localized excitations. Here we first show that the parental temperature $T^*$ modifies the GM/NGM ratio in $D(\omega)$. In particular, the phonon attenuation is reflected in a parental temperature dependency of the exponent $s(T^*)$ in the low-frequency power law $D(\omega) \sim \omega^{s(T^*)}$, with $2 \leq s(T^*) \leq 4 $. Secondly, by comparing $s(T^*)$ with $s(p)$, i. e., the same quantity obtained by pinning \mttp{a} $p$ particle fraction, we suggest that $s(T^*)$ reflects the presence of dynamical heterogeneous regions of size $\xi^3 \propto p$. Finally, we provide an estimate of $\xi$ as a function of $T^*$, finding a mild power law divergence, $\xi \sim (T^* - T_d)^{-\alpha/3}$, with $T_d$ the dynamical crossover temperature and $\alpha$ falling in the range $\alpha \in [0.8,1.0]$.

Full text

Introduction.

At small enough frequencies $\omega$, the density of states $D(\omega)$ of a three dimensional glass follows Debye’s law $D(\omega)\sim\omega^{2}$. This is because at large enough length scale glasses are continuum media and thus phonons dominate the low frequency spectrum Kittel (2005). However, compared with crystals, glasses show thermodynamic anomalies at low temperatures. For instance, the thermal conductivity $\kappa(T)$ scales with $T^{2}$ Zeller and Pohl (1971) instead $T^{3}$, as predicted by Debye’s law Kittel (2005). Moreover, also the specific heat $C_{v}$ below $1$ K deviates from Debye’s law acquiring a linear dependency on $T$ Zeller and Pohl (1971). Remarkably, these anomalies are shared by a broad class of glassy systems providing evidences of universality.

As it has been noticed in Ref. Gurarie and Chalker (2003), in disordered media and small $\omega$, $D(\omega)$ takes contributions from both, extended Goldstone bosons, e. g., phonons in structural glasses or spin-waves in Heisenberg spin glasses, and non-Goldstone modes, i. e., excitations that are not generated by the spontaneous symmetry breaking of a continuous symmetry. The Goldstone contribution gives rise to the Debye spectrum {{\color[rgb]{0,0,0}D(\omega)}}\sim\omega^{d-1}, with $d$ the number of spatial dimensions. The non-Goldstone sector is still soft, i.e., normal modes whose density of states vanishes as a power law {{\color[rgb]{0,0,0}D(\omega)}}\sim\omega^{s} Gurarie and Chalker (2003), but it is populated by localized modes. Only in the last few years, thanks to the possibility of eliminating Goldstone bosons from the low-energy spectrum Baity-Jesi et al. (2015); Lerner et al. (2016) or discriminating non-extended modes from the extended ones Mizuno et al. (2017), it has been possible to observe numerically the non-Goldstone sector in numerical simulations obtaining that, in agreement with arguments suggested in Ref. Gurarie and Chalker (2003), non-Goldstone modes give a contribution to $D(\omega)$ that scales with $s=4$.

In a previous work, we showed that a population of soft-localized modes with {\color[rgb]{0,0,0}D(\omega)}\sim\omega^{s(p)} and $2\leq s(p)\leq 4$ emerges in the low-frequency spectrum of a three dimensional model of glass when a fraction $p$ of particles are randomly frozen Angelani et al. (2018). In particular, the value of the effective exponent $s(p)$ starts from $s$=2 at $p$=0 and approaches $4$ above a threshold $p_{th}$ value that is of the order of $50\%$ of frozen particles.

In this paper we study the properties of the vibrational density of states of a model of glass in its inherent states, i. e., configurations that minimize the potential energy at $T=0$. Inherent states have been obtained after a fast quench from equilibrium configurations at the parental temperature $T^{*}$. In agreement with Ref. Lerner and Bouchbinder (2017), the slope in the tail of $D(\omega)$ depends on the parental temperature $T^{*}$. Moreover, as a first result, we observe a progressive attenuation of the Debye spectrum in favor of the non-Debye one as $T^{*}$ approaches from above the dynamical crossover temperature $T_{d}$.

It turns out that the Debye scaling $s$=2 holds at $T^{*}\gg T_{d}$, $s$ increases by decreasing $T^{*}$, and it saturates to $s=4$ right above $T_{d}$. This crossover between Debye to non-Debye is accompanied by a progressively localization of the normal modes below the Boson peak. We will show that it is possible to relate the suppression of extended excitations with the proliferation of spatially heterogeneous regions. We observe that the behavior of $s(T^{*})$ mirrors that of $s(p)$ observed in Ref. Angelani et al. (2018), indicating an increase in the size of frozen heterogeneous regions as $T^{*}$ decreases. This finding suggests a way for measuring the size of these regions.

In particular, we are able to perform a mapping between the properties of the inherent states of the randomly pinned system at high parental temperatures, i. e., $(T^{*}=\infty,p)$, with the inherent structures of the same system at low temperatures without frozen particles. i. e., $(T^{*},p=0)$. Defining $s(T)\equiv s(T^{*},p=0)$ and $s(p)\equiv s(T^{*}=\infty,p)$ and looking at the solution of $s(T^{*})\!=\!s(p)$, we will show that the resulting curve $p(T^{*})$ provides an estimate for a correlation length $\xi$, being $\xi\equiv\xi_{pin}=(pN/\rho)^{1/3}$. We are then able to extract the behavior of $\xi_{pin}$ as a function of $T^{*}$. It turns out that $\xi_{pin}$ is compatible with a power law divergence at $T_{d}$, $\xi_{pin}^{3}=(T^{*}-T_{d})^{-\alpha}$ and $\alpha$ not far from one.

It is worth noting that pinned particles have been intensively employed for gaining insight into the structural and dynamical properties of glassy materials in both, analytical models Cammarota and Biroli (2012); Franz (2006); Szamel and Flenner (2013) and numerical simulations Scheidler et al. (2002); Biroli et al. (2008); Ozawa et al. (2015); Kob et al. (2012); Brito et al. (2013); Nagamanasa et al. (2015); Karmakar and Parisi (2013); Chakrabarty et al. (2015); Kob and Coslovich (2014); Kob and Berthier (2013); Kim (2003). Here, in parallel with the study of $D(\omega)$ as a function of $T^{*}$, we use the method of random pinning to gain insight into the $T^{*}$ dependence of the growing correlation lengths in structural glasses through the low-frequency spectrum of $D(\omega)$.

Methods.

We consider a standard $50\!:\!50$ binary mixture composed of $N\!=\!N_{A}\!+\!N_{B}$ spherical particles at density $\rho\!=\!N/L^{3}\!=\!1$ Bernu et al. (1987); Grigera and Parisi (2001); Ninarello et al. (2017). The system is enclosed in a cubic box of side $L$ where periodic boundary conditions are considered. Particle radii are $\sigma_{A}$ and $\sigma_{B}$ with $\sigma_{A}/\sigma_{B}\!=\!1.2$ and $\sigma_{A}\!+\!\sigma_{B}\!\equiv\!\sigma\!=\!1$ Grigera and Parisi (2001). Details about numerical simulations can be found in Supplemental Material (SM). We consider hybrid Brownian/Swap Monte Carlo simulations that combine the numerical integration of the equations of motion with Swap Montecarlo moves Grigera and Parisi (2001). For computing dynamical properties, i.e., the dynamical temperature $T_{d}$, the 4-point susceptibility $\chi_{4}(t)$ Lačević et al. (2003a); Biroli and Bouchaud (2004), and the 4-point correlation $S_{4}(q,t)$ Kob et al. (1997); Lačević et al. (2003a); Flenner et al. (2014); Flenner and Szamel (2015), we consider the Brownian evolution of thermodynamically stable configurations obtained through hybrid Swap/Brownian dynamics for system sizes $N=10^{3},20^{3}$. The minimization of the mechanical energy has been performed through the Limited-memory Broyden-Fletcher-Goldfarb-Shanno algorithm Bonnans et al. (2006). We then compute the eigenvalues of the dynamical matrix, i. e., the Hessian matrix $\mathbf{M}$, and thus we obtain the spectrum of the harmonic oscillations around the inherent structure.

We also considered configurations equilibrated at high temperatures, i. e., $T^{*}\gg T_{d}$, where a finite number of particles $pN$, with $p\in[0,1[$, are maintained frozen during the minimization (see Ref. Angelani et al. (2018) for details). The eigenfrequencies are $\omega_{\kappa}^{2}=\lambda_{\kappa}$, with $\lambda_{\kappa}$ the $\kappa-$th eigenvalue of $\mathbf{M}$.

We focus our attention on the cumulative $F(\omega)$=$\int_{0}^{\omega}d\omega^{\prime}\,D(\omega^{\prime})$ of the density of states $D(\omega)$=$\mathcal{N}^{-1}\sum_{\kappa}\delta(\omega-\omega_{\kappa})$, with $\mathcal{N}$ the number of non-zero modes. We study the localization property of the normal-mode $\omega$ through the inverse of the participation ratio $\mathcal{R}(\omega)\equiv\sum_{i}|{\mathbf{e}}_{i}(\omega)|^{4}/\left(\sum_{i}|{\mathbf{e}}_{i}(\omega)|^{2}\right)^{2}$ where ${\mathbf{e}}_{i}(\omega)$ is the eigenvector of the mode $\omega$ Bell and Dean (1970). Let ${\mathbf{r}}\!\equiv\!({\mathbf{r}}_{1},...,{\mathbf{r}}_{N})$ be a configuration of the system thermalized at temperature $T^{*}$. We indicate with ${\mathbf{r}}^{0}\!\equiv\!({\mathbf{r}}_{1}^{0},...,{\mathbf{r}}_{N}^{0})$ the configuration that minimizes the mechanical energy. We also quantify the effect of the parental temperature on the inherent configuration through the distribution $P(\Delta r)=N^{-1}\sum_{i}\delta(\Delta r-\Delta r_{i})$, with $\Delta r_{i}\equiv|{\mathbf{r}}_{i}-{\mathbf{r}}_{i}^{0}|$, i. e., the total displacement covered by the particle $i$ for reaching the inherent configuration ${\mathbf{r}}^{0}$ starting from ${\mathbf{r}}$, i. e., the thermally equilibrated one. Finally, we compare the behavior of $\xi_{pin}$ with the dynamical correlation length $\xi_{dyn}$ that is obtained by fitting the four-point correlation function $S_{4}(q,\tau_{4})$ to an Ornstein-Zernike expression Kob et al. (1997). The structural relaxation time $\tau_{4}$ has been evaluated looking at the peak of the $\chi_{4}(t)$ susceptibility Lačević et al. (2003); Biroli et al. (2013). Details about the computation are provided in SM.

Results.

Let us start with discussing the effect of the parental temperature $T^{*}$ on the cumulative function. $F(\omega)$ is shown in Fig. (1a) for different $T^{*}$ and system size $N=10^{3}$. Approaching the dynamical temperature, i.e. $T^{*}/T_{d}\to 1$, the exponent of the low-frequency power law $F(\omega)\sim\omega^{s(T^{*})+1}$ increases as temperature decreases departing from the Debye value $s=2$ to higher values. A dependency of $s$ on both, the protocol adopted for cooling down the system and $T^{*}$, has been observed also in Ref. Lerner and Bouchbinder (2017). As shown in Fig. (2a), $s+1\to 5$ as $T^{*}\to T_{d}$. The exponents have been computed fitting the tail of $F(\omega)$ below the Boson peak Buchenau et al. (1984); Binder and Kob (2011) with a power law. Since the Boson peak is populated by extended modes, we select the low-frequency sector through the {{\color[rgb]{0,0,0}\mathcal{R}}} value of the mode $\omega$. As one can see in the inset of Fig. (1a), below $\omega\sim 0.04$, $\mathcal{R}$ grows as frequency decreases. Moreover, in that region, $\mathcal{R}$ grows with decreasing $T^{*}$ (the arrow goes in the direction of decreasing temperatures), indicating that low-frequency modes become more localized as temperature decreases. The situation is different above $\omega\sim 0.04$ where $\mathcal{R}$ approaches the $1/N$ limit. The increasing in $\mathcal{R}(\omega)$ and the behavior $F(\omega)\sim\omega^{1+s(T^{*})}$ on lowering frequency is consistent with the presence of soft-localized modes.

In order to gain insight into the nature of the rearrangements made by the system for reaching the inherent configuration ${\mathbf{r}}^{0}$, we have then computed the distribution $P(\Delta r)$ that is shown in Fig. (1b). The distribution becomes peaked at lower and lower $\Delta r$ values as temperature decreases indicating that particles in configurations at lower temperature turn to be more caged during minimization. The red curves in Fig. (1b) are fits to $\gamma\Delta r^{2}e^{-\Delta r^{2}/\sigma_{1}^{2}}$, with $\gamma$ and $\sigma_{1}$ adjustable parameters. One can notice the presence of non-Gaussian tails at high temperatures that progressively disappear for $T^{*}\to T_{d}$. The non-Gaussian tails indicate that particles travel long distances for reaching the optimal configuration when the parental configuration is taken at high $T^{*}$. As $T^{*}$ decreases towards $T_{d}$, particles turn to be more caged by their neighbors and they thus perform small uncorrelated displacements to find the closer inherent structure. To be more quantitative, we have also computed the true variance $\sigma_{2}$ of the distribution $P(\Delta r)$, $\sigma_{2}(T^{*})$(see Fig. (2b,2d)).

A similar phenomenology is observed looking at the system at high $T^{*}$ but including a fraction of frozen particle during the research of the inherent structure Angelani et al. (2018). In particular, when the concentration of frozen particles is large enough, moving particles are caged by the non-moving ones. In Fig. (1c) the cumulative $F(\omega)$ is shown at $T^{*}/T_{d}=3.12$ by varying the fraction of frozen particles $p$, that increases from left to right. We have thus computed the distributions $P(\Delta r)$ when a fraction $p$ of particles is maintained frozen and the results are shown in Fig. (1d). Again, the red curves are the gaussian fit. The behavior of $P(\Delta r)$ with increasing the fraction of frozen particles $p$ is qualitatively the same obtained with decreasing temperature in the un-pinned system. This is a strong indication that the same crossover from a spectrum dominated by soft-extended modes to soft-localized modes takes place in both protocols, i.e, decreasing parental temperature or increasing the fraction of pinned particles. The advantage of introducing randomly frozen particles lies in the fact that controlling $p$ we are also controlling the typical size $\xi_{pin}$ of frozen regions. In order to define $\xi$, we notice that the number of frozen particles $N_{p}=pN$ is naturally proportional to the volume of frozen particles $V_{p}$, thus $\xi^{3}\propto pN$.

In order to make quantitative progresses, we look at the curves $s(T^{*})$ and $s(p)$ obtained from $F(\omega)$, as well as at $\sigma_{1}(T^{*})$ and $\sigma_{1}(p)$ (or, equivalently, at $\sigma_{2}(T^{*})$ and $\sigma_{2}(p)$) obtained from $P(\Delta r)$. As already noticed $s(T^{*})$ increases with decreasing $T^{*}$: this behavior is reported in Fig. (2a). The dashed-red line is a fit to logistic curve. Similarly, the behaviour of $\sigma_{2}(T^{*})$ as a function of $T^{*}$ is reported in Fig. (2b). In panels (c) and (d) of the same figure we show the same observables as a function of the fraction of frozen particles $p$ for configurations thermalized well above the dynamical temperature, i. e., $T^{*}/T_{d}=3.12$.

We can thus provide a quantitative estimate of the behavior of $\xi_{pin}$ as a function of $T^{*}$ mapping the properties of the pinned system into the properties of the thermal system. In particular, we assume that, in a system where a fraction $p$ of particles are frozen randomly in space, one introduces a correlation length $\xi_{pin}\equiv\left(pN/\rho\right)^{1/3}$. For inferring the correlation length $\xi_{pin}$ in the real system, i. e., without artificially frozen particles, we invert the relation $\mathcal{O}(T^{*},p=0)=\mathcal{O}(T^{*}\gg T_{d},p)$, where $\mathcal{O}$ is a generic observable, and thus we obtain a function $p(T^{*})$ that allows to measure $\xi_{pin}^{3}(T^{*})$. The green dashed arrows in Fig. (2a, 2c) give a pictorial representation of the mapping we employ to infer $\xi_{pin}$ choosing as an observable the exponent $s$. The results of our analysis are shown in Fig. (3)-a for system sizes $N=10^{3},12^{3}$. Diamonds are obtained considering the exponents of the power laws $s(T^{*})$ and $s(p)$ as observable $\mathcal{O}$ for mapping, i.e. $\mathcal{O}\equiv s$. Circles refer to the true variance of the distribution $P(\Delta r)$, $\mathcal{O}\equiv\sigma_{2}$. Ttriangles are obtained considering the parameter $\sigma_{1}$ from the fit of $P(\Delta r)$ to a gaussian distribution ($\mathcal{O}\equiv\sigma_{1}$). The dashed curves are the power law $\xi_{pin}^{3}\sim(T^{*}-T_{d})^{-\alpha}$. The exponent $\alpha$ has been computed considering $\mathcal{O}=s$ and data set $N=10^{3}$ (cyan symbols), $N=12^{3}$ (red symbols). We then obtain $\alpha_{max}=1.0$ and $\alpha_{min}=0.8$ for $N=12^{3},10^{3}$, respectively, indicating that the exponent $\alpha$ varies in the range $\alpha\in[0.8,1.0]$. As one can appreciate, different observables $\mathcal{O}$ provide estimates for $\xi$ that are consistent with the same mild power-law divergence as $T^{*}$ decreases towards $T_{d}$.

Dynamical heterogeneities are fingerprint patterns of glassy dynamics Berthier and Biroli (2011). They suggest the existence of a dynamical correlation length $\xi_{dyn}$ that can be estimated through multi-point correlation functions Biroli and Bouchaud (2004); Biroli et al. (2006); Lačević et al. (2003); Karmakar et al. (2009); Biroli et al. (2013); Flenner et al. (2014); Flenner and Szamel (2015). We thus compute $\xi_{dyn}$ for system sizes $N=10^{3},20^{3}$ and compare it with $\xi_{pin}$ computed before. The result is shown in Fig. (3b), the green circles are $\xi_{dyn}$ for $N=10^{3}$, cyan circles refer to $N=20^{3}$. $\xi_{dyn}$ has been normalized with the value of $\xi$ at high temperature, i. e., $\xi_{\infty}\equiv\xi(T^{*}\gg T_{d})$. As one can appreciate, both the correlation length $\xi_{pin}$ and $\xi_{dyn}$ show a mild growth that is compatible with $(T^{*}-T_{d})^{-1/3}$.

Discussion.

In this paper, we have explored the properties of the low-frequency excitations in a three-dimensional model glass obtained by fast quench from a well equilibrated super-cooled liquid configuration at $T=T^{*}$. The values of $T^{*}$ spanned from high temperatures down to the dynamical transition temperature $T_{d}$. We have shown that quasi-localized soft modes progressively populate the low-frequency spectrum. The density of states of these glassy modes follows a scaling law $\omega^{s(T^{*})}$ with $2\leq s(T^{*})\leq 4$. In particular, far away from the dynamical transition, the low-frequency spectrum below the Boson peak is well described by Debye’s law, i. e., $s(T^{*})=2$. As $T^{*}$ decreases, $D(\omega)$ at small $\omega$ is still power law with an exponent that is temperature dependent and deviates from Debye’s law. In particular, $s$ starts to increase its value and $s(T^{*})\to 4$ for $T^{*}\to T_{d}$. As shown here and also before in Ref. Angelani et al. (2018), the same quasi-logistic growth of $s$ is observed when, instead of varying the parental temperature, we introduce a fraction $p$ of frozen particles. In particular, the spectrum of the low-energy excitations below the Boson peak remains gapless and progressively deviates from Debye’s law following a scaling $\omega^{s(p)}$, with $2\leq s(p)\leq 4$. In this case, $s(p)$ increases as $p$ increases and $s\to 4$ above a threshold value $p_{th}\sim 0.5$.

The emerging phenomenology is consistent with a picture of heterogeneous regions where particles experienced different mobilities Kob et al. (1997); Lačević et al. (2003b); Berthier and Biroli (2011). Considering a three dimensional system, an estimate of the typical linear size $\xi$ of these heterogeneous regions in the pinned system is provided by $p^{1/3}$. This argument for the scaling of $\xi$ together with our estimate of $p(T)$ leads to an estimate of $\xi(T)$ that is compatible with the inhomogeneous Mode-Coupling theory discussed in Ref. Biroli et al. (2006) where $\xi_{dyn}\sim(T-T_{d})^{-\nu}$ with $\nu=1/4$.

It has been shown in Refs. Karmakar et al. (2012); Gutiérrez et al. (2015) that the statistical properties of the lowest eigenfrequency of $D(\omega)$ can be employed to define a static length scale whose behavior is consistent with the point-to-set length Biroli et al. (2008, 2013); Karmakar et al. (2009); Berthier and Kob (2012). Our study shows that low-frequency modes in $D(\omega)$ not only regulate the growing of a static length, and thus the changing in the thermodynamic properties of the system, but also the growing of dynamic heterogeneous patterns.

In conclusion, we showed that $D(\omega)$ provides useful information about both, the structural properties of the glassy state through its inherent structures, and the dynamical properties of the corresponding supercooled equilibrium configurations at $T^{*}$. As a consequence, from $D(\omega)$ we can extract important information about the correlation length of the heterogeneous regions in supercooled liquids, which, most likely, are the ultimate origin of the instability giving rise to the non-Goldstone modes Schirmacher (2006); Schirmacher et al. (2007, 2008); Marruzzo et al. (2013a, b); Tomaras and Schirmacher (2013); Schirmacher et al. (2015). A recent experiment showed that $D(\omega)$ results to be modified by natural hyperaging Pogna et al. . As a future direction, it would be interesting to investigate numerically $D(\omega)$ in the aging regime for understanding the relation between non-Debye spectrum and aging in glasses.

Acknowledgments.

GP acknowledges the financial support of the Simons Foundation (Grant No. 454949) and the European Research Council (ERC) (grant agreement No [694925]). GP and MP were also supported by ADINMAT, WIS-Sapienza.