# Microscopic observation of two-level systems in a metallic glass model

Felix Mocanu, Ludovic Berthier, Simone Ciarella, Dmytro Khomenko, David

Reichman, Camille Scalliet, Francesco Zamponi

► To cite this version:

Felix Mocanu, Ludovic Berthier, Simone Ciarella, Dmytro Khomenko, David Reichman, et al. Microscopic observation of two-level systems in a metallic glass model. The Journal of Chemical Physics, 2023, 158 (1), pp.014501. 10.1063/5.0128820. hal-03920514

HAL Id: hal-03920514

https://cnrs.hal.science/hal-03920514v1

Submitted on 19 Mar 2025

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
Microscopic observation of two-level systems in a metallic glass model

Felix C. Mocanu,¹ Ludovic Berthier,²,³ Simone Ciarella,¹ Dmytro

Khomenko,⁴,⁵ David R. Reichman,⁴ Camille Scalliet,⁶ and Francesco Zamponi¹

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

²Yusuf Hamied 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, 3000 Broadway, New York, NY 10027, USA

⁵Dipartimento di Fisica, Sapienza Università di Roma, P.le A. Moro 2, I-00185, Rome, Italy

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

(*felix-cosmin.mocanu@phys.ens.fr)

(Dated: January 5, 2023)

The low-temperature quasi-universal behavior of amorphous solids has been attributed to the existence of spatially-localized tunneling defects found in the low-energy regions of the potential energy landscape. Computational models of glasses can be studied to elucidate the microscopic nature of these defects. Recent simulation work has demonstrated the means of generating stable glassy configurations for models that mimic metallic glasses using the swap Monte Carlo algorithm. Building on these studies, we present an extensive exploration of the glassy metabasins of the potential energy landscape of a variant of the most widely used model of metallic glasses. We carefully identify tunneling defects and reveal their depletion with increased glass stability. The density of tunneling defects near the experimental glass transition temperature appears to be in good agreement with experimental measurements.

# I. INTRODUCTION

The mechanical and thermal properties of glassy systems at cryogenic temperatures are determined by their low-energy excitations. In particular, at temperatures below 1 K, where quantum effects are important, the model of tunneling two-level systems (TLS) [1, 2] and its extensions [3] has proven to be remarkably successful in the prediction of the linear dependence in temperature *T* of the specific heat, the *T*² dependence of the thermal conductivity [4, 5], and the plateau in ultrasonic sound attenuation [6, 7]. Despite its success, alternative explanations for these experimental observations have been proposed [8, 9]. In particular, the interactions between TLS [10], aspects related to the scattering of phonons by TLS [11–14], and the collective TLS dynamics [15] remain the subject of debate [16]. Computational models can in principle be used to establish the relative merits and the ultimate validity of different hypotheses [17–22].

Beyond answering fundamental questions about the glassy state, the existence and nature of low-energy excitations in non-crystalline solids are of significant practical interest as well. In the case of superconducting circuits used for quantum computation, the tunneling TLS present in the amorphous dielectric layers of these devices are thought to be a major source of noise and decoherence [23]. Similarly, the relationship of TLS defects to internal friction is an important factor for reducing the optical losses in the coatings used for glass mirrors that are part of the complex assemblies of gravitational wave detectors [24, 25].

A particularly important question is what controls the density of TLS, which in turn determines the thermal and transport properties of the glass. Previous work has identified the fictive temperature (*T*<sub>f</sub>), which is “the temperature at which

the glass sample would find itself in equilibrium if suddenly brought there from a given state” [26] and thus efficiently encodes the thermal history of the glass, its stability and position in the potential energy landscape, as a crucial determinant of the TLS density. In fact, a reduction of roughly two orders of magnitude was observed in experiments on evaporated silicon thin films [27], on vapor-deposited indomethacin [28], and in numerical simulations of polydisperse soft-spheres [21], where *T*<sub>f</sub> was decreased from above (hyperquenched glass) to below (ultrastable glass) the experimental glass transition temperature *T*<sub>g</sub>.

However, the numerical results for continuously polydisperse systems [21] indicate a TLS density significantly larger than experiments, which suggests that other parameters might be relevant in its determination. Here, we address this question by taking advantage of a recent extension of the swap Monte Carlo algorithm [29, 30] to computer models of metallic glasses [31], which feature both a reduced polydispersity and more realistic attractive interactions. We find the same reduction of the TLS density as a function of *T*<sub>f</sub> as in previous work. In addition, we find an overall depletion with respect to previous simulations of highly polydisperse soft spheres. Our work thus suggests that the use of a more realistic yet relatively simple model is sufficient to bring down the absolute density of TLS, in closer agreement with experimental measurements of typical materials [2, 32]. The TLS density is not a strongly universal quantity, being determined by both the microscopic parameters of the interaction potentials and the thermal history of the specific glass sample. Physically, our results confirm that increasing glass stability considerably decreases the density of TLS, thus leading to an expected decrease in dissipation in amorphous solids at low temperatures and a concomitant improvement in the material attributes as-
sociated with a variety of practical applications.

Lastly, we revisit previous numerical protocols for TLS de-
termination [21] and provide detailed insight on how the mea-
sured density of TLS depends on the exploration protocol,
and on the properties of energy minima and barriers inside
glass metabasins. Finally, we investigate the glass vibrational
modes in order to analyze the low-frequency quasi-localized
modes (QLM) found at the harmonic level. We find a weak
correlation between the density of QLM and TLS, distinct
from the proportionality suggested recently [33]. Our results
emphasize the diversity of low-energy excitations governing
the behavior of glasses at low and cryogenic temperatures.

Our manuscript is organized as follows. In Sec. II A we present the numerical tools used to isolate and observe TLS. In Sec. III we describe the results of the potential energy landscape exploration. In Sec. IV we characterise the statistical properties of the TLSs and their temperature evolution. We discuss our results in Sec. V.

II. NUMERICAL METHODS

A. General strategy

Our ultimate goal is to detect and characterize numerically
TLS in a model metallic glass, and investigate how their num-
ber and properties evolve with glass preparation. In exper-
iments, TLS are naturally excited at low temperature where
quantum effects become relevant. However, simulating the
quantum dynamics of glasses containing thousands of parti-
cles at very low temperature is prohibitively difficult, and one
should find alternative approaches to identify TLS.

Our aim is therefore to estimate n(E), i.e. the number
of tunneling double-well potentials (transition paths connect-
ing distinct potential energy minima) with associated quantum
splitting smaller than E, counted per atom and per glass sam-
ple. The TLS model [1, 2] postulates the small-E scaling

$$
n(E) \approx n_0 E + O(E^2). \tag{1}
$$

The plateau value $n_0$ reached by $n(E)/E$ at small $E$ allows one to estimate $n_{TLS}$, the number of active TLS per atom in a typical glass sample at temperature $T_Q$ where experiments are performed, typically around 1 K,

$$
n_{TLS} = n_0 k_B T_Q . \tag{2}
$$

The quantity $n_{TLS}$ deduced from $n(E)$ from Eqs. (1) and (2) is the key factor controlling physical properties at low temperatures.

We provide a brief summary of our numerical strategy to
measure *n(E)* for glasses characterised by various fictive tem-
peratures, while more detailed explanations and analysis will
be presented in the following sections. There are four impor-
tant steps:

1. Prepare glass samples of various degrees of stability, or fictive temperatures. To do so, we generate $N_g$

independent equilibrium configurations in the super-
cooled liquid phase at different temperatures T<sub>f</sub> rang-
ing from the mode-coupling crossover temperature T<sub>met</sub>
down to the experimental glass transition temperature
T<sub>g</sub> using the swap Monte Carlo algorithm [29]. We
then form glasses by rapidly quenching the configu-
rations to a lower “exploration” temperature T<sub>exp</sub> us-
ing conventional molecular dynamics (MD). Therefore,
our glassy states at T<sub>exp</sub> would correspond to equilib-
rium supercooled liquid configurations if brought back
to T<sub>f</sub>, which allows us to identify T<sub>f</sub> with Tool’s fictive
temperature [26]. Varying T<sub>f</sub> is similar to varying the
cooling or deposition rate in experiments. Each of the
N<sub>g</sub> swap-generated equilibrium configurations defines a
“glass sample” or a “glass metabasin.”

2. Explore using MD the potential energy landscape of each glass metabasin [34, 35], which contains multiple potential energy minima, or inherent structures (IS) [36–38], separated by energy barriers. The low exploration temperature $T_{exp}$ is chosen to completely suppress particle diffusion over the simulated timescales, thus confining the exploration to a single glass metabasin defined by the initial configuration. At regular intervals during the dynamical trajectory simulated at $T_{exp}$, the potential energy is minimized to generate inherent structures. We call $n_{IS}$ the number of distinct IS sampled over a single glass metabasin, and $N_{IS}$ that sampled over all $N_g$ glasses.

3. The pairs of inherent structures visited consecutively in the trajectory are candidate double-well (DW) potentials if the transition is observed at least once in both directions. For those we compute the minimum-energy path (MEP) connecting these pairs of IS using the nudged elastic band (NEB) method [39, 40]. We obtain a library of $N_{DW}$ double-well potentials.

4. The minimum energy path provides an effective one-dimensional energy profile with two minima along which the quantum splitting, decay rate and tunneling matrix elements are estimated via a one-dimensional Schrödinger equation [17, 21]. We call $N_{DW}(E)$ the total number of double wells with a quantum splitting lower than $E$, with $N_{DW}(E \to \infty) = N_{DW}$. The quantity $N_{DW}(E \sim k_B T_Q)$ is directly related to the number of two-level systems in the glass we wish to estimate.

It is useful to decompose the number *n(E)* of low-energy
excitations per atom and per glass sample as

$$
n(E) = \frac{N_{DW}(E)}{N N_g} = \frac{1}{N} \times \frac{N_{IS}}{N_g} \times \frac{N_{DW}}{N_{IS}} \times \frac{N_{DW}(E)}{N_{DW}}, \quad (3)
$$

where *N* is the number of atoms. The standard TLS model postulates that DW potentials originate from strongly localized atomic motions with a small, finite density. The number of DW in a glass is thus extensive, yielding a finite *n(E)*. Physically, in order to be active at a given temperature *T*<sub>Q</sub> over an observation time *t*<sub>w</sub>, the TLS need not only to have the relevant energy splitting *E* ~ *k*<sub>B</sub>*T*<sub>Q</sub>, but also to have a tunneling
rate $\Gamma$ such that $\Gamma t_w \gg 1$. Because the TLS model postulates a flat distribution of $\log \Gamma$ at small $\Gamma$, the effective number of active TLS grows as $\log t_w$ at large $t_w$ [1, 2, 41]. In summary, according to the TLS model we should expect a finite $n(E) \propto k_B T_Q \log t_w$ when $N$ and $N_g$ both diverge.

Given that our TLS detection protocol is different from the experimental one, we need to discuss how each ratio in Eq. (3) is expected to behave when its denominator diverges, i.e. when we increase statistics. Let us first consider the ratio $N_{DW}(E)/N_{DW}$. This is the cumulative histogram of quantum splittings, which then converges to a finite limit when $N_{DW} \to \infty$. The TLS model then predicts that the remaining ratio $N_{DW}/(N_g \times N)$ converges to a finite value for large $N$ and $N_{IS}$. However, we discuss in Appendix A that this ratio, as measured in our numerical protocol, may have very different behavior depending on the nature of TLS and their interactions. In particular, because we work at fixed (and not very large) $N$, and because the $N_g$ glass samples are independent, the total number of sampled IS scales as $N_{IS} \propto N_g$ and the ratio $N_{IS}/N_g$ then converges to a finite value when $N_g \to \infty$ at finite $N$. But, as we will see, this value is so large (naturally scaling as $\exp N$) that we are not able to enumerate all the IS in a reasonable exploration time, and as a result $N_{IS}/N_g$ depends (albeit quite weakly) on exploration time in the classical MD simulation at temperature $T_{exp}$. Predicting the large scale behavior of $N_{DW}/N_{IS}$ is more intricate, see Appendix A for a discussion; we observe in our exploration protocol that this ratio remains of order one, and depends very weakly on the chosen exploration time.

Overall, the determination of $n_{TLS}$ weakly depends on time in both the experimental measurements where quantum dynamics is at play, and in our simulations involving classical exploration of the landscape. However, the physical origin is quite distinct in both situations, and this renders a direct quantitative comparison with experiments delicate.

In the following, we discuss in more details each step of the construction of the DW library discussing in particular in a more detailed way the behavior of the different terms in Eq. (3).

## B. Metallic glass model

We simulate a ternary mixture of Lennard-Jones particles (TLJ) in three dimensions. Two particles i and j at positions $\mathbf{r}_i$ and $\mathbf{r}_j$, separated by a distance $r_{ij} = |\mathbf{r}_i - \mathbf{r}_j|$ interact with a Lennard-Jones (LJ) pairwise potential $V(r_{ij})$, given by

$$V(r_{ij}) = 4\epsilon_{ij} \left[ \left( \frac{\sigma_{ij}}{r_{ij}} \right)^{12} - \left( \frac{\sigma_{ij}}{r_{ij}} \right)^{6} \right] + S(r_{ij}), \quad (4)$$

if the particle pair distance $r_{ij}/\sigma_{ij} = x_c$ is within a cutoff distance, $x_c < 2.5$, and have no interaction if their pair distance exceeds that.

To remove any resulting discontinuities from this truncation, we consider the smoothing polynomial [42]

$$S(r_{ij}) = 4\epsilon_{ij} \left[ C_0 + C_2 \left( \frac{r_{ij}}{\sigma_{ij}} \right)^2 + C_4 \left( \frac{r_{ij}}{\sigma_{ij}} \right)^4 \right]. \quad (5)$$

The choice of coefficients: $C_0 = 10/x_c^6 - 28/x_c^{12}$, $C_2 = 48/x_c^{14} - 15/x_c^8$ and $C_4 = 6/x_c^{10} - 21/x_c^{16}$ ensures the continuity of the potential $V(r_{ij})$ and of its first two derivatives at the cutoff $x_c$.

We consider three distinct types of particles, A (large), B (small), and C (medium), in a ratio A:B:C=4:1:1. The unit length is $\sigma_{AA}$ and the unit energy is $\epsilon_{AA}$, noted $\sigma$ and $\epsilon$ for simplicity. The interaction parameters are $\epsilon_{AB} = 1.5\epsilon$, $\epsilon_{AC} = 0.9\epsilon$, $\epsilon_{BB} = 0.5\epsilon$, $\epsilon_{BC} = 0.84\epsilon$ and $\epsilon_{CC} = 0.94\epsilon$ for the energies; and $\sigma_{AB} = 0.8\sigma$, $\sigma_{AC} = 1.25\sigma$, $\sigma_{BB} = 0.88\sigma$, $\sigma_{BC} = 1.0\sigma$ and $\sigma_{CC} = 0.75\sigma$ for the ranges. All particles have the same mass $m$. The unit time is $\tau_{LJ} = \sqrt{m\sigma^2/\epsilon}$.

This ternary mixture is an extension of the well-known Kob-Andersen model (KA) [43] with a third particle type (C) ensuring resistance against crystallization and increased efficiency for equilibration with the swap Monte Carlo (MC) algorithm, as detailed in Ref. [31] where the model is referred to as KA2. We simulate $N$ particles at number density $\rho = N/L^3 = 1.35$ in a cubic box of linear size $L$ with periodic boundary conditions. This particle density is higher than the commonly investigated value $\rho = 1.2$ [43]. Indeed, the pairwise attraction gives rise to a liquid-gas spinodal at low enough temperature, which may intersect the glass transition line [44], leading to a gas-glass instability [45]. The choice $\rho = 1.35$ ensures the stability of all the studied glasses, as revealed by positive values for the equilibrium pressure of the liquid down to the lowest temperatures investigated, $T = 0.488$. Our model size $L \sim 9.6\sigma$ is large enough to avoid finite-size effects on isolated TLS, but small enough to make it unlikely to observe more than one TLS in a single configuration.

## C. Supercooled liquid dynamics and glass preparation

Once the details of the model have been established, one has to estimate the relevant temperatures that govern its glassy dynamics. These are the onset temperature $T_o$, the mode-coupling temperature $T_{mct}$, and the laboratory glass transition temperature $T_g$.

To this aim, standard molecular dynamics has been carried out using the LAMMPS code [46]. A time-step of $dt = 5 \times 10^{-3} \tau_{LJ}$ together with a Nosé-Hoover thermostat was used to sample the canonical (NVT) ensemble.

A hybrid swap MC/MD dynamics was also carried out using a modified version of LAMMPS, as detailed in [47]. Based on previous work [31], we have found that short blocks of 10 MD steps with $dt = 5 \times 10^{-3} \tau_{LJ}$ interspersed with blocks of 2100 ($1.75 \times N$) attempted particle swaps produce the optimal speed-up of the dynamics of this system in the deeply supercooled regime. We recall that both the standard and swap dynamics sample the same equilibrium thermody-
<figure><img src="image_2.png" /><figcaption>FIG. 1. Angell plot of relaxation time $\tau_{\alpha}$ versus inverse temperature measured in physical molecular dynamics (MD, green) and equilibration particle-swap dynamics (Swap, magenta). The onset $T_o \simeq 1.0$ (grey), and mode-coupling crossover $T_{mct} \simeq 0.62$ (orange) temperatures are located with dashed lines. The parabolic fit (purple) is used to extrapolate the relaxation time data at lower temperature. The shaded region between a VFT and parabolic law extrapolation for $T_g$ locates the “glass ceiling”.</figcaption></figure>

nematics of the system, provided proper decorrelation can be achieved within the simulation time window [29].

In order to monitor structural relaxation during standard and swap MC dynamics, we have estimated the structural relaxation time $\tau_{\alpha}$ from the self-intermediate scattering function $F_s(q = 7.2, t = \tau_{\alpha}) = 1/e$. The relaxation time $\tau_{\alpha}$ is plotted against inverse temperature in the manner of an Angell plot in Fig. 1.

The onset of glassy dynamics is located at $T_o \sim 1.0$, where the relaxation time data of MD dynamics starts to deviate from a log-linear Arrhenius dependence valid at higher temperatures. To obtain $T_{mct}$, we fit the relaxation times to a power-law behavior, $\tau_{\alpha} \propto (T - T_{mct})^{-\gamma}$, in a region of moderate supercooling, which results in an empirical estimate of $T_{mct} \sim 0.62$. Finally, the laboratory glass transition temperature $T_g$ is the temperature at which $\tau_{\alpha}(T_g) = 10^{12} \tau_o$, in microscopic units. Simulations of such long time scales are not tractable, so this temperature is obtained from extrapolations of the available data. Two extrapolations have been used: the Vogel-Fulcher-Tammann law [48] and the parabolic law [49]. These extrapolations define a range for the laboratory glass transition temperature of $T_g \sim 0.46 - 0.51$ which we refer to as a “glass ceiling”, because experiments cannot easily access temperatures lower than $T_g$. Swap dynamics accelerates the relaxation by up to 8 orders of magnitude at $T_g$, and enables the preparation of deeply supercooled configurations of the TLJ model, equilibrated in the liquid state near the glass ceiling. This allows us to explore for the present model of a metallic glass a broad range of fictive temperatures.

### D. Landscape exploration via classical molecular dynamics

To explore the energy landscape we used an in-house standard MD code in order to have better control over the exploration workflow. The equations of motion were integrated with a smaller timestep of $dt = 2.5 \times 10^{-3} \tau_{LJ}$ to ensure a good quality of the structures sampled and to adequately distinguish between distinct IS. Initial velocities are chosen from the Maxwell-Boltzmann distribution at $T_{exp} = 0.4 < T_g$, and temperature is then kept constant by a Berendsen thermostat [50]. This very low exploration temperature is meant to confine the exploration to the glass metabasin selected by the initial conditions, while still allowing the system to cross enough energy barriers within the metabasin. We have carefully checked that, because standard MD is essentially arrested at those temperatures (Fig. 1), no diffusion is observed. The mean-squared displacement reaches a plateau over a microscopic time and does not grow above it. Hence, the glass metabasin exploration is fully driven by the thermal vibrations of the particles at $T_{exp}$, and no significant change in the solid structure is detected.

Once every $\tau_{cg}/dt$ MD steps, we use the current configuration as initial state to minimize the energy via a conjugate gradient algorithm, to obtain an inherent structure. We repeat this procedure $N_{cg}$ times until a large library of IS is obtained, reaching a simulation time $\tau_{total} = N_{cg} \tau_{cg}$.

Once a library of IS has been constructed, we need to select IS pairs as likely candidates to be connected by a path having a double-well potential shape. To this aim, we select pairs of IS that appear consecutively in the exploration dynamics, and record the number of transitions between IS $\alpha$ and IS $\beta$ in a matrix $T_{\alpha \to \beta}$, where the labeling is ordered by energy (see Fig. 2). If a transition occurs between $\alpha$ and $\beta$ in both directions at least once ($T_{\alpha \to \beta} \ge 1$ and $T_{\beta \to \alpha} \ge 1$), we consider that the pair $(\alpha, \beta)$ is a good candidate double-well potential and we attempt to evaluate the MEP between them.

Once a library of candidate transitions is obtained from the minima, the MEP is found using the string method of Ref. [51], which is implemented in LAMMPS as a modification of the nudged-elastic band (NEB) approach [40]. We used $N_{images} = 64$ images of the system to interpolate between the two IS. The images are connected by harmonic springs with a spring constant $\kappa = 0.125 \epsilon\sigma^{-2}$. The optimization of the string method was done using the FIRE minimizer [52, 53], while interpolation of the MEP was carried out with a cubic spline. This should result in an error scaling as $O(N_{images}^{-4})$ [54]. The location of the transition state is identified using a climbing image approach [55]. Finally, we check whether the converged MEP contains any intermediate minimum. If not, we preserve it as a ‘good’ double-well potential. Those with intermediate minima can be recurrently refined by splitting into elementary double-well potentials, if these intermediate transitions are not already sampled in the dynamics.

The parameters and results of the landscape exploration are summarized in Table I. The data is split between two datasets. Dataset 1 is a coarse sampling with a long minimization period $\tau_{cg}$ and a relatively small number of independent glass
<table><caption>TABLE I. Summary of the parameters used and the obtained statistics in the landscape exploration of the metallic glass model. Two independent data sets are used. The ratios are evaluated at the end of the exploration protocol $t = \tau_{total}$.</caption><thead><tr><th></th><th colspan="3">Dataset 1</th><th colspan="4">Dataset 2</th></tr></thead><tbody><tr><td>T<sub>f</sub></td><td>0.617</td><td>0.558</td><td>0.509</td><td>0.617</td><td>0.558</td><td>0.509</td><td>0.488</td></tr><tr><td>N<sub>g</sub></td><td>16</td><td>64</td><td>256</td><td>128</td><td>256</td><td>512</td><td>1024</td></tr><tr><td>τ<sub>cg</sub></td><td>11</td><td>17</td><td>43</td><td>2</td><td>3</td><td>18</td><td>31</td></tr><tr><td>N<sub>cg</sub></td><td>160000</td><td>40000</td><td>10000</td><td>100000</td><td>100000</td><td>100000</td><td>100000</td></tr><tr><td>τ<sub>total</sub></td><td>1760000</td><td>680000</td><td>430000</td><td>200000</td><td>300000</td><td>1800000</td><td>3100000</td></tr><tr><td>N<sub>IS</sub></td><td>36247</td><td>23287</td><td>30180</td><td>70326</td><td>88112</td><td>119159</td><td>117780</td></tr><tr><td>N<sub>DW</sub></td><td>10039</td><td>8359</td><td>12117</td><td>11841</td><td>20692</td><td>41342</td><td>38027</td></tr><tr><td>N<sub>TLS</sub></td><td>40</td><td>20</td><td>41</td><td>34</td><td>61</td><td>137</td><td>125</td></tr><tr><td>N<sub>IS</sub>/N<sub>g</sub></td><td>2.27 × 10<sup>3</sup></td><td>3.64 × 10<sup>2</sup></td><td>1.18 × 10<sup>2</sup></td><td>5.49 × 10<sup>2</sup></td><td>3.44 × 10<sup>2</sup></td><td>2.33 × 10<sup>2</sup></td><td>1.15 × 10<sup>2</sup></td></tr><tr><td>N<sub>DW</sub>/N<sub>IS</sub></td><td>2.77 × 10<sup>-1</sup></td><td>3.59 × 10<sup>-1</sup></td><td>4.01 × 10<sup>-1</sup></td><td>1.48 × 10<sup>-1</sup></td><td>2.35 × 10<sup>-1</sup></td><td>3.47 × 10<sup>-1</sup></td><td>3.23 × 10<sup>-1</sup></td></tr><tr><td>N<sub>TLS</sub>/N<sub>DW</sub></td><td>3.98 × 10<sup>-3</sup></td><td>2.39 × 10<sup>-3</sup></td><td>3.38 × 10<sup>-3</sup></td><td>2.87 × 10<sup>-3</sup></td><td>2.96 × 10<sup>-3</sup></td><td>3.31 × 10<sup>-3</sup></td><td>3.29 × 10<sup>-3</sup></td></tr></tbody></table>

metabasins $N_g$, at three preparation temperatures: $T_f = 0.509$, 0.558, 0.617. Dataset 2 is a finer and more extensive sampling, with a short minimization period $\tau_{cg}$, a large number of independent glass metabasins $N_g$ and covering four temperatures, $T_f = 0.488$, 0.509, 0.558, 0.617. The use of two data sets allowed us to validate our choices for the many parameters involved in the construction of a library of double-wells.

### III. TUNNELING STATES IN THE POTENTIAL ENERGY LANDSCAPE

#### A. Statistics of potential energy minima

Landscape exploration of large glassy systems is an arduous task as the number of minima is expected to scale exponentially in the number of particles [37, 38]. Equilibrium configurations prepared with swap Monte Carlo at a temperature $T_f$ ($T_g < T_f < T_{mct}$) are used to start long MD runs ($\sim 10^5$ steps) during which the system is quenched to a lower temperature ($T_{exp} < T_f$).

We can assign to each IS $\alpha$ its energy

$$e_{\alpha} = \frac{1}{N} \sum_{i<j} V_{ij} (|\mathbf{r}_{i}^{(\alpha)} - \mathbf{r}_{j}^{(\alpha)}|), \quad (6)$$

and compute the Euclidean distance between a pair $(\alpha, \beta)$ of IS as

$$d_{\alpha\beta} = \sqrt{\frac{1}{N} \sum_{i=1}^{N} |\mathbf{r}_{i}^{(\alpha)} - \mathbf{r}_{i}^{(\beta)}|^2}, \quad (7)$$

where $\mathbf{r}_{i}^{(\alpha)}$ is the position of the $i$th particle in the $\alpha$th IS.

Fig. 2 gives a pictorial overview of the results of the exploration process for four distinct glass metabasins at the lowest $T_f = 0.488$. Each row corresponds to a glass metabasin, ordered from top to bottom by increasing number of IS sampled during exploration. From left to right, the following quantities are displayed: the matrices of the energy differences, of transitions, and of Euclidean distances between IS, and the probability density function (PDF) of the IS energies $e_\alpha$ observed during exploration. The IS in the matrices are ordered

by increasing energy $e_\alpha$. The pairs of minima corresponding to TLS are indicated by red circles in all the matrices, they typically are found near the diagonal, which means that the two minima forming a TLS are close in energy.

The most striking aspects that emerge from these figures are (i) the large heterogeneity of the number of IS (we find between $n_{IS} \sim 6$ and $n_{IS} \sim 600$ IS in a single metabasin) observed even at this low $T_f$ near the estimated $T_g$, (ii) the fact that individual IS are clustered both in energy and in Euclidean distance, as shown, e.g., by the block structure of the $d_{\alpha\beta}$ matrix [56–58], (iii) the fact that the number of transitions recorded during exploration scales roughly as the number of minima, which is apparent from the sparse and linear scaling nature of non-zero elements of the transition matrix, and (iv) the fact that low energy IS are sampled repeatedly within the metabasin and the PDF decays with increasing IS energy, as expected from a Boltzmann distribution.

To obtain more quantitative insight, we first perform exploration runs with $\tau_{cg} = dt$, i.e. minimizing the potential energy at every MD step. The number of distinct IS per glass is counted, and a distribution of persistence times ($\tau_p$) is inferred from the resulting time series: here $\tau_p$ is the total time during which the energy minimization always ends in the same IS, before a new IS is found in the next step. The cumulative distribution function (CDF) of $\tau_p$, for different $T_f$, is shown in Fig. 3a. The shape of the CDF indicates that the distribution of $\tau_p$ is bimodal, with a first peak around $\tau_p \sim dt$ (i.e. a single MD step) and a second one around $\tau_p \sim 10^2 dt$. A more careful analysis of the IS time series reveals that the peak at short $\tau_p$ corresponds to processes during which the system jumps from a low-energy IS to a higher energy IS, where it stays for one or two steps, before transiting back to the low-energy IS. These processes are associated with strongly asymmetric double-well potentials, which do not give rise to TLS.

Fig. 3b shows the logarithm of the number of distinct IS obtained during exploration of a given glass metabasin, $n_{IS}(t)$, averaged over glasses, as a function of the number of steps for different $T_f$. We chose to average the logarithm of $n_{IS}$ in order to make sure that the average would not be dominated by rare metabasins with many IS. We however found that taking the logarithm of the average, i.e. $\log[N_{IS}(t)/N_g]$, yields similar
<figure><img src="image_3.png" /><figcaption>FIG. 2. Heterogeneity of the potential energy landscape of four different glass metabasins labelled in red on the left hand side. Every row of panels corresponds to a glass prepared at $T_f = 0.488$ and sampled at $T_{exp} = 0.4$. In each glass, distinct inherent structures (IS) are labelled with greek letters $\alpha = 0, 1, \dots, n_{IS} - 1$ in order of increasing energy $e_\alpha$. Each column of panels corresponds to a different observable, labeled at the top in blue. The first three columns correspond to matrices where values are evaluated for every pair of minima in a given glass. A consistent color map (viridis) is used for all matrices, with purple highlighting a low value and yellow a large one. First column (a,e,i,m): matrix of IS energy difference $|e_\alpha - e_\beta|$ in units of $10^{-3}\epsilon$. Second column (b,f,j,n): matrix $T_{\alpha \to \beta}$ counting the number of observed transitions from $\alpha$ to $\beta$, shown in a log-scale. Third column (c,g,k,o): matrix of Euclidean distance $d_{\alpha\beta}$ in units of $\sigma$. Fourth column (d,h,l,p): probability density function of IS energies. Red circles indicate the transitions that correspond to tunneling two-level systems.</figcaption></figure>

results.

We find a sub-linear power-law dependence of the number of minima on the exploration time, i.e.

$$ \langle \log(n_{IS}) \rangle \propto \beta \log(t/\tau_0) \approx \log[N_{IS}(t)/N_g], \quad (8) $$

with an exponent $\beta \in [0.4, 0.8]$ depending on the preparation temperatures considered. As $T_f$ decreases, the typical persistence time increases, and the total number of minima per glass, at a given time $t$, decreases. This reflects the already well-documented trend that glass metabasins tend to become much simpler in more stable glasses [56].

The results of Fig. 3a suggest that one can make the exploration process more efficient by performing the energy minimization (that is computationally costly) after every $\tau_{cg}/dt$ MD steps only, instead of after every step. We have chosen $\tau_{cg}$ in our final data production runs for landscape exploration as a compromise between accuracy and efficiency. A low value of $\tau_{cg}$ results in a more accurate count of minima and transitions, with fewer cases where intermediate minima need to be resolved. However, it can be computationally demanding to minimize too frequently as new minima appear in significant numbers only after $\sim 100-1000$ steps. The values $\tau_{cg}$
<figure><img src="image_4.png" /><figcaption>FIG. 3. (a) Cumulative distribution function of the persistence time $\tau_p$ between distinct consecutive IS in the landscape exploration when minimising every MD step at $T_{exp} = 0.4$. (b) Number of newly found IS versus time, geometrically averaged over 88 distinct glasses. Colors code for the glass preparation temperature $T_f$.</figcaption></figure>

employed are provided in Table I.

We do not observe any saturation of the total number of IS, $N_{IS}(t)$, with growing exploration time $t$. As a result, we cannot perform an exhaustive search of all the IS within a glass metabasin, as can be expected for a system of mesoscopic size (1200 atoms). Yet, the average number of minima found in each glass basin, $N_{IS}/N_g$, is strongly depleted with increased glass stability, at each fixed $t$. We can conclude that more stable glasses have fewer IS, provided the exploration time of the metabasin exceeds $t \sim \tau_{LJ}$. The time unit $\tau_{LJ}$ seems to be close to the characteristic time for moving outside of the basin of a single IS. In the real system at 1 K one might also expect that only a subset of minima are visited. However, in that case, exploration is dominated by tunneling through barriers and not by crossing them, as in our classical simulations.

## B. Identification of double-well potentials

We now want to understand whether our landscape exploration protocol is able to identify all the relevant double well potentials. To this aim, for each DW potential, we measure the first time $t$ at which the corresponding IS pair passes all the filters to be considered a candidate transition: both IS have been included in the library, and the condition $T_{\alpha \to \beta} \ge 1$

<figure><img src="image_5.png" /><figcaption>FIG. 4. (a) Number of inherent structures (green), double-wells (orange) and two-level systems (purple) sampled in glasses ($T_f = 0.488$) as a function of exploration time. The number of IS grows as a power-law while the growth of DW is slower at long times. TLS follows the same time evolution as DW. (b) Ratio of the number of DW to that of IS as a function of exploration time. The two panels share the horizontal axis.</figcaption></figure>

and $T_{\beta \to \alpha} \ge 1$ are met. We then compute the number of detected DW and TLS (defined in Sec. III C) as a function of $t$, shown in Fig. 4a for glasses prepared at the lowest temperature $T_f = 0.488$. Similar results are obtained at all temperatures. We observe that the growth of the number of detected transitions with exploration time is slower than that of the number of IS, and is somehow intermediate between a power-law with a small exponent and a logarithmic asymptotic behavior. The growth in the number of TLS qualitatively matches that of the DW.

Assuming logarithmic behavior, comparison with Eq. (8) gives $N_{DW}(t) \propto \log t \propto \log N_{IS}(t)$, which suggests that we are able to identify the proper elementary excitations of our system, as discussed in Appendix A, see the discussion after Eq. (A2). A faster growth of $N_{DW}$ with $N_{IS}$, which remains compatible with our data, would suggest that interactions between elementary excitations play a role. This result also implies that even if we are unable to reach a proper saturation of the IS library with exploration time, we are instead able to achieve a much more satisfactory saturation of the DW library. This suggests that the new IS discovered at large times correspond to combinations of already detected excitations.

We also note that the connectivity of the explored IS, as encoded in the ratio $N_{DW}/N_{IS}$ shown in Fig. 4b, remains of order one and has a mild dependence on the exploration time. Note that this mild time-dependence is still compatible with
<figure><img src="image_6.png" /><figcaption>FIG. 5. Schematic representation of a double-well potential illustrating the approximate MEP obtained numerically (orange) and its extrapolation beyond the two minima (blue), along with the forward barrier $V_a$, asymmetry $\Delta V$, barrier height $V_b = V_a - \Delta V/2$ (black arrows), and the curvature of the potential-energy surface at the lowest minimum $\omega_1$, highest minimum $\omega_2$, and transition state $\omega_{TS}$ (blue curved arrows). The quantum splitting $E$ is the difference between the first two quantum energy levels $\epsilon_1$ and $\epsilon_2$ (dashed red).</figcaption></figure>

$N_{DW}$ being logarithmic in time, and $N_{IS}$ being a power-law with a small exponent. We provide the equivalent of Fig. 4 for the polydisperse soft sphere model in the Appendix, Fig. 12.

### C. Estimation of quantum tunneling

Two-level systems within our glassy models correspond to tunneling DW potentials with a low quantum splitting. The temperature scale below which quantum effects are important is $T_Q \sim 1$ K in experiments and can be obtained from comparing the interparticle distance with the thermal wavelength, leading to:

$$T_Q = \frac{2\pi\hbar^2}{m\sigma^2 k_B} \qquad (9)$$

The MEP obtained via a converged optimization of the string method, with no intermediates, results in a 1D potential energy profile, as schematized in Fig. 5. By convention, the lowest energy minimum is taken to be the one on the left. From this simplified 1D potential, several quantities can be defined: the forward barrier, $V_a$, is the energy difference between the transition state and the lowest energy minimum; the asymmetry, $\Delta V$, is the difference in energy between the two minima, and the barrier height can be obtained from these two quantities: $V_b = V_a - \frac{1}{2}\Delta V$.

As a first approximation, we reduce the tunneling problem to the effective one dimensional potential of the MEP. Considering a normalized reaction coordinate $\xi \in [0, 1]$ which corresponds to the fraction of the Euclidean distance $d$ along the minimum energy path, we can write the resulting Schrödinger equation as

$$-\frac{\hbar^2}{2md^2\epsilon}\partial_{\xi}^2\psi(\xi) + V(\xi)\psi(\xi) = \epsilon\psi(\xi). \qquad (10)$$

An interpolated MEP obtained from the NEB calculation is discretized by $d\xi = d/N_\xi$, where $N_\xi = 2000$, using a cubic spline. The Laplacian is evaluated with a 5-point finite difference stencil. We use $\hbar = 1$ to evaluate the quantum properties of the DW. This leads to a dimensionless effective mass parameter $\tilde{m}$ which controls the “quantumness” of the problem, as first introduced by Vineyard [59],

$$\tilde{m} = m \frac{\epsilon \sigma^2}{\hbar^2} \qquad (11)$$

The MEP and the effective mass $\tilde{m}$ define a 1D potential for the DW excitation. We assume that the tunneling problem can be treated in a 1D approximation when the relevant classical path, the MEP, is nearly independent from the others [18]. We have calculated the eigenvalues and eigenvectors of the Hessian along the MEP to check that orthogonally to the MEP, the dynamics is harmonic and independent from the reaction coordinate. The results are summarised in the Appendix and in Fig. 13.

The $\tilde{m}$ parameter depends on the choice of units. Nevertheless, it can be tuned over a wide range of values ($10^2 - 10^5$) without changing the qualitative behavior of the quantum splitting distribution obtained from our simulations [21]. In this case, for the two sets of units detailed in Sec. IV B, namely argon (Ar) and nickel-phosphorous (NiP), the value of $\tilde{m}$ is 1200 and 5250, respectively. The corresponding values of $k_B T_Q$ are 0.005$\epsilon$ (for Ar) and 0.0012$\epsilon$ (for NiP).

The five smallest eigenvalues $\epsilon_1, \dots, \epsilon_5$ are then evaluated using ARPACK [60]. The splitting is given by the first two energy levels, $E = \epsilon_2 - \epsilon_1$. We consider as tunneling TLS all the filtered and tunneling DW that have a splitting $E < T_Q$. Following our notations, their total number is $N_{DW}(T_Q)$.

We can also attempt to map the DW profile into an actual two-level system. 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 (12)$$

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 (13)$$

where $\omega_1$ and $\omega_2$ are the characteristic frequencies of the two minima of the double-well potential, see Fig. 5. The decay rate of the TLS can be evaluated in the Wentzel–Kramers–Brillouin (WKB) approximation as

$$\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], \qquad (14) \\ p(x) = \sqrt{2m(\epsilon_2 - V(x))}.$$

The tunneling matrix element $\Delta_0$ can be obtained from the effective one-dimensional potential via the WKB approximation as

$$\Delta_0 \approx \bar{\epsilon} \exp \left[ -\frac{1}{\hbar} \int_a^b |p(x)| dx \right], \quad \bar{\epsilon} = \frac{\epsilon_1 + \epsilon_2}{2}. \qquad (15)$$
<figure><img src="image_7.png" /><figcaption>FIG. 6. (a) Scatter plot of the quantum splitting $E$ versus $\sqrt{\Delta^2 + \Delta_0^2}$. Colors code for different preparation temperatures $T_f$. (b) Cumulative distribution function of the tunneling matrix element $\Delta_0$. Inset: decay rate $\Gamma$ versus $\Delta_0$, following a $\Gamma \propto \Delta_0^2$ scaling (line). Cumulative distribution functions of (c) the diagonal splitting $\Delta$ and (d) the quantum splitting $E$. Dashed lines indicate $T_Q = 0.0012$. The linear behavior in (c,d) below $T_Q$ directly validates the hypothesis of the TLS model. Data in all panels were calculated using $\tilde{m} = 5250$, corresponding to NiP units.</figcaption></figure>

The integration limits $a$ and $b$ correspond to values of the effective reaction coordinate $x$ where a certain energy level crosses the double-well potential curve at either side of the barrier. This level is $\epsilon_2$ for the calculation of $\Gamma$, and $\tilde{\epsilon}$ for $\Delta_0$.

As such, the quantities $\Delta$ and $\Delta_0$ are expected to approximately obey the relationship $E \approx \sqrt{\Delta^2 + \Delta_0^2}$. The quantum splitting $E$ is plotted against $\sqrt{\Delta^2 + \Delta_0^2}$ in Fig. 6a, showing the expected correlation at all temperatures. Small deviations from $E \approx \sqrt{\Delta^2 + \Delta_0^2}$ are observed at low splittings $E$. These deviations could appear due to the breakdown in the various approximations that are made to calculate $E$, $\Delta$ and $\Delta_0$. The magnitude of the deviations becomes important as $E \sim 1/N$ suggesting that finite-size effects could influence the absolute magnitude of very low splittings.

Fig. 6b shows the cumulative probability distribution function $F(\Delta_0)$ of the tunneling matrix element $\Delta_0$. The TLS model predicts a flat PDF for $\log\Delta_0$, which corresponds to $F(\Delta_0) \sim \log\Delta_0$ at small $\Delta_0$. Our data instead suggest that $F(\Delta_0) \sim \Delta_0^\phi$ with a very small exponent $\phi \approx 0.1$, corresponding to a PDF $p(\Delta_0) \sim \Delta_0^{\phi-1}$. Such behavior was discussed in more detail elsewhere [2, 61].

The decay rate $\Gamma$ is also expected to approximately obey $\Gamma \propto \Delta_0^2 E$ due to Landau-Zener tunneling [62], which is seen in the inset of Fig. 6b. Finally, the CDFs of $\Delta$ and $E$ are shown in Fig. 6c and Fig. 6d, respectively. Below $T_Q$, both the CDFs of $\Delta$ and $E$ are linear, as expected in the TLS model, and change very little with $T_f$.

The ratio of tunneling TLS to generic double-wells, $N_{DW}(T_Q)/N_{DW}$, seems to be roughly constant $\sim 0.3\%$, at all temperatures. Because the function $N_{DW}(E)/N_{DW}$ converges to a stable finite limit when $N_{DW} \to \infty$, the ratio $N_{DW}(T_Q)/N_{DW}$ is also stable in both datasets presented in Table I, indicating that it is not particularly sensitive to the details of landscape exploration.

# IV. MICROSCOPIC PROPERTIES OF TWO-LEVEL SYSTEMS

## A. Statistical properties of two-level systems

The data from the extensive exploration of the energy landscape is summarized in Fig. 7, which presents the statistics of different observables of tunneling DW potentials in relation to
<figure><img src="image_8.png" /><figcaption>FIG. 7. Statistics of tunneling double-well potentials versus the estimated quantum splitting E: (a) asymmetry ΔV, (b) barrier height Vb, (c) tunneling decay rate Γ, and (d) participation ratio PR. All values were calculated with a reduced mass $\tilde{m} = 5250$, corresponding to NiP units. Colors code for different preparation temperatures Tf. Solid lines correspond to average values at a given E obtained via locally weighted regression.</figcaption></figure>

their estimated quantum splitting E and the preparation temperature Tf of the glass. The DW asymmetry, ΔV, is shown in Fig. 7a where one can notice a threshold value ΔV ~ 10⁻², below which the distribution of quantum splittings becomes much broader. Asymmetry between energy minima seems to be on average higher for higher preparation temperature Tf.

The statistics of the energy barrier, Vb, versus E are shown in Fig. 7b. Our dataset probes a range of two orders of magnitude in barrier heights between sampled tunneling DW. In general, DW with a low splitting E have a relatively high barrier Vb. DW sampled in lower Tf glasses typically have a higher barrier than those found in glasses prepared at high Tf.

The statistics of the tunneling decay rate Γ (≈ Δ₀² × E) are shown in Fig. 7c. The distribution of tunneling decay rates Γ becomes narrower as the tunnel splitting E decreases. On average the decay rate Γ decreases with the quantum splitting E. Roughly one third of the sampled tunneling DW decay in more than 1 hour.

Our numerical study provides information on the microscopic properties of the transition sampled. From a structural point of view, the degree of localization can be evaluated from the atomic displacement $d_{\alpha\beta}^i = |\mathbf{r}_i^{(\alpha)} - \mathbf{r}_i^{(\beta)}|$ of atom i between

the two IS α and β, by calculating the participation ratio (PR)

$$PR = \frac{[\sum_i (d_{\alpha\beta}^i)^2]^2}{\sum_i (d_{\alpha\beta}^i)^4}, \qquad (16)$$

defined such that $1 \le PR \le N$. The participation ratio indicates the number of atoms involved in a transition.

The statistics of PR versus E are shown in Fig. 7d. Most of the displacements observed in our database are localized, and they have an average participation ratio PR ~ 8. Even the most delocalized DW that pass the filtering procedure do not have a PR much larger than 100. We see from Fig. 7d that on average, the PR decreases with decreasing Tf. This can be rationalized in terms of displacements that tend to be more local in more stable glasses.

However, examples of delocalized DW (PR ~ 100) with a low E do exist, but they are found only in glasses prepared at high Tf ≈ Tmct. Such delocalized excitations correspond to very small individual displacements and have a small barrier. This is a similar pattern to what was observed also in the atomic tunneling of complex crystalline defects such as kinks in metallic Cu [63].

We have also looked at the scatter plot of Euclidean distances versus the quantum splitting which is shown in Fig. 14
<figure><img src="image_9.png" /><figcaption>FIG. 8. Visualization of atomic displacements between the two minima forming a double-well potential (DW) with low quantum splitting. Only the particles that displace the most are shown, and the rest are made transparent and faded in the background for clarity. (a) A localized two-level system in which a single particle moves in an almost frozen structure, *PR* ~ 1. (b) A typical two-level system, with *PR* ~ 7. (c) A delocalized two-level system with *PR* ~ 83; these cases are rare and mostly occur in less stable glasses *T*<sub>f</sub> = 0.617.</figcaption></figure>

in the Appendix. The patterns of real-space atomic displacements between the configurations of the two IS forming a DW potential are illustrated in Fig. 8. Atoms are shown with a radius proportional to their interaction range $\sigma$ and colored by $d_{\alpha\beta}^i$, the magnitude of their displacement between the two minima. In order to ease visualization only the particles that displace the most are shown, with the rest faded into the background for clarity. For clarity, only the first particles up to PR, to the nearest integer, are highlighted. In general, displacements in our database of DW are highly localized and involve only several atoms. We nevertheless can distinguish three classes of TLS, similar to a recent study on amorphous Si [64]. The most localized DW potentials have a $PR \sim 1$ such as the one shown in Fig. 8a. These would correspond to defect hopping within the material. A more typical TLS is shown in Fig. 8b, it has a $PR \sim 7$. It is still one atom that moves significantly, but this triggers the rearrangement of nearest neighbours as well. Although very rare, we also do find delocalized excitations, see Fig 8c, where nearly 10% of the particles move a relatively small distance throughout the whole structure. These tend to occur in less stable, higher $T_f \approx T_{mct} \approx 0.62$ glasses.

## B. Depletion of tunneling two-level systems with increasing glass stability

The standard TLS model predicts a plateau $n(E)/E \to n_0$ as $E \to 0$, where $n(E)$ is the cumulative distribution of quantum splittings. The results obtained from our energy landscape exploration of the ternary LJ model are shown in Fig. 9a. In agreement with the tunneling model, we observe a plateau for $E < T_Q$ (vertical dashed) at all $T_f$, and a peak at $E \sim 5 \times 10^{-2}$. In addition to the low-*E* plateau behavior observed at each temperature, we identify a clear depletion of $n(E)/E$ as $T_f$ decreases. This agrees with the previous computational estimation of the TLS density as a function of glass stability in

<figure><img src="image_10.png" /><figcaption>FIG. 9. (a) Cumulative distribution of the quantum splitting *E* for all tunneling double-well potentials, as a function of glass stability $T_f$ ($\tilde{m} = 5250$). Plateau values indicate $n_0$. (b) Density of two-level systems $n_0$ in ternary Lennard-Jones (TLJ) and polydisperse soft-sphere (PSS) glasses as a function of stability, as encoded in $\tau_\alpha(T_f)/\tau_o$. Data for two sets of units: NiP (circle, $\tilde{m} = 5250$), and Ar (square, $\tilde{m} = 1200$). PSS data from Ref. [21]. To compare the models we have fixed the exploration time $t = 10^3 \tau_{LJ}$ to estimate $n_0$.</figcaption></figure>
polydisperse soft-spheres (PSS) [21].

In order to compare the TLJ data with the PSS model we have plotted $n_0/\tilde{m}$ versus the stability of the glass, as encoded in the relaxation time ratio $\log(\tau_\alpha(T_f)/\tau_o)$, in Fig. 9b. In this representation, glass stability increases from left to right. The relaxation time $\tau_\alpha(T_f)$ is either directly measured, or estimated by extrapolating the data to lower temperature using a parabolic law, which was shown to perform well such extrapolation [65]. The $n_0$ values for PSS are from Ref. [21], while we have refined our $\tau_\alpha$-extrapolation using the latest data obtained via long MD simulations of PSS [66]. The range of glass stabilities explored in the ternary model is roughly two orders of magnitude smaller than in PSS. This is due to the particle-swap algorithm thermalizing more efficiently continuously polydisperse mixtures compared to ternary mixtures. To make a meaningful comparison of the defects in both models, we have estimated the TLS density $n_0$ after a fixed exploration time of $t = 10^3 \tau_{LJ}$. This is the longest exploration time that can be probed from our simulations, in both models, at all $T_f$. This restriction has reduced the range of $n_0$ values showed previously for the PSS model [21]. Nevertheless, at this fixed exploration time, the two models show a similar depletion of two-level systems, with the TLJ model showing a slightly steeper depletion. As evidenced from Table I, this depletion of tunneling defects seems to be driven by a reduction in the number of IS available in a glass metabasin as $T_f$ decreases. The curves are shown for different values of $\tilde{m}$ to probe quasi-universality by emulating different material properties within the same model. This confirms that the precise value of $\tilde{m}$ does not affect our main conclusion.

The choices of $\tilde{m}$ prompt a short discussion on the units. Two sets of physical units have been investigated in order to make a connection to experimental observations. The first one corresponds to physical units that mimic pairwise interactions between argon (Ar) atoms [18], and has the following parameters: $\sigma = 3.405 \times 10^{-10}$ m, $\epsilon/k_B = 125.2$ K and a mass of $m = 6.634 \times 10^{-26}$ kg. The unit of time is $\tau_{LJ} = 2.1$ ps. In these units, the numerically estimated glass transition temperature is $T_g \approx 59$ K, $T_Q \approx 0.658$ K and $\tilde{m} \approx 1200$. For the preparation temperatures $T_f = 0.488, 0.509, 0.558$ and $0.617$ one gets $n_0^{\text{sim}} \sim 0.04, 0.09, 0.14$ and $0.674 \epsilon^{-1} \sigma^{-3}$ which corresponds to $n_0^{\text{exp}} \sim 4.05 \times 10^4, 1.19 \times 10^4, 2.61 \times 10^4$ and $8.95 \times 10^4$ J$^{-1}$ m$^{-3}$.

The second set of units, used to mimic the average pairwise interactions in a NiP glass [19, 67], consists of the following parameters: $\sigma = 2.21 \times 10^{-10}$ m, $\epsilon/k_B = 934$ K and an average mass per particle of $m = 9.266 \times 10^{-26}$ kg. The unit of time is $\tau_{LJ} = 2.5$ ps. In these units, the glass transition temperature is $T_g \approx 438$ K, $T_Q \approx 1.118$ K and $\tilde{m} \approx 5250$. The defect density is $n_0^{\text{sim}} \sim 0.08, 0.19, 0.23$ and $1.84 \epsilon^{-1} \sigma^{-3}$ which corresponds to $n_0^{\text{exp}} \sim 4.24 \times 10^4, 1.11 \times 10^4, 1.53 \times 10^4$ and $1.17 \times 10^4$ J$^{-1}$ m$^{-3}$.

We find a modest variation in the absolute value of the tunneling defect density with $\tilde{m}$. The difference between the defect density of the two models (PSS and TLJ) at the same glass stability is similar in magnitude to the difference between the defect density calculated for different values of $\tilde{m}$ and the data

is insufficient to establish any clearer pattern.

A reduction of the number of defects could be expected in the TLJ model due to the change from a polydisperse system in which every particle is different, to one with only three types of indistinguishable particles (A,B,C). An additional influence could come from the attractive interactions, absent in the PSS model, which have been recently suggested to have an important effect on the elastic properties of glassy solids. In particular, a reduction of the density of quasi-localized modes (QLM) was reported, as the attractive part, or “stickiness” of the pair-potential increases [68]. The decrease in QLM density can be up to an order of magnitude and was compared to the reduction in defects observed during thermal annealing. We defer the discussion of the behavior of QLM in these models and their potential relationship to $n_0$ to the next section.

### C. Comparison to quasi-localized harmonic modes

To begin evaluating the vibrational modes of the TLJ glass configurations, we first estimate their elastic moduli at different $T_f$ from athermal quasistatic deformation simulations. With an initial IS obtained through energy minimization from an equilibrium configuration, a suitable cycle of small deformations is applied, each deformation being followed by an energy minimization, until a stress-strain curve is obtained. The shear modulus $G$ is obtained from athermal quasi-static shear simulations (AQS) under Lees-Edwards boundary conditions, while the bulk modulus $K$ is obtained from quasistatic deformation with periodic boundary conditions and a hydrostatic strain. The elastic moduli are obtained from the elastic regimes of the corresponding stress-strain curves,

$$G = \frac{\sigma_{xy}}{\gamma_{xy}}, \quad K = -V \frac{\Delta p}{\Delta V}. \qquad (17)$$

Once the moduli have been determined, we calculate the transverse $c_t = \sqrt{G/(mp)}$ and longitudinal $c_l = \sqrt{(K + \frac{4G}{3})/(mp)}$ sound velocities, and from them the frequency $\omega_t$ of the first phonon,

$$\omega_t = \frac{2\pi}{L} c_t, \qquad (18)$$

and the Debye frequency,

$$\omega_D = \left( \frac{18\pi^2\rho}{2c_t^{-3} + c_l^{-3}} \right)^{\frac{1}{3}}. \qquad (19)$$

We have calculated and diagonalized the dynamical matrix to obtain the vibrational modes and their frequencies. This was done for the ensemble of structures in Table I, as well as for larger 12000-atom TLJ configurations at $T_f = 0.509, 0.558$ and $0.617$, using ensembles of 1000 configurations at each $T_f$. Non-phononic vibrational modes can be seen below the first phonon peak at $\omega_t$, when the model size is small enough to avoid significant hybridization with phonons. It is typically
<figure><img src="image_11.png" /><figcaption>FIG. 10. (a) Integrated density of states, $I(\omega)/\omega^5$ in large $N = 12000$ ternary Lennard-Jones (TLJ) glasses prepared at different $T_f$. Dashed lines represent the $A_4/5$ fits of $I(\omega)/\omega^5$ just below $\omega_t$. (b) $A_4$ versus $n_0$ for $N = 1200$ (green) and larger $N = 12000$ (orange) TLJ glasses (orange circles) and the polydisperse soft sphere (PSS) configurations from Ref. [21] (purple). The $A_4$ and $\omega_D$ values for the PSS model are from Ref. [69].</figcaption></figure>

assumed that the density of states $D(\omega)$ scales as $\omega^4$ in the low-frequency limit [70].

In Fig. 10a we show the integrated density of states $I(\omega)$, scaled by the expected $\omega^5$ for the larger 12000-atom models at each $T_f$. Provided $D(\omega) \propto A_4\omega^4$, we would expect a plateau in $I(\omega)/\omega^5$ just below $\omega_t$. We find this to be the case, however, the plateau (dashed) persists only over a limited range of frequency, below which the scaling seems to change, in agreement with recent results on the behavior of very-low frequency modes in the standard KA model [71]. The low-frequency behaviour of $D(\omega)$ is currently under intense scrutiny, and this discussion would require much more statistics than we currently possess. Our goal is simpler, as we only wish to understand whether the evolution of localized harmonic modes in $D(\omega)$ and of TLS are strongly correlated.

Moving on to the quantitative relationship between QLM and TLS, we report in Fig 10b the adimensional quantity $A_4\omega_D^5$ versus the density of TLS obtained with NiP units. Again the TLS density $n_0$ is compared after a fixed exploration time of $10^3\tau_{LJ}$ for both models. Since both quantities decrease with decreasing $T_f$, they necessarily appear correlated in such a representation, but this of course does not imply any causal relationship [72]. In particular, our data seems to preclude a direct proportionality between the two quantities of the kind

predicted in Ref. [33]. This is especially true for the more realistic metallic glass model. Interestingly, while the density of TLS of the two models is very similar for glasses of equivalent stability, we find instead an increased density of QLM in the TLJ model. These observations suggest that generic QLM are poor predictors of TLS even though both families of localized excitations are similarly sensitive to the stability of the glass.

V. CONCLUSIONS

Our analysis confirms that TLS are depleted with increasing glass stability in a second, more realistic and physically distinct glass-forming model [21]. The number of two-level systems converges relatively quickly at $\sim 0.3\%$ of the number of the double-well potentials identified during the classical exploration of the potential energy landscape. As the glass preparation temperature $T_f$ decreases, we find a slowly increasing ratio of the number of double wells over that of minima, which reflects a locally more connected energy landscape. The significant depletion of two-level systems with glass stability is in the end driven by the drastic reduction in the number of minima in a typical glass metabasin as $T_f$ decreases. The number of minima as a function of exploration time increases slowly as a sublinear power-law for all temperatures considered, with an exponent that decreases with decreasing $T_f$. This means that while the depletion of two-level systems remains robust, an exhaustive counting of both minima and defects remains elusive, even when exploration is strongly confined to a well-defined metabasin.

The exploration of low energy glass metabasins reveals a significant degree of heterogeneity. In particular, we observe very large fluctuations, of several orders of magnitude, in the number of minima present in independent glass basins, even as $T_f \approx T_g$. An analysis of the metabasins from the perspective of connected pairs of IS shows a hierarchical structure of the energy landscape, as previously observed in Refs. [56–58]. The observation of these extremely strong sample-to-sample fluctuations in glass samples deserves further investigation [73, 74].

Comparing the continuously polydisperse soft-sphere system to a more realistic ternary model like TLJ, some important differences emerge as well. The absolute number of two-level systems, as estimated from the $n_0$ values for the TLJ model, are a bit closer to those obtained in real experiments, but remain larger by roughly an order of magnitude. While in the polydisperse model every particle is different, the ternary model contains particles of different types that are indistinguishable, which should decrease the possible number of defects, hence that of minima. From this point of view, glassy materials of increased chemical complexity should naturally exhibit more defects, while monodisperse or elemental counterparts of similar kinetic stability are likely to have fewer low-energy excitations. We observe only a modest effect though, which seems consistent with the experimentally observed quasi-universality of TLS density.

Another factor for the reduced density of two-level systems is the presence of attractive interactions and their influence
on barrier heights. Control over the strength of attractive interactions may help modulate the properties of glassy materials through the reduction of low-energy excitations, but trade-offs might come in to play with other types of defects such as quasi-localized modes. Other interaction potentials could be studied as well, in order to better understand the role played by many-body interactions, anisotropy and dimensionality.

When attempting to perform a quantitative comparison of the number of TLS observed in experiments and in computer simulations, one should keep in mind that in experiments landscape exploration is driven by quantum tunneling at ~ 1 K, which leads to a characteristic log(t) dependence of the TLS density, while in our computer study the landscape exploration is driven by classical thermal fluctuations at temperatures much higher than 1 K, leading to a power-law growth of the number of minima with time. A direct comparison is then probably hindered by this important difference.

We have also explored the nature of the relationship between quasi-localized modes and two-level systems. We find the density of two-level systems to be correlated with that of quasi-localized modes, but we do not find them to be proportional or related by any causal link. At very low frequency, modes may no longer obey an $\omega^4$ scaling, although large uncertainties are present in our data. Interestingly, while the two models have a similar density of TLS, the density of QLM seems to increase for the TLJ model. This observation suggests that there is a diversity of defects present in glasses at low temperatures, which still evade an exhaustive classification.

As perspectives for further numerical studies we would like to delineate three important research directions. First, to make our main conclusion linking glass stability to the depletion of the density $n_0$ of TLS even stronger, some of the approximations used to estimate $n_0$ should be removed from the simulations. The most ambitious task would be to explore the landscape and obtain the quantum splitting for candidate TLS, from a fully quantum-mechanical calculation using path-integral methods [75, 76]. A second task would be to improve the computational tools used in this work to significantly increase the library of TLS through use of techniques in enhanced sampling and machine learning. Third, the discovery that glass metabasins contain such a large number of minima organised in a way that depends on glass stability, although hinted by previous work, requires a dedicated study to better quantify this evolution and understand how it impacts the physical properties of glassy and viscous liquid states.

## ACKNOWLEDGMENTS

We thank E. Flenner, G. Folena, M. Ozawa, J. Sethna, S. Elliott, G. Ruocco and W. Schirmacher for useful discussions. This work was granted access to the HPC resources of MesoPSL financed by the Region Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d'Avenir supervised by the Agence Nationale pour la Recherche. This project received funding from the European Research Council (ERC) under the Euro-

pean Union's Horizon 2020 research and innovation program, Grant No. 723955 - GlassUniversality (FZ), and from the SIMons Foundation (#454933, LB, #454955, FZ, #454951 DR). CS acknowledges support from the Herchel Smith Fund and Sidney Sussex College, Cambridge.

## DATA AVAILABILITY STATEMENT

The data that support the findings of this study such as the statistics of sampled double-well potentials and equilibrium configurations of the ternary Lennard-Jones model are openly available at http://doi.org/10.5281/zenodo.7117711.

## Appendix A: Large-scale behavior of $N_{DW}/N_{IS}$

We discuss here the behavior of $N_{DW}/N_{IS}$ in the limit of a large number of samples and particles. Let us consider, as an example, a large system having K independent localized excitations (i.e. DW potentials associated to the displacement of a few atoms). Note that within the TLS model, K should be proportional to the number of atoms N provided N is large enough. Each DW can be represented by a one-dimensional local reaction coordinate $x_i$, with an associated potential taken, for the sake of illustration, to be of the form

$$v_i(x_i) = \frac{1}{4!} x_i^4 + \frac{1}{2} \kappa_i x_i^2 - h_i x_i. \quad (A1)$$

If $\kappa_i^3 + 9h_i^2/8 < 0$, this potential has two local minima, leading to a DW with energy splitting $\Delta V_i(\kappa_i, h_i)$. Let us introduce a spin variable $\sigma_i = 0, 1$ if $x_i$ is in the absolute, respectively local, minimum of $v_i(x_i)$. If the excitations are diluted and thus do not interact, the total energy of the system can be written as

$$V(\{x_i\}) = \sum_{i=1}^{K} v_i(x_i) \quad \Rightarrow \quad E_{IS}(\{\sigma_i\}) = \sum_{i=1}^{K} \Delta V_i \sigma_i. \quad (A2)$$

Such a system has $N_{IS} = 2^K$ local minima, corresponding to all possible combinations of the $\sigma_i$, with energy $E_{IS}(\{\sigma_i\})$. On the other hand, we have $N_{DW} = K$, with each IS being connected to K others: in fact, the MEP associated to a composite excitation of the form $(\sigma_i = 0, \sigma_j = 0) \rightarrow (\sigma_i = 1, \sigma_j = 1)$ would be decomposed into elementary events, e.g. $(\sigma_i = 0, \sigma_j = 0) \rightarrow (\sigma_i = 0, \sigma_j = 1) \rightarrow (\sigma_i = 1, \sigma_j = 1)$, leading to the appearance of an intermediate local minimum in $(\sigma_i = 0, \sigma_j = 1)$ in the MEP. This very simple argument thus suggests a relationship of the kind $N_{DW} \approx \log N_{IS}$.

In summary, for $K \propto N$ independent excitations in each glass, and for $N_g$ glasses, we would have in total $N_{IS} \sim N_g \exp(K) \sim N_g \exp(N)$, $N_{DW} \sim N_g K \sim N_g N$, and $N_{DW}/(NN_g)$ would then converge to a finite value, as in the TLS model.

Yet, elementary excitations deform the solid matrix, resulting in long-range elastic interactions, which would lead to a more accurate representation of the form (see e.g. [73] and
references therein)

$$V(\{x_i\}) = \sum_{i=1}^{K} v_i(x_i) - \sum_{i<j} J_{ij}x_i x_j. \quad (A3)$$

For weak enough coupling $J_{ij}$, the system still has $N_{IS} = 2^K$ local minima, but the DW profile of a single elementary excitation now depends on the state of all the other excitations in the system. As a result the MEP path connecting a given pair of IS is not necessarily decomposed into elementary transitions, and the number of DW can be as large as $N_{DW} \sim N_{IS}^2$. Depending on the nature of elementary excitations and their interactions, one would then observe $\log N_{IS} \lesssim N_{DW} \lesssim N_{IS}^2$, and thus very different large scale behavior of the ratio $N_{DW}/N_{IS}$. Extracting the relevant elementary excitations in this situation will obviously be cumbersome.

We can tentatively interpret the stability dependence of our numerical data as a smooth crossover between interacting DW for poorly stable glasses with a large concentrations of excitations to non-interacting ones in stable glasses where excitations become more dilute.

## Appendix B: Data for polydisperse soft spheres

In Fig. 11 we show data for the scaling of inherent structures and DWs numbers with time, for the polydisperse soft-sphere glasses studied in [21]. We can see that the number of IS grows as a power-law, with scaling exponents changing from 0.65 to 0.17 when increasing the glass stability, which is qualitatively similar to the scaling of $N_{IS}(t)$ in TLJ glasses, shown in Fig. 3b. These exponents seem compatible with those obtained for the TLJ system.

<figure><img src="image_12.png" /><figcaption>FIG. 11. Number of IS per glass averaged geometrically averaged over all glasses for polydisperse soft spheres. Dashed lines are guides for the eye, showing power-law scaling $\sim t^\beta$, with $\beta = 0.17, 0.28, 0.65$ for $T_f = 0.062, 0.07, 0.092$ respectively.</figcaption></figure>

In Fig. 12 we show the number of detected DW as a function of time for PSS glasses. The number of DW grows slower than the number of IS, and scaling seems slower than a power-law. All these features are the same as those observed for TLJ

glasses, as shown in Fig. 4. A significant difference however is that at a similar glass stability and length of sampling there are significantly more DW per IS in the PSS model.

<figure><img src="image_13.png" /><figcaption>FIG. 12. Landscape exploration of polydisperse soft spheres prepared at $T_f = 0.07$. (a) Number of minima (green), double-wells (orange) and TLS (purple) sampled as a function of time. (b) Ratio of number of DW to number of IS as a function of time. The two panels share the horizontal axis.</figcaption></figure>

## Appendix C: Relationship between the Hessian and the minimum-energy path

We have studied the eigenvectors and eigenvalues of the Hessian for configurations along the MEP and these are shown in Fig. 13. We find a single eigenvalue (the lowest one) becomes negative along the MEP, indicating a first-order saddle point. Only several of the lowest frequencies change appreciably along the MEP. We also projected the first six eigenvector tangent to the MEP, and we find that the lowest mode contributes the most along the path, while contributions from several others become important only when approaching the minima.

This validates the picture of a network of minima where each pair that corresponds to a TLS is connected by a 1D path, the MEP. Orthogonal to this path the dynamics is quasi-harmonic and is determined by a set of frequencies that are nearly independent from the reaction coordinate.
<figure><img src="image_14.png" /><figcaption>FIG. 13. Eigenvalues and eigenvectors of the Hessian along the MEP of a typical TLS. (a) The first ten eigenvalues (in units of $\epsilon$) along the MEP, the lowest eigenvalue is shown in purple while the tenth mode is shown in yellow. (b) Projection of the first six eigenvectors on the local tangent to the MEP. The first eigenvector is shown in blue, while the next are shown in orange, green, red, purple and brown, respectively.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>FIG. 14. Euclidean distance *d* versus splitting *E*. Colors code for different preparation temperatures *T*<sub>*f*</sub>. Average values of *d* at a given *E* are shown as a solid line.</figcaption></figure>

**Appendix D: Euclidean distance between minima and quantum splitting**

The Euclidean distance between minima is a key quantity for tunneling DWs as it is strongly related to the tunneling rate. The scatter plot of Euclidean distance *d* versus the quantum splitting *E* is shown in Fig 14. The Euclidean distance between IS tends to be larger at higher *T*<sub>*f*</sub>.
[1] P. W. Anderson, B. I. Halperin, and C. M. Varma, Anomalous low-temperature thermal properties of glasses and spin glasses, *The Philosophical Magazine: A Journal of Theoretical Experimental and Applied Physics* **25**, 1 (1972).

[2] W. A. Phillips, Two-level states in glasses, *Reports on Progress in Physics* **50**, 1657 (1987).

[3] M. A. Ramos and U. Buchenau, Low-temperature thermal conductivity of glasses within the soft-potential model, *Physical Review B* **55**, 5749 (1997).

[4] R. C. Zeller and R. O. Pohl, Thermal Conductivity and Specific Heat of Noncrystalline Solids, *Physical Review B* **4**, 2029 (1971).

[5] V. Lubchenko and P. G. Wolynes, The origin of the boson peak and thermal conductivity plateau in low-temperature glasses, *Proceedings of the National Academy of Sciences* **100**, 1515 (2003).

[6] P. Esquinazi, ed., * Tunneling Systems in Amorphous and Crystalline Solids* (Springer, Berlin, Heidelberg, 1998).

[7] U. Buchenau, G. D'Angelo, G. Carini, X. Liu, and M. A. Ramos, Sound absorption in glasses, *Reviews in Physics* **9**, 100078 (2022).

[8] V. Lubchenko and P. G. Wolynes, Intrinsic quantum excitations of low temperature glasses, *Phys. Rev. Lett.* **87**, 195901 (2001).

[9] A. J. Leggett and D. C. Vural, “ Tunneling Two-Level Systems” Model of the Low-Temperature Properties of Glasses: Are “Smoking-Gun” Tests Possible?, *The Journal of Physical Chemistry B* **117**, 12966 (2013).

[10] A. L. Burin and Y. Kagan, On the nature of the universal properties of amorphous solids, *Physics Letters A* **215**, 191 (1996).

[11] D. Zhou and A. J. Leggett, Random nonlinear infinite-level-system model for amorphous solid phonon echo and saturation phenomena, arXiv:1510.05538 (2015).

[12] D. Zhou, Universal ratio of TTLS-phonon coupling constants in low-temperature amorphous solids, *Journal of Physics: Condensed Matter* **32**, 055704 (2019).

[13] H. M. Carruzzo and C. C. Yu, Why phonon scattering in glasses is universally small at low temperatures, *Phys. Rev. Lett.* **124**, 075902 (2020).

[14] H. M. Carruzzo, A. Bilmes, J. Lisenfeld, Z. Yu, B. Wang, Z. Wan, J. R. Schmidt, and C. C. Yu, Distribution of two-level system couplings to strain and electric fields in glasses at low temperatures, *Physical Review B* **104**, 134203 (2021).

[15] C. Artiaco, F. Balducci, and A. Scardicchio, Signatures of many-body localization in the dynamics of two-level systems in glasses, *Phys. Rev. B* **103**, 214205 (2021).

[16] V. Lubchenko and P. G. Wolynes, The Microscopic Quantum Theory of Low Temperature Amorphous Solids, in *Advances in Chemical Physics*, edited by S. A. Rice (John Wiley & Sons, Inc., Hoboken, NJ, USA, 2008) pp. 95–206.

[17] A. Heuer and R. J. Silbey, Microscopic description of tunneling systems in a structural model glass, *Phys. Rev. Lett.* **70**, 3911 (1993).

[18] F. Demichelis, G. Viliani, and G. Ruocco, Properties of the double well potential and relaxation processes in a model glass, *PhysChemComm* **2**, 20 (1999).

[19] J. Reinisch and A. Heuer, How Cooperative Are the Dynamics in Tunneling Systems? A Computer Study for an Atomic Model Glass, *Journal of Low Temperature Physics* **137**, 267 (2004).

[20] T. Damart and D. Rodney, Atomistic study of two-level systems in amorphous silica, *Phys. Rev. B* **97**, 014201 (2018).

[21] D. Khomenko, C. Scalliet, L. Berthier, D. R. Reichman, and F. Zamponi, Depletion of Two-Level Systems in Ultrastable Computer-Generated Glasses, *Physical Review Letters* **124**, 225901 (2020).

[22] A. Kumar, I. Procaccia, and M. Singh, Density of quasi-localized modes in athermal glasses, *Europhysics Letters* **135**, 66001 (2021).

[23] C. Müller, J. H. Cole, and J. Lisenfeld, Towards understanding two-level-systems in amorphous solids: insights from quantum circuits, *Reports on Progress in Physics* **82**, 124501 (2019).

[24] R. Birney, J. Steinlechner, Z. Tornasi, S. MacFoy, D. Vine, A. S. Bell, D. Gibson, J. Hough, S. Rowan, P. Sortais, S. Sproules, S. Tait, I. W. Martin, and S. Reid, Amorphous Silicon with Extremely Low Absorption: Beating Thermal Noise in Gravitational Astronomy, *Physical Review Letters* **121**, 191101 (2018).

[25] J. Steinlechner, I. W. Martin, A. S. Bell, J. Hough, M. Fletcher, P. G. Murray, R. Robie, S. Rowan, and R. Schnabel, Silicon-Based Optical Mirror Coatings for Ultrahigh Precision Metrology and Sensing, *Physical Review Letters* **120**, 263602 (2018).

[26] A. Q. Tool, Relation between inelastic deformability and thermal expansion of glass in its annealing range, *Journal of the American Ceramic society* **29**, 240 (1946).

[27] D. R. Queen, X. Liu, J. Karel, T. H. Metcalf, and F. Hellman, Excess specific heat in evaporated amorphous silicon, *Phys. Rev. Lett.* **110**, 135901 (2013).

[28] T. Pérez-Castañeda, C. Rodríguez-Tinoco, J. Rodríguez-Viejo, and M. A. Ramos, Suppression of tunneling two-level systems in ultrastable glasses of indomethacin, *Proceedings of the National Academy of Sciences* **111**, 11275 (2014).

[29] A. Ninarello, L. Berthier, and D. Coslovich, Models and Algorithms for the Next Generation of Glass Transition Studies, *Physical Review X* **7**, 021039 (2017).

[30] L. Berthier and D. R. Reichman, Modern computational studies of the glass transition, arXiv:2208.02206 (2022).

[31] A. D. S. Parmar, M. Ozawa, and L. Berthier, Ultrastable Metallic Glasses In Silico, *Physical Review Letters* **125**, 085505 (2020).

[32] J. F. Berret and M. Meißner, How universal are the low temperature acoustic properties of glasses?, *Zeitschrift für Physik B Condensed Matter* **70**, 65 (1988).

[33] W. Ji, Toward understanding the depletion of two-level systems in ultrastable glasses, arXiv:2112.10105 (2021).

[34] B. Doliwa and A. Heuer, Energy barriers and activated dynamics in a supercooled lennard-jones liquid, *Physical review E* **67**, 031506 (2003).

[35] R. A. Denny, D. R. Reichman, and J.-P. Bouchaud, Trap models and slow dynamics in supercooled liquids, *Physical review letters* **90**, 025503 (2003).

[36] F. H. Stillinger and T. A. Weber, Hidden structure in liquids, *Physical Review A* **25**, 978 (1982).

[37] F. Sciortino, Potential energy landscape description of supercooled liquids and glasses, *Journal of Statistical Mechanics: Theory and Experiment* **2005**, P05015 (2005).

[38] A. Heuer, Exploring the potential energy landscape of glass-forming systems: from inherent structures via metabasins to macroscopic transport, *Journal of Physics: Condensed Matter* **20**, 373101 (2008).

[39] H. Jónsson, G. Mills, and K. W. Jacobsen, Nudged elastic band method for finding minimum energy paths of transitions, in *Classical and Quantum Dynamics in Condensed Phase Simulations*, edited by B. J. Berne, G. Ciccotti, and D. F. Coker (World
Scientific, Singapore, 1998).

[40] G. Henkelman and H. Jónsson, Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points, *The Journal of Chemical Physics* **113**, 9978 (2000).

[41] M. T. Loponen, R. C. Dynes, V. Narayanamurti, and J. P. Garno, Measurements of the time-dependent specific heat of amorphous materials, *Physical Review B* **25**, 1161 (1982).

[42] R. Gutiérrez, S. Karmakar, Y. G. Pollack, and I. Procaccia, The static lengthscale characterizing the glass transition at lower temperatures, *EPL (Europhysics Letters)* **111**, 56009 (2015).

[43] W. Kob and H. C. Andersen, Testing mode-coupling theory for a supercooled binary Lennard-Jones mixture I: The van Hove correlation function, *Physical Review E* **51**, 4626 (1995).

[44] S. Sastry, Liquid Limits: Glass Transition and Liquid-Gas Spinodal Boundaries of Metastable Liquids, *Physical Review Letters* **85**, 590 (2000).

[45] V. Testard, L. Berthier, and W. Kob, Influence of the Glass Transition on the Liquid-Gas Spinodal Decomposition, *Physical Review Letters* **106**, 125702 (2011).

[46] S. Plimpton, Fast Parallel Algorithms for Short-Range Molecular Dynamics, *Journal of Computational Physics* **117**, 1 (1995).

[47] L. Berthier, E. Flenner, C. J. Fullerton, C. Scalliet, and M. Singh, Efficient swap algorithms for molecular dynamics simulations of equilibrium supercooled liquids, *Journal of Statistical Mechanics: Theory and Experiment* **2019**, 064004 (2019).

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

[49] Y. S. Elmatad, D. Chandler, and J. P. Garrahan, Corresponding states of structural glass formers. ii, *The Journal of Physical Chemistry B* **114**, 17113 (2010).

[50] D. Frenkel and B. Smit, *Understanding Molecular Simulation From Algorithms to Applications* (Academic Press, 1996).

[51] W. E. W. Ren, and E. Vanden-Eijnden, String method for the study of rare events, *Physical Review B* **66**, 052301 (2002).

[52] E. Bitzek, P. Koskinen, F. Gähler, M. Moseler, and P. Gumbsch, Structural Relaxation Made Simple, *Physical Review Letters* **97**, 170201 (2006).

[53] J. Guénolé, W. G. Nöhring, A. Vaid, F. Houllé, Z. Xie, A. Prakash, and E. Bitzek, Assessment and optimization of the fast inertial relaxation engine (fire) for energy minimization in atomistic simulations and its implementation in lammps, *Computational Materials Science* **175**, 109584 (2020).

[54] W. E. W. Ren, and E. Vanden-Eijnden, Simplified and improved string method for computing the minimum energy paths in barrier-crossing events, *The Journal of Chemical Physics* **126**, 164103 (2007).

[55] G. Henkelman, B. P. Uberuaga, and H. Jónsson, A climbing image nudged elastic band method for finding saddle points and minimum energy paths, *The Journal of Chemical Physics* **113**, 9901 (2000).

[56] C. Scalliet, L. Berthier, and F. Zamponi, Nature of excitations and defects in structural glasses, *Nature Communications* **10**, 5102 (2019).

[57] Q. Liao and L. Berthier, Hierarchical landscape of hard disk glasses, *Physical Review X* **9**, 011049 (2019).

[58] C. Artiaco, P. Baldan, and G. Parisi, Exploratory study of the glassy landscape near jamming, *Phys. Rev. E* **101**, 052605 (2020).

[59] G. H. Vineyard, Frequency factors and isotope effects in solid state rate processes, *Journal of Physics and Chemistry of Solids* **3**, 121 (1957).

[60] R. B. Lehoucq, D. C. Sorensen, and C. Yang, *ARPACK Users' Guide: Solution of Large-Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods* (Society for Industrial and Applied Mathematics, 1998).

[61] G. Frossati, J. leG Gilchrist, J. C. Lasjaunias, and W. Meyer, Spectrum of low-energy dipolar states in hydrated vitreous silica, *Journal of Physics C: Solid State Physics* **10**, L515 (1977).

[62] V. Jakšić and J. Segert, On the Landau–Zener formula for two-level systems, *Journal of Mathematical Physics* **34**, 2807 (1993).

[63] T. Vegge, J. P. Sethna, S.-A. Cheong, K. W. Jacobsen, C. R. Myers, and D. C. Ralph, Calculation of quantum tunneling for a spatially extended defect: The dislocation kink in copper has a low effective mass, *Phys. Rev. Lett.* **86**, 1546 (2001).

[64] C. Lévesque, S. Roorda, F. Schiettekatte, and N. Mousseau, Internal mechanical dissipation mechanisms in amorphous silicon, *Physical Review Materials* **6**, 123604 (2022).

[65] M. Ozawa, C. Scalliet, A. Ninarello, and L. Berthier, Does the adam-gibbs relation hold in simulated supercooled liquids?, *The Journal of Chemical Physics* **151**, 084504 (2019).

[66] C. Scalliet, B. Guiselin, and L. Berthier, Thirty Milliseconds in the Life of a Supercooled Liquid, *Physical Review X* **12**, 041028 (2022).

[67] T. A. Weber and F. H. Stillinger, Local order and structural transitions in amorphous metal-metalloid alloys, *Physical Review B* **31**, 1954 (1985).

[68] K. González-López, M. Shivam, Y. Zheng, M. P. Ciamarra, and E. Lerner, Mechanical disorder of sticky-sphere glasses. I. Effect of attractive interactions, *Physical Review E* **103**, 022605 (2021).

[69] L. Wang, A. Ninarello, P. Guan, L. Berthier, G. Szamel, and E. Flenner, Low-frequency vibrational modes of stable glasses, *Nature Communications* **10**, 26 (2019).

[70] G. Kapteijns, E. Bouchbinder, and E. Lerner, Universal Non-phononic Density of States in 2D, 3D, and 4D Glasses, *Physical Review Letters* **121**, 055501 (2018).

[71] L. Wang, L. Fu, and Y. Nie, Density of states below the first sound mode in 3D glasses, *The Journal of Chemical Physics* **157**, 074502 (2022).

[72] D. Khomenko, D. R. Reichman, and F. Zamponi, Relationship between two-level systems and quasilocalized normal modes in glasses, *Physical Review Materials* **5**, 055602 (2021).

[73] G. Folena and P. Urbani, Marginal stability of soft anharmonic mean field spin glasses, *Journal of Statistical Mechanics: Theory and Experiment* **2022**, 053301 (2022).

[74] L. Berthier, Self-induced heterogeneity in deeply supercooled liquids, *Phys. Rev. Lett.* **127**, 088002 (2021).

[75] M. Marchi and D. Chandler, Path-integral calculation of the tunnel splitting in aqueous ferrous–ferric electron transfer, *The Journal of Chemical Physics* **95**, 889 (1991).

[76] C. L. Vaillant, D. J. Wales, and S. C. Althorpe, Tunneling splittings from path-integral molecular dynamics using a Langevin thermostat, *The Journal of Chemical Physics* **148**, 234102 (2018).