# Depletion of Two-Level Systems in Ultrastable Computer-Generated Glasses

Dmytro Khomenko, Camille Scalliet, Ludovic Berthier, David Reichman, Francesco Zamponi

► To cite this version:

Dmytro Khomenko, Camille Scalliet, Ludovic Berthier, David Reichman, Francesco Zamponi. Depletion of Two-Level Systems in Ultrastable Computer-Generated Glasses. Physical Review Letters, 2020, 124 (22), 10.1103/PhysRevLett.124.225901. hal-02738157

HAL Id: hal-02738157

https://hal.science/hal-02738157

Submitted on 2 Jun 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Depletion of two-level systems in ultrastable computer-generated glasses

Dmytro Khomenko*,¹,² Camille Scalliet*,³,\* Ludovic Berthier,⁴,⁵ David R. Reichman,² and Francesco Zamponi¹

¹Laboratoire de Physique de l'Ecole Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université de Paris, 75005 Paris, France

²Department of Chemistry, Columbia University, New York, NY 10027, USA

³DAMTP, Centre for Mathematical Sciences, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, United Kingdom

⁴Department of Chemistry, University of Cambridge, Lensfield Road, Cambridge CB2 1EW, United Kingdom

⁵Laboratoire Charles Coulomb (L2C), Université de Montpellier, CNRS, 34095 Montpellier, France

(Dated: June 17, 2020)

Amorphous solids exhibit quasi-universal low-temperature anomalies whose origin has been ascribed to localized tunneling defects. Using an advanced Monte Carlo procedure, we create in silico glasses spanning from hyperquenched to ultrastable glasses. Using a multidimensional path-finding protocol, we locate tunneling defects with energy splittings smaller than $k_B T_Q$, with $T_Q$ the temperature below which quantum effects are relevant ($T_Q \approx 1$ K in most experiments). We find that as the stability of a glass increases, its energy landscape as well as the manner in which it is probed tend to deplete the density of tunneling defects, as observed in recent experiments. We explore the real-space nature of tunneling defects, and find that they are mostly localized to a few atoms, but are occasionally dramatically delocalized.

The theory of low-temperature properties of perfect crystals stands as one of the most profound early tests of the power of quantum statistical mechanics. In particular, Debye's calculation of the observed $T^3$ behavior of the low-temperature specific heat highlighted the importance of long wavelength phonons as low energy excitations in ordered solids [1]. Given that the wavelength of populated phonon modes around $T \sim 1$ K is significantly longer than the interparticle distance in a solid, it came as a major surprise in 1971 when Zeller and Pohl [2] measured large deviations from the expected Debye behavior of the specific heat and the thermal conductivity of vitreous silica, selenium and germanium-based glasses. An explanation for this puzzling observation was almost immediately put forward by Anderson, Halperin and Varma [3], and Phillips [4]. They posited that the disorder intrinsic to amorphous solids causes their energy landscape to have many minima. Rare, nearly degenerate, adjacent local minima support tunneling defects or two-level systems (TLS) with energy splittings of the order of 1 K, which provide a large excess contribution to the specific heat and a new mode of scattering that determines the thermal conductivity. In the subsequent decades, the behavior described by Zeller and Pohl was observed in numerous other amorphous materials, and the TLS theory has withstood essentially all experimental tests [5–9]. Despite this great progress, the microscopic real-space structure of the tunneling defects remains debated, as do the factors that determine their density and distribution in amorphous solids [10–22].

A powerful platform for addressing these issues is the use of computer simulation to prepare amorphous materials in silico and to interrogate the simulated energy landscape for TLS [23]. This program was initiated by Stillinger and Weber [24, 25], then carried out more extensively by Heuer and Silbey [26–29], nearly three decades ago. Limited by the computational power and algorithms available back then, they created computer glasses with cooling rates roughly nine orders of magnitude larger than in laboratory settings. They were able to locate only a handful of TLS with the requisite tunnel splittings, necessitating uncontrolled extrapolations. The situation then improved incrementally [30–33], although the algorithmic ability to simulate glasses which are cooled in an experimentally realistic way has remained out of reach until very recently. This limited greatly our microscopic understanding of the universal anomalous thermal behavior of low-temperature glasses from a computational viewpoint.

In this work we leverage the remarkable ability of the swap Monte Carlo algorithm to produce in silico amorphous materials with fictive temperatures that range from those found in typical experiments to the significantly slower rates found in recent vapor deposition studies [34]. We find a dramatic depletion of active TLS (those with a tunnel splitting $\sim 1$ K) with decreasing quench rate, as found in recent experiments [35–39], with the notable exception of old amber glasses [40]. We use a state-of-the-art reaction path-finding protocol [41, 42] to efficiently locate double-well potentials in the multidimensional potential energy landscape, yielding a direct sampling of tunneling states with sufficient statistics to avoid any extrapolation. We determine the degree of

<footnote>\*Equal contribution to the work</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 1. (a) Glasses are prepared at equilibrium (black line) at temperatures $T_f = 0.092, 0.07, 0.062$ (bullets), from hyperquenched to ultrastable. We follow their potential energy after rapid quenches (colored lines). (b) Sketch of the potential energy landscape. Double-well potentials are detected with molecular dynamics simulations at $T_{MD} = 0.04$ (blue). (c) A detected double-well potential $V(\xi)$. The classical asymmetry $\Delta V$, activation energy $V_a$, energy barrier $V_b = V_a - \Delta V/2$, the energy levels, $\mathcal{E}_1$ and $\mathcal{E}_2$, of the ground-state doublet, and the tunnel splitting, $E = \mathcal{E}_2 - \mathcal{E}_1$, are illustrated.</figcaption></figure>

localization of individual TLS, providing a detailed, real-space understanding of how atoms participate in tunneling motion and how the thermal exploration of the energy landscape in well-annealed amorphous materials determines the density of tunneling centers.

**Glass preparation** – Past works on the landscape of low-temperature glasses focused on simple models for real materials such as NiP [26], Argon [30], and silica [32, 33]. Our goal is to understand how glass preparation affects the density of TLS, within a single model. We study a polydisperse mixture of particles interacting via an inverse power law potential [34]. Our choice is motivated by the fact that low-temperature anomalies are observed in glasses, regardless the material. Given the diversity of models, we choose one for which swap Monte Carlo enables maximally efficient thermalization on the computer over a range of temperatures at least as wide as in experiments [34].

We provide minimal details on the system and measures of equilibration (see SI for details). We simulate a non-additive polydisperse mixture of $N = 1500$ particles of mass $m$. Two particles $i$ and $j$ separated by a distance $r_{ij}$ interact via the potential

$$v(r_{ij}) = \epsilon \left( \frac{\sigma_{ij}}{r_{ij}} \right)^{12} + \epsilon F(r_{ij}/\sigma_{ij}) \quad (1)$$

only if $r_{ij} < r_{cut} = 1.25\sigma_{ij}$, $\sigma_{ij}$ being the non-additive interaction length scale. The function $F$ is a fourth-order polynomial which guarantees continuity of the potential up to the second derivative at $r_{cut}$. We characterize the physical classical dynamics of the model using molecular dynamics (MD) with energies and lengths expressed in units of $\epsilon$ and the average diameter $\sigma$, respectively. Times measured during MD simulations are expressed in units of $\sqrt{\epsilon/m\sigma^2}$. Number density is set to $\rho = 1$. The relaxation time $\tau_\alpha$ of the equilibrium fluid is measured from the self-intermediate scattering function $F_s(k = 7.0, \tau_\alpha) = e^{-1}$. The onset of glassy dynamics, signaled by deviations from Arrhenius behavior, takes place at $T_o = 0.18$, where $\tau_\alpha(T_o) \equiv \tau_o$, and the mode-coupling crossover temperature is $T_{MCT} = 0.104$ [34].

We analyze in silico glasses by preparing fully equilibrated configurations at various preparation temperatures $T_f$ using the swap Monte Carlo algorithm, which utilizes the exchange of particles’ diameters in addition to standard translational moves, leading to a massive thermalization speed-up. Our implementation is that of Ref. [43]. The configurations are then rapidly cooled to lower temperatures using regular molecular dynamics. Therefore, $T_f$ represents the “temperature at which the glass would find itself in equilibrium if suddenly brought to it from its given state,” which is precisely the definition of Tool’s fictive temperature [44]. The temperature $T_f$ characterizes the degree of stability of the glasses, see Fig. 1(a). In experiments, $T_f$ would be determined by the cooling rate [45, 46], or the substrate temperature in a vapor deposition experiment [47–51]. We present results for glasses in wide range of stabilities: poorly annealed (hyperquenched) glasses ($T_f = 0.092$ where $\log(\tau_\alpha/\tau_o) = 4.9$, slightly below $T_{MCT}$), liquid-cooled experimental glasses ($T_f = 0.07 \simeq T_g$, where $\log(\tau_\alpha/\tau_o) = 10.7$), and ultrastable glasses ($T_f = 0.062$, where $\log(\tau_\alpha/\tau_o) = 14.8$). To obtain statistically significant results, we analyze $N_g$ independent samples ($N_g = 200, 50, 15$ for increasing $T_f$).

**Landscape exploration** – We identify transitions between nearby minima, or double well potentials (DWPs) in the glasses. Briefly, starting from the configurations equilibrated at $T_f$, we run MD simulations at $T_{MD} = 0.04$, which is sufficiently low to confine each glass in a single metabasin, but high enough that the system can rapidly visit distinct minima (inherent structures) within the metabasin [23], see Fig. 1(b). Details are given in the SI.
<figure><img src="image_3.png" /><figcaption>FIG. 2. Probability distribution functions of DWP parameters as a function of glass preparation temperature $T_f$: (a) asymmetry $\Delta V$, (b) energy barrier $V_b$, (c) distance $d$ normalized by $\sqrt{PR}$, which characterizes the typical individual displacement of particles that participate actively in a double-well transition, and (d) participation ratio $PR$. Lines are a guide for the eye.</figcaption></figure>

By sampling the inherent structures along MD trajectories [30], we obtain a library of visited minima, as well as the pairs of them that are dynamically connected. We use the isoconfigurational ensemble [52]: for each of the $N_g$ independent configurations, we run up to 200 MD simulations, each initialized with different velocities. The number and duration of isoconfigurational runs is large enough for the probability distribution function (pdf) of the inherent structures potential energy, and the number of transitions, to reach stationarity. While we reach convergence of the pdfs, all possible minima are not sampled, and their numbers increase with additional runs. We however sample a significant number of minima: 13252, 26898, 848698 for $T_f = 0.062, 0.07, 0.092$, respectively. As shown below, this is enough to directly determine the density of tunneling TLS.

We select transitions between adjacent minima as described in the SI. We compute the minimum energy path connecting each pair of minima using a climbing image Nudged Elastic Band (NEB) algorithm [41, 42], which ensures the accurate determination of the saddle point, and provides a smooth energy profile. Occasionally, especially for $T_f = 0.092$, the NEB energy profile contains intermediate minima. In such cases, we apply an iterative method to split the energy profile into distinct DWPs, which are then analyzed similarly to the other ones.

We parametrize a DWP by its potential energy $V(\xi)$ along the minimum energy path, with $\xi = 0, 1$ corresponding to the two minima (we arbitrarily choose $\xi = 0$ for the deepest minimum), see Fig. 1(c). A DWP is characterized by its asymmetry $\Delta V = V(1) - V(0)$, energy barrier $V_b = V_a - \Delta V/2$, where $V_a$ is the activation energy,

<figure><img src="image_4.png" /><figcaption>FIG. 3. Cumulative distribution of energy splitting $n(E)$ divided by $E$, from hyperquenched to ultrastable glasses (top to bottom). The values $\tilde{m}$ are chosen for comparison with real materials. The plateau at small $E$ affords a direct determination of the TLS density $n_0$, which is suppressed by two orders of magnitude as glass preparation is varied.</figcaption></figure>

and the distance $d$ between minima. The distance is calculated along the reaction coordinate given by the NEB, as the sum of Euclidean distances between images of the system: $d^2 = \sum_{i,\mu} d_{i,\mu}^2$, where $d_{i,\mu}$ is the displacement of particle $i$ in direction $\mu = x, y, z$. The participation ratio, $PR = d^4 / (\sum_{i,\mu} d_{i,\mu}^4)$, characterizes the number of particles involved in the transition.

We present in Fig. 2 the statistics of the DWP parameters. The pdfs for $T_f = 0.062$ and 0.07 agree quantitatively, within noise, while we observe an evolution for $T_f = 0.092$. In particular, the pdfs of asymmetries and energy barriers are almost exponential in all glasses. The mild dependence of these tails on $T_f$ may stem from the fact that the sampling temperature $T_{MD}$ sets a limit on the DWPs that can be detected, independently of $T_f$ (see SI for the effect of $T_{MD}$). While the distribution of energy asymmetry is exponential at high energies $\Delta V$, it becomes flat at low energies where TLS are found (see Fig. S6). The pdfs of distances (not shown) and participation ratios vary more significantly between $T_f \le 0.07$ and $T_f = 0.092$. Since $d \propto \sqrt{PR}$, an increase of $PR$ will affect the distribution of $d$. We instead study the pdf of $d/\sqrt{PR}$. This quantity can be interpreted as an average displacement of the particles that participate in the transition. On average, the number of particles involved in DWPs is larger in poorly annealed glasses, while the displacements of individual particles remain comparable. To our knowledge, the dependence of the quench rate on DWPs classical parameters has not yet been reported. Note that TLS typically correspond to DWPs with very small $\Delta V$ and relatively large $V_b$. The tunnel splitting stems from non-trivial correlations between the classical parameters [30], thus the pdfs of classical parameters alone are not informative on quantum tunneling (see SI).

**Density of two-level systems** – At sufficiently low temperatures, thermal activation over the energy barrier $V_b$
is suppressed, and quantum tunneling becomes important [53]. In our analysis, we reduce the problem to a one dimensional effective Schrödinger equation along the reaction coordinate. Following Vineyard [54], the effective mass remains $m$, with a reaction coordinate $x \in [0, d]$. Using the normalized variable $\xi = x/d$, and scaling energies by $\epsilon$, we obtain

$$-\frac{\hbar^2}{2md^2\epsilon} \partial_{\xi}^2 \Psi(\xi) + V(\xi)\Psi(\xi) = \mathcal{E}\Psi(\xi), \quad (2)$$

where the “quantumness” of the problem is controlled by the dimensionless mass $\tilde{m} = m\sigma^2\epsilon/\hbar^2$ (see Fig. 3). In general, the Laplacian should take into account curvature effects, which we neglect here. The potential $V(\xi)$ obtained from the NEB is defined in $\xi \in [0, 1]$. To solve Eq. (2), we extrapolate it outside this interval (see SI).

We solve Eq. (2) for all DWPs. We compute the first two energy levels, $\mathcal{E}_1$ and $\mathcal{E}_2$, of the double well and define the tunnel splitting $E = \mathcal{E}_2 - \mathcal{E}_1$. We illustrate in Fig. 1(c) the two energy levels and tunnel splitting of a DWP. The tunnel splitting $E$ is the relevant parameter for low-temperature properties. The transitions that occur by quantum tunneling have a tunnel splitting $E \sim T$ [6]. These particular DWPs are called tunneling two-level systems (TLS).

We characterize the distribution of TLS using a cumulative distribution of tunnel splittings $n(E)$, defined as the number of DWPs with tunnel splitting smaller than $E$, normalized by the number of particles $N$ in the glass, and the number of independent samples $N_g$. In TLS theory, $n(E)$ can be expanded as $n(E) \simeq n_0E + O(E^2)$ for small $E$, the specific heat at low temperature is linear with $T$, and $n_0$ enters the prefactor [3, 6].

In order to estimate the density $n_0$ of TLS and its dependence on glass stability, we plot $n(E)/E$ as a function of tunnel splitting $E$ in Fig. 3. All curves indicate a saturation to a plateau value, $n_0$, at low $E$. The existence of a plateau value demonstrates our ability to directly estimate the density of TLS $n_0$ without any uncontrolled extrapolation. Data on tunneling rates [55] and distribution of tunneling matrix elements $\Delta_0$ are presented in the SI.

Our key result is that the TLS density $n_0$ (as estimated by $n(E)/E$ in the range $10^{-3} - 10^{-4}$) decreases by two orders of magnitude from hyperquenched to ultrastable glasses. To our knowledge, this constitutes the first numerical evidence for a significant suppression of TLS with increasing glass stability.

**Microscopic nature of TLS** – How many particles are involved in the tunneling motion of a TLS? [30–33] We analyzed how the participation ratio of transitions correlates with the tunnel splitting $E$. We find that the participation ratio of low-temperature active TLS with tunnel splittings $E \sim 10^{-3} - 10^{-4}$ varies from 1 to 200 (see SI). The higher participation ratios ($PR \sim 200$) are found in hyperquenched glasses, while in ultrastable

<figure><img src="image_5.png" /><figcaption>FIG. 4. Snapshots of TLS with low tunnel splitting $E$ for $T_f = 0.092$ and $\tilde{m} = 30000$. (a) $PR \approx 126$ with $E = 8.9 \times 10^{-5}$. (b) $PR \approx 1.6$ and similarly low $E = 5.4 \times 10^{-5}$. The size and color of particles are proportional to their displacement between the initial and final configurations of the TLS, normalized to the highest displacement.</figcaption></figure>

glasses the participation ratio rarely exceeds $\sim 30$. We provide systematic numerical evidence that TLS active at low temperature are typically very localized, but occasionally associated with collective excitations. We provide two snapshots in Fig. 4, corresponding to a collective TLS (left) and a very localized TLS (right) identified in a hyperquenched glass.

**Discussion** – Our study of tunneling TLS in a simple computer model demonstrates their importance to understand low-temperature glass anomalies. We show that the density $n_0$ of TLS directly controls the linear temperature dependence of the specific heat at low temperatures. Several recent works advocated the idea that quantized low-frequency harmonic modes alone could explain this behavior [10, 18, 56–60]. These soft modes are known for our glasses [61], but we find that their contribution to the low-temperature specific heat is subdominant (see SI), suggesting that the specific heat of structural glasses is dominated by tunneling TLS, as originally proposed in [3, 6].

To relate our data to experiments we convert simulation units into physical ones. The temperature scale below which quantum effects become important is obtained by comparing the thermal wavelength to the interparticle distance: $T_Q = \frac{2\pi\hbar^2}{m\sigma^2k_B}$. We consider DWPs with $E < k_B T_Q$ as low-temperature active TLS, and their total number is $n_{active} = NN_g n(E = k_B T_Q)$.

A detailed analysis on experimental comparisons is presented in the SI. We first consider Argon parameters: $\sigma = 3.4 \times 10^{-10}$m, $\epsilon/4 = 1.65 \times 10^{-21}$J, $m = 6 \times 10^{-26}$kg [30]. This gives $T_g \sim 32$K, $T_Q \sim 0.73$K, and $\tilde{m} \sim 4000$. For this choice, we estimate from Fig. 3 $n_0^{sim} \sim 4, 0.4, 0.04$ for increasing glass stability. This gives $n_0^{exp} \sim 10^{49}, 10^{48}, 10^{47}$ J$^{-1}$m$^{-3}$. Active TLS have $E < k_B T_Q = 0.0015\epsilon$ and we find $n_{active} = 1008, 291, 61$ such TLS for $T_g = 0.092, 0.07, 0.062$, respectively.
A second choice motivated by NiP metallic glasses [26] would be to use Nickel as a reference, for which $\sigma = 2.21 \times 10^{-10}$m, $\epsilon = 6.14 \times 10^{-20}$J, $m = 1.02 \times 10^{-25}$kg [62]. In this case, we have $T_g \sim 298$ K, $T_Q \sim 0.9$ K, and $\tilde{m} \sim 30000$. For this value of $\tilde{m}$, we find $n_0^{sim} \sim 60, 6, 0.6$ for $T_f = 0.092, 0.07, 0.062$, yielding $n_0^{exp} \sim 10^{50}, 10^{49}, 10^{48}$ J$^{-1}$m$^{-3}$. Active TLS have $E < k_B T_Q = 0.0002\epsilon$ and we find $n_{active} = 248, 46, 28$ such TLS for $T_g = 0.092, 0.07, 0.062$, respectively.

Experimentally, a value of $n_0 \sim 10^{46}$J$^{-1}$m$^{-3}$ is reported [6, 7]. Our estimation for $T_f = 0.07$ is larger by a factor $\sim 10^2 - 10^3$ (Argon and Nickel, respectively). It is difficult to rationalize this discrepancy but we can offer several hypothesis. One possibility is that we include in our estimates all DWPs detected at temperature $T_{MD} = 0.04 \gg T_Q$, while experimentally, the glass is directly brought to $T_Q$ and only a small fraction of TLS that lie at the bottom of the energy metabasin would be excited. Furthermore, not all TLS would tunnel on the relevant timescales: it is known that $n_0 \sim \log(\tau)$ where $\tau$ is the observation time at $T_Q$ [6]. This should persist up to the timescale of complete exploration of the energy landscape. Our exploration protocol at $T_{MD} \gg T_Q$ artificially sets $\tau$ larger than this cutoff (see SI for the tunneling rates of TLS). Another explanation could be that our model is too simple to describe real molecular materials, for example network glasses. Since we analyze a single model, the fundamental question of universality in TLS density remains unanswered. Analyzing different glass-forming models will be crucial to answer this question.

The reduction of $n_0$ by two orders of magnitude when moving from hyperquenched to ultrastable glasses is robust and in good agreement with recent experiments [36, 39]. We show that for a given glass-forming model, glass preparation affects dramatically the density of TLS. Our results demonstrate that glass ultrastability (rather than potential anisotropy of the vapor-deposited samples) is responsible for the depletion of TLS.

We thank M. Ediger, D. Rodney and W. Ji for discussions. This project received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program, Grant No. 723955 – GlassUniversality (FZ), and Grant No. 740269 – ADNeSP (Prof. M. E. Cates), and from the SIMONS Foundation (#454933, LB, #454955, FZ, #454951 DR). CS acknowledges the Fondation l'Oréal.

<footnote>* cs2057@cam.ac.uk</footnote>

<footnote>[1] N. W. Ashcroft and N. D. Mermin, *Solid State Physics* (Thomson Learning, 1976).</footnote>

<footnote>[2] R. C. Zeller and R. O. Pohl, Phys. Rev. B **4**, 2029 (1971).</footnote>

<footnote>[3] P. W. Anderson, B. I. Halperin, and C. M. Varma, *Philosophical Magazine* **25**, 1 (1972).</footnote>

<footnote>[4] W. Phillips, *Journal of Low Temperature Physics* **7**, 351 (1972).</footnote>

<footnote>[5] M. T. Loponen, R. C. Dynes, V. Narayanamurti, and J. P. Garno, *Physical Review B* **25**, 1161 (1982).</footnote>

<footnote>[6] W. A. Phillips, *Reports on Progress in Physics* **50**, 1657 (1987).</footnote>

<footnote>[7] J. F. Berret and M. Meißner, *Zeitschrift für Physik B Condensed Matter* **70**, 65 (1988).</footnote>

<footnote>[8] A.-M. Boiron, P. Tamarat, B. Lounis, R. Brown, and M. Orrit, *Chemical physics* **247**, 119 (1999).</footnote>

<footnote>[9] A. L. Burin, J. M. Leveritt III, G. Ruyters, C. Schötz, M. Bazrafshan, P. Fassl, M. von Schickfus, A. Fleischmann, and C. Enss, *EPL (Europhysics Letters)* **104**, 57006 (2013).</footnote>

<footnote>[10] Y. M. Galperin, V. Karpov, and V. Kozub, *Advances in Physics* **38**, 669 (1989).</footnote>

<footnote>[11] S. N. Coppersmith, *Phys. Rev. Lett.* **67**, 2315 (1991).</footnote>

<footnote>[12] A. J. Leggett, *Physica B: Condensed Matter* **169**, 322 (1991).</footnote>

<footnote>[13] A. Burin and Y. Kagan, *Physics Letters A* **215**, 191 (1996).</footnote>

<footnote>[14] D. R. Reichman, P. Neu, and R. J. Silbey, *Molecular Crystals and Liquid Crystals* **Science and Technology. Section A. Molecular Crystals and Liquid Crystals** **291**, 65 (1996).</footnote>

<footnote>[15] V. Lubchenko and P. G. Wolynes, *Phys. Rev. Lett.* **87**, 195901 (2001).</footnote>

<footnote>[16] V. Gurarie and J. T. Chalker, *Phys. Rev. B* **68**, 134207 (2003).</footnote>

<footnote>[17] V. Lubchenko and P. G. Wolynes, *Proceedings of the National Academy of Sciences* **100**, 1515 (2003).</footnote>

<footnote>[18] D. A. Parshin, H. R. Schober, and V. L. Gurevich, *Phys. Rev. B* **76**, 064206 (2007).</footnote>

<footnote>[19] D. C. Vural and A. J. Leggett, *Journal of Non-Crystalline Solids* **357**, 3528 (2011).</footnote>

<footnote>[20] A. J. Leggett and D. C. Vural, *The Journal of Physical Chemistry B* **117**, 12966 (2013), pMID: 23924397.</footnote>

<footnote>[21] D. Zhou and A. J. Leggett, arXiv:1510.05538 (2015).</footnote>

<footnote>[22] H. M. Carruzzo and C. C. Yu, *Phys. Rev. Lett.* **124**, 075902 (2020).</footnote>

<footnote>[23] A. Heuer, *Journal of Physics: Condensed Matter* **20**, 373101 (2008).</footnote>

<footnote>[24] F. H. Stillinger and T. A. Weber, *Physical Review A* **25**, 978 (1982).</footnote>

<footnote>[25] T. A. Weber and F. H. Stillinger, *Physical Review B* **32**, 5402 (1985).</footnote>

<footnote>[26] A. Heuer and R. J. Silbey, *Phys. Rev. Lett.* **70**, 3911 (1993).</footnote>

<footnote>[27] A. Heuer and R. J. Silbey, *Phys. Rev. B* **49**, 1441 (1994).</footnote>

<footnote>[28] D. Dab, A. Heuer, and R. Silbey, *Journal of Luminescence* **64**, 95 (1995).</footnote>

<footnote>[29] A. Heuer and R. J. Silbey, *Phys. Rev. B* **53**, 609 (1996).</footnote>

<footnote>[30] F. Demichelis, G. Viliani, and G. Ruocco, *PhysChemComm* **2**, 20 (1999).</footnote>

<footnote>[31] J. Reinisch and A. Heuer, *Phys. Rev. B* **70**, 064201 (2004).</footnote>

<footnote>[32] J. Reinisch and A. Heuer, *Phys. Rev. Lett.* **95**, 155502 (2005).</footnote>

<footnote>[33] T. Damart and D. Rodney, *Phys. Rev. B* **97**, 014201 (2018).</footnote>

<footnote>[34] A. Ninarello, L. Berthier, and D. Coslovich, *Physical Review X* **7**, 021039 (2017).</footnote>
[35] D. R. Queen, X. Liu, J. Karel, T. H. Metcalf, and F. Hellman, Phys. Rev. Lett. 110, 135901 (2013).

[36] T. Pérez-Castañeda, C. Rodríguez-Tinoco, J. Rodríguez-Viejo, and M. A. Ramos, Proceedings of the National Academy of Sciences 111, 11275 (2014).

[37] X. Liu, D. R. Queen, T. H. Metcalf, J. E. Karel, and F. Hellman, Phys. Rev. Lett. 113, 025503 (2014).

[38] M. D. Ediger, Proceedings of the National Academy of Sciences 111, 11232 (2014).

[39] D. R. Queen, X. Liu, J. Karel, H. C. Jacks, T. H. Metcalf, and F. Hellman, Journal of Non-Crystalline Solids 426, 19 (2015).

[40] T. Pérez-Castañeda, R. J. Jiménez-Riobóo, and M. A. Ramos, Phys. Rev. Lett. 112, 165901 (2014).

[41] H. Jónsson, G. Mills, and K. W. Jacobsen, in *Classical and Quantum Dynamics in Condensed Phase Simulations*, edited by B. J. Berne, G. Ciccotti, and D. F. Coker (World Scientific, Singapore, 1998).

[42] G. Henkelman, B. P. Uberuaga, and H. Jónsson, The Journal of chemical physics 113, 9901 (2000).

[43] L. Berthier, E. Flenner, C. J. Fullerton, C. Scalliet, and M. Singh, Journal of Statistical Mechanics: Theory and Experiment 2019, 064004 (2019).

[44] A. Q. Tool, Journal of the American Ceramic society 29, 240 (1946).

[45] M. D. Ediger, C. A. Angell, and S. R. Nagel, The Journal of Physical Chemistry 100, 13200 (1996).

[46] A. Cavagna, Physics Reports 476, 51 (2009).

[47] S. F. Swallen, K. L. Kearns, M. K. Mapes, Y. S. Kim, R. J. McMahon, M. D. Ediger, T. Wu, L. Yu, and S. Satija, Science 315, 353 (2007).

[48] K. L. Kearns, S. F. Swallen, M. D. Ediger, T. Wu, Y. Sun, and L. Yu, The Journal of Physical Chemistry B 112, 4934 (2008).

[49] M. Tylinski, Y. Chua, M. Beasley, C. Schick, and M. Ediger, The Journal of chemical physics 145, 174506 (2016).

[50] M. D. Ediger, The Journal of Chemical Physics 147, 210901 (2017).

[51] J. Ràfols-Ribé, M. Gonzalez-Silveira, C. Rodríguez-Tinoco, and J. Rodríguez-Viejo, Physical Chemistry Chemical Physics 19, 11089 (2017).

[52] A. Widmer-Cooper, P. Harrowell, and H. Fynewever, Phys. Rev. Lett. 93, 135701 (2004).

[53] M. J. Gillan, Journal of Physics C: Solid State Physics 20, 3621 (1987).

[54] G. H. Vineyard, Journal of Physics and Chemistry of Solids 3, 121 (1957).

[55] A. Andreassen, D. Farhi, W. Frost, and M. D. Schwartz, Phys. Rev. D 95, 085011 (2017).

[56] L. Gil, M. A. Ramos, A. Bringer, and U. Buchenau, Phys. Rev. Lett. 70, 182 (1993).

[57] D. A. Parshin, Phys. Rev. B 49, 9400 (1994).

[58] W. Ji, M. Popović, T. W. J. de Geus, E. Lerner, and M. Wyart, Phys. Rev. E 99, 023003 (2019).

[59] S. Franz, T. Maimbourg, G. Parisi, and A. Scardicchio, Proceedings of the National Academy of Sciences 116, 13768 (2019).

[60] M. Baggioli, R. Milkus, and A. Zaccone, arXiv:1903.05237 (2019).

[61] L. Wang, A. Ninarello, P. Guan, L. Berthier, G. Szamel, and E. Flenner, Nature communications 10, 26 (2019).

[62] W. Ching and C. C. Lin, in *Amorphous Magnetism II* (Springer, 1977) pp. 469–477.
Supplemental Material: Depletion of two-level systems in ultrastable
computer-generated glasses

Dmytro Khomenko*,¹ Camille Scalliet*,² Ludovic Berthier,³,⁴ David R. Reichman,⁵ and Francesco Zamponi¹

¹Laboratoire de Physique de l'Ecole Normale Supérieure, ENS, Université PSL,
CNRS, Sorbonne Université, Université de Paris, 75005 Paris, France

²Department of Applied Mathematics and Theoretical Physics,
University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, United Kingdom

³Department of Chemistry, University of Cambridge,
Lensfield Road, Cambridge CB2 1EW, United Kingdom

⁴Laboratoire Charles Coulomb (L2C), Université de Montpellier, CNRS, 34095 Montpellier, France

⁵Department of Chemistry, Columbia University, New York, NY 10027, USA

(Dated: June 17, 2020)

I. MODEL

We study a three-dimensional, non-additive, continuously polydisperse mixture of particles interacting via the pair potential Eq. (1) (main text). The particle diameters $\sigma_i$ are drawn from the normalized distribution $P(0.73 < \sigma < 1.62) \sim 1/\sigma^3$, with average diameter $\langle\sigma\rangle = \int P(\sigma)d\sigma = 1$ (simply denoted as $\sigma$ in the main text). We employ a non-additive cross-diameter $\sigma_{ij} = 0.5(\sigma_i + \sigma_j)(1 - 0.2|\sigma_i - \sigma_j|)$. This model is efficiently simulated with the swap Monte Carlo algorithm. The choice of particle dispersity and non-additivity make the homogeneous fluid robust against fractionation and crystallization at all temperatures numerically explorable [S1]. As in the main text, classical quantities are given in units of $\epsilon, \sigma, m$.

We characterize the dynamic slowdown of the supercooled liquid using molecular dynamics (MD) simulations. The temperature evolution of the relaxation time $\tau_\alpha$ provides three temperatures relevant to this work. The onset of glassy dynamics, signaled by deviations from Arrhenius behavior, takes place at $T_o = 0.18$, where $\tau_\alpha(T_o) = \tau_o$. The mode-coupling crossover temperature $T_{MCT} = 0.104$ is measured by fitting the relaxation time data to $\tau_\alpha \sim (T - T_{MCT})^{-\gamma}$ in a limited window of relatively high temperatures. At $T_{MCT}$, the dynamics is four orders of magnitude slower than at the onset temperature, $\tau_\alpha(T_{MCT}) \simeq 10^4\tau_o$. We take $T_{MCT}$ as a proxy for the computer glass transition, below which standard MD fail to equilibrate the homogeneous fluid. The experimental glass transition $T_g = 0.067$ is located by fitting the relaxation time data to a parabolic law, using $\tau_\alpha(T_g) = 10^{12}\tau_o$ [S1].

This glass-forming model is efficiently simulated at equilibrium with the swap Monte Carlo algorithm. We employ the swap Monte Carlo algorithm implemented in the LAMMPS package, with the optimal parameters, as described in Ref. [S2]. Supercooled liquid configurations can be generated down to $T = 0.062$, i.e. below the experimental glass transition $T_g$. In this work, we focus on configurations prepared in equilibrium conditions at temperatures $T_f = 0.092, 0.07, 0.062$. For the two lowest $T_f$ values, standard MD dynamics initialized from an equilibrium configuration is completely arrested: no diffusion is observed and the system is trapped within a glass metabasin. For the higher $T_f = 0.092$, some diffusion is observed in equilibrium; but if the system is rapidly cooled at lower temperatures, once again no diffusion is observed and a glass state is obtained. Borrowing from experimental conventions [S3], we call this the “fictive temperature” of the glasses, as discussed in the main text.

II. LANDSCAPE EXPLORATION

The amorphous configurations generated at $T_f$ are first thermalized to $T_{MD} = 0.04$ using a Berendsen thermostat. At this temperature, no diffusion is observed for all glasses. The system is then simulated in the NVE ensemble, using an integration time step of $dt = 0.01$. Configurations along the MD trajectory are used as the starting point for energy minimization via a conjugate gradient algorithm, which brings them to their inherent structure (IS), i.e. the nearest local minimum in the potential energy landscape (PEL). Note that contrary to the MD simulations, the minimization procedure does not conserve the total energy. We minimize the MD configurations obtained every 20, 10, 5 time steps for $T_f = 0.062, 0.07, 0.092$. The high frequency of minimization is chosen to identify nearby local minima, separated by an energy barrier.

From the MD simulations, we obtain a time series of inherent structures. We are interested in transitions between two different inherent structures, identified by comparing two consecutive minima. We recorded 70970, 130859, and 1593359 transitions between inherent structures, for $T_f = 0.062, 0.07, 0.092$ respectively. We wish to characterize the potential energy barriers corresponding to the transitions. This analysis is computationally costly. Given the large number of transitions detected, we chose to analyze transitions detected at least 4 times, regardless of the direction
<figure><img src="image_6.png" /><figcaption>FIG. S1. Probability distribution of the potential energy barriers $V_b$ sampled at temperatures $T_{MD} = 0.01$ and 0.04.</figcaption></figure>

of the transition ($A \rightarrow B$ or $B \rightarrow A$). The number of transitions analyzed is equal to 14195, 23535, and 117361, for $T_f = 0.062, 0.07, 0.092$ respectively.

In order to investigate the influence of the temperature $T_{MD}$ on the characteristics of the double well potentials identified, we present in Fig. S1 the probability distribution function of potential energy barriers identified with $T_{MD} = 0.01$, and 0.04 (used in the main text). The sampling temperature $T_{MD}$ influences the tail of the distribution only, which decays as $\exp(-V_b/T_{MD})$. We conclude that the sampling of relevant TLS is not affected by a variation of $T_{MD}$ within a reasonable interval.

### III. TUNNEL SPLITTING

For each analyzed transition, the Nudged Elastic Band (NEB) algorithm outputs a one-dimensional potential defined for the reduced coordinate $\xi$, between the two minima only. We run the algorithm using 40 images of the system. We need to extrapolate the potential in order to obtain a full double well potential. We tested various extrapolation schemes, such as parabolic fitting of the minima, and mirroring around each minima, defining $V(-\xi) = V(\xi)$ for $\xi < \xi_a$ and $V(2 - \xi) = V(\xi)$ for $\xi > \xi_a$, where $\xi_a$ is the coordinate of the potential maximum: $V(\xi_a) = V(0) + V_a$. In the main text, we used a linear extrapolation of the reaction path obtained with the NEB. Let us denote $\mathbf{r}_1$ and $\mathbf{r}_2$ the coordinates of the particles in the first two images of the system along the reaction path ($\mathbf{r}_1$ is an energy minimum). We measure the potential energy of the configuration, starting from $\mathbf{r}_1$, and moving in the direction $\mathbf{r}_1 - \mathbf{r}_2$. We perform a similar extrapolation at the other minimum. We show in Fig. 1 (main text) a double-well potential obtained from the NEB algorithm (blue part), and by linear extrapolation of the reaction path (black part). We have compared all methods and found that while each scheme gives a slightly different potential, the statistics of tunnel splittings remains unaffected by our choice. In the main text, we use the most physical scheme, namely the linear extrapolation of the reaction path.

Once the classical potential $V(\xi)$ is obtained by extrapolation as discussed above, we solve numerically the Schrödinger Eq. (2) (main text) using a standard Python package. Note that, in general, the Laplacian term should take into account curvature effects along the reaction coordinate, $\nabla^2 = \partial^2 + (\det g)^{-1/2} \partial_\xi [(\det g)^{1/2} \partial_i] + (\det g)^{-1/2} \partial_i [(\det g)^{1/2} \partial_\xi]$, where $g$ is a metric tensor and $\xi_i$ are coordinates orthogonal to $\xi$. For simplicity, we neglect these effects and use the standard second derivative along the reaction coordinate.

We present in Fig. S2 a scatter plot of the tunnel splitting $E$ as a function of participation ratio $PR$ and energy barrier $V_b$ for DWPs identified in glasses from hyper-quenched (top) to ultrastable (bottom). In Fig. S2 (a,c,e) we observe that in the relevant range of $E \lesssim 10^{-3}$ the value of $PR$ can be as large as ~200 for $T_f = 0.092$ and ~30 for $T_f = 0.062$. In Fig. S2 (b,d,f) we observe that in the same relevant range, the barrier is always $V_b \gtrsim 10^{-2} \gg E$, indicating that the relevant DWP are indeed TLS. We also checked (not shown) that in the same range, one always has $E_3 - E_1 \gg E$, where $E_3$ is the third energy level.

Finally, we checked that for most TLS in the relevant range, the wavefunctions of the first two levels are sufficiently delocalized over the two wells, indicating that tunnelling is active. However, the barrier $V_b$ has a relatively wide
<figure><img src="image_7.png" /><figcaption>FIG. S2. Scatter plot of the participation ratio PR (a,c,e) and the potential energy barrier Vb (b,d,f) versus the tunnel splitting E, of double-well potentials in glasses prepared at T<sub>f</sub> = 0.092 (a-b), T<sub>f</sub> = 0.07 (c-d), and T<sub>f</sub> = 0.062 (e-f). The data for DWP<sub>s</sub> found in the same glass sample (there are N<sub>g</sub> of them) are presented with the same color. The tunnel splittings are computed using an adimensional mass m<sub>tilde</sub> = 30000.</figcaption></figure>

distribution (see Fig. S2), and there exist TLS with large barriers and hence very small tunnelling matrix elements. For those TLS, the wavefunctions are almost localized, indicating that tunnelling is highly suppressed. These TLS would be frozen in experimental conditions. The wide distribution of V<sub>b</sub> is known to be responsible for a logarithmic growth of n<sub>0</sub> with observation time τ, n<sub>0</sub> ∝ log τ [S4], which might also provide an explanation for why our simulations overestimate n<sub>0</sub> with respect to experiments (see main text).

IV. DIMENSIONAL SCALING ANALYSIS FOR UNIT CONVERSION

The number of TLS per particle in a given glass sample with a tunnel splitting less than E, for E → 0, is n₀E. Hence, n₀ has the dimensions of an inverse energy, expressed in units of ϵ⁻¹. However, because in our simulation units the number density is ρ = 1, n₀ is also the number of TLS per unit volume, in units of ϵ⁻¹σ⁻³, i.e.

$$n_0^{\text{exp}} = n_0^{\text{sim}} \times \epsilon^{-1} \sigma^{-3}. \qquad (\text{S1})$$

For Argon, ϵ⁻¹σ⁻³ ~ 3.85 × 10⁴⁸ J⁻¹m⁻³, while for Nickel ϵ⁻¹σ⁻³ ~ 1.51 × 10⁴⁸ J⁻¹m⁻³, which allows one to convert our numerical results for n<sub>0</sub><sup>sim</sup> into experimental values for these two materials.
<figure><img src="image_8.png" /><figcaption>FIG. S3. The TLS density $n_0$, as extrapolated from $n(E)/E$ in the limit $E \to 0$, divided by $\tilde{m}$, as a function of glass stability, encoded by the value of $\log(\tau_\alpha/\tau_o)$ for the three values of $T_f$ and different $\tilde{m}$. A reduction of $n_0$ by two orders of magnitude is robustly observed, independently of $\tilde{m}$. Furthermore, we find $n_0 \propto \tilde{m}$ independently of $T_f$.</figcaption></figure>

Because the glass transition temperature in the simulation is $T_g = 0.067$, the corresponding glass transition temperature in physical units is

$$T_g = 0.067 \times \frac{\epsilon}{k_B}. \qquad (S2)$$

Finally, the temperature $T_Q$ at which quantum effects become relevant is that at which the thermal wavelength equals the interparticle distance. Since our simulation density is $\rho = 1/\sigma^3$, we get

$$T_Q = \frac{2\pi\hbar^2}{m\sigma^2 k_B} = \frac{2\pi}{\tilde{m}} \frac{\epsilon}{k_B}, \qquad (S3)$$

and we recall that $\tilde{m} = m\sigma^2\epsilon/\hbar^2$ is the dimensionless mass that appears in the Schrödinger equation. Finally, note that the relevant (active) TLS are those with $E < k_B T_Q$, and their total number in our simulation is given by

$$\begin{aligned} n_{active} &= NN_g n(E = k_B T_Q) \\ &\sim NN_g n_0^{sim} \frac{k_B T_Q}{\epsilon} = 2\pi NN_g \frac{n_0^{sim}}{\tilde{m}}. \end{aligned} \qquad (S4)$$

We find that $n_0/\tilde{m}$ is roughly constant for a given glass stability $T_f$, as shown in Fig. S3. This results implies that the choice of $\tilde{m}$ can be quite arbitrary, within a reasonable range. In particular, $\tilde{m}$ cannot be too small otherwise the condition $T_Q \ll T_g$ would be violated.

## V. VIBRATIONAL AND TLS CONTRIBUTIONS TO THE SPECIFIC HEAT

Here, we compare the contributions to specific heat coming from TLS and from harmonic vibrations around an inherent structure. To obtain the TLS contribution, for a given $T_f$ and $\tilde{m}$, we collect all the TLS found in all the $N_g$ glasses, with splitting $E_i$, and compute [S4, S5]

$$C_{TLS} = \frac{1}{NN_g} \sum_i \frac{x_i^2}{\cosh(x_i)^2}, \quad x_i = \frac{E_i}{2k_B T}. \qquad (S5)$$

To obtain the vibrational contribution, we considered a single glass prepared at $T_f = 0.062$, and all the $N_{IS}$ inherent structures found within that glass. For each inherent structure, we diagonalize the Hessian matrix of the classical potential, to obtain a set of eigenvalues (spring constants) $\kappa_\alpha$. Each of these provides a quantum harmonic oscillator contribution to the vibrational specific heat, which is

$$C_{vib} = \frac{1}{NN_{IS}} \sum_{\alpha} \frac{x_{\alpha}^2}{\sinh(x_{\alpha})^2}, \quad x_{\alpha} = \frac{\hbar\omega_{\alpha}}{2k_B T}, \quad \omega_{\alpha} = \sqrt{\frac{\kappa_{\alpha}}{m}}. \qquad (S6)$$
<figure><img src="image_9.png" /><figcaption>FIG. S4. (a) TLS contribution to the specific heat per particle, as a function of T. Dotted-dashed lines are the asymptotic low-temperature scaling, $C_{TLS} \sim T$, for $\tilde{m} = 30000$ (for $\tilde{m} = 5000$ the statistics is not sufficient for a good extrapolation). (b) Vibrational contribution to the specific heat per particle, as a function of T, for $T_f = 0.062$, averaged over inherent structures in a typical glass metabasin.</figcaption></figure>

Note that the finite size of the system imposes a cutoff on the low-frequency Debye behavior. In fact, the lowest frequency found in our system is $\hbar\omega_\alpha/k_B = 0.0038$ (in simulation units) for $\tilde{m} = 5000$, and $\hbar\omega_\alpha/k_B = 0.0016$ for $\tilde{m} = 30000$. To avoid this problem we used data of larger system with $N = 192000$ particles [S6] and extrapolate the quadratic region of the density of states to the $\omega \to 0$ limit.

Yet, from Fig. S4 we conclude that the vibrational contribution decays much faster than the TLS one upon lowering temperature, and that already at $T \sim 10^{-4}$ we have $C_{vib} \ll C_{TLS}$ for all considered values of $\tilde{m}$ and $T_f$.

## VI. TUNNELING RATES

In this section we compute the tunneling rates of the TLS identified in our simulations. This allows us to determine if the TLS with energy splitting $E \lesssim k_B T_Q$, which we considered as “active” at $T \lesssim T_Q$, can actually tunnel over experimental timescales.

The characteristic time to observe a transition from the higher energy level to the lower one can be estimated from the decay rate $\Gamma$, computed within the Wentzel-Kramers-Brillouin (WKB) approximation [S7]

$$ \Gamma = \frac{1}{m} \left[ \int_{0}^{a} \frac{dx}{p(x)} \right]^{-1} \exp \left[ -\frac{2}{\hbar} \int_{a}^{b} |p(x)| dx \right], \quad p(x) = \sqrt{2m(\mathcal{E}_2 - V(x))}, \qquad (S7) $$

where $[0, a]$ and $[a, b]$ are the width of the higher energy well and of the energy barrier, respectively, both measured at the energy level $\mathcal{E}_2$. This construction is sketched in Fig. S5(a). We present in Fig. S5(b-c) the scatter plots of the tunneling rate $\Gamma$ of a DWP as a function of its energy splitting $E$. Each panel corresponds to a different glass stability, from ultrastable (b) to hyperquenched (d).

To compare these results to experimental timescales, we consider the parameters for NiP ($\tilde{m} = 30000$), for which $\epsilon = 6.14 \times 10^{-20}$ J. Together with $\hbar = 1.05 \times 10^{-34}$ m$^2$kg/s, we get a natural frequency unit $\epsilon/\hbar \sim 10^{14}$ Hz for $\Gamma$. Most TLS with $E < 10^{-3}$, which could tunnel around $T_Q$, have rates $\Gamma \ge 10^{-14}$, which corresponds to $\Gamma \ge 1$ Hz in physical units. This means that most TLS identified as active based on their energy splitting would actually tunnel over experimental timescales.

## VII. DISTRIBUTION OF TUNNELING MATRIX ELEMENT $\Delta_0$

In the standard TLS model, the Hamiltonian of a single tunneling state takes the form

$$ H = \frac{1}{2} \begin{pmatrix} \Delta & \Delta_0 \\ \Delta_0 & -\Delta \end{pmatrix} \qquad (S8) $$

in the localized representation. The diagonal splitting of the TLS can be estimated by

$$ \Delta = \Delta V + \hbar \frac{\omega_2 - \omega_1}{2}, \qquad (S9) $$
<figure><img src="image_10.png" /><figcaption>FIG. S5. (a) Definition of the parameters a and b used in the WKB calculation of tunneling rates, Eq. (S7). (b-d) Transition rates Γ as a function of energy splitting E for different glass stabilities, T<sub>f</sub> = 0.062 (b), 0.07 (c), 0.092 (d). The data for DWPs found in the same glass sample (there are N<sub>g</sub> of them) are presented with the same color. The rates are computed with NiP parameters (m̃ = 30000).</figcaption></figure>

where ω<sub>1,2</sub> are the classical frequencies in the two energy minima [S5]. The tunneling matrix element is Δ<sub>0</sub> = 2⟨ψ<sub>L</sub>|H|ψ<sub>R</sub>⟩, where |ψ<sub>R,L</sub>⟩ are the wave functions of the left and right wells. While this definition has been widely used in the theoretical literature, it is not obvious how to define the localized basis |ψ<sub>R,L</sub>⟩, and therefore not practical to obtain a numerical evaluation of Δ<sub>0</sub>. We used instead the WKB approximation to estimate the tunneling matrix element [S4]

$$
\Delta_0 \approx \bar{\mathcal{E}} \exp \left[ -\frac{1}{\hbar} \int_a^b |p(x)| dx \right], \quad \bar{\mathcal{E}} = \frac{\mathcal{E}_1 + \mathcal{E}_2}{2}, \qquad (S10)
$$

where here a, b, and p(x) are derived from Ē instead of E<sub>2</sub> as in the calculation of Γ. Note that the most important
contribution, in both cases, comes from the exponential term. We show in Fig. S6(a) the resulting cumulative density
of diagonal elements Δ (in the same format as Fig. 3 in the main text), and probability distribution function P(Δ<sub>0</sub>)
of tunneling elements, computed for each T<sub>f</sub>. Our numerical results are consistent with the TLS model prediction
P(Δ<sub>0</sub>) ~ 1/Δ<sub>0</sub>. The density of diagonal elements is constant at small Δ, also in agreement with the TLS model.
However, the condition Δ<sub>0</sub> ≪ Δ is not always satisfied in our data, and as a result the true splitting E is generally
larger than Δ, which results in a decrease of the plateau values of n(E)/E compared to n(Δ)/Δ.

SUPPLEMENTARY REFERENCES

[S1] Andrea Ninarello, Ludovic Berthier, and Daniele Coslovich, “Models and algorithms for the next generation of glass transition studies,” Physical Review X 7, 021039 (2017).

[S2] Ludovic Berthier, Elijah Flenner, Christopher J Fullerton, Camille Scalliet, and Murari Singh, “Efficient swap algorithms for molecular dynamics simulations of equilibrium supercooled liquids,” Journal of Statistical Mechanics: Theory and Experiment 2019, 064004 (2019).
<figure><img src="image_11.png" /><figcaption>FIG. S6. (a) Cumulative distribution of the diagonal splitting $\Delta$, divided by $\Delta$ for different glass stabilites. The data is shown in a way similar to Fig. 3 of the main article. (b) Probability distribution of the off-diagonal element $\Delta_0$ for different glass stabilities. Our results are consistent with the TLS model prediction $P(\Delta_0) \sim 1/\Delta_0$. Data obtained for $\tilde{m} = 30000$ (NiP parameters).</figcaption></figure>

[S3] Arthur Q Tool, “Relation between inelastic deformability and thermal expansion of glass in its annealing range,” Journal of the American Ceramic society 29, 240–253 (1946).

[S4] W A Phillips, “Two-level states in glasses,” Reports on Progress in Physics 50, 1657–1708 (1987).

[S5] P. W. Anderson, B. I. Halperin, and C. M. Varma, “Anomalous low-temperature thermal properties of glasses and spin glasses,” Philosophical Magazine 25, 1–9 (1972).

[S6] Lijin Wang, Andrea Ninarello, Pengfei Guan, Ludovic Berthier, Grzegorz Szamel, and Elijah Flenner, “Low-frequency vibrational modes of stable glasses,” Nature communications 10, 26 (2019).

[S7] Anders Andreassen, David Farhi, William Frost, and Matthew D Schwartz, “Precision decay rate calculations in quantum field theory,” Phys. Rev. D 95, 085011 (2017).