{'fraction_non_alphanumeric': 0.034198561776586305, 'fraction_numerical': 0.013385983661492202, 'mean_word_length': 4.361561619889288, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Microscopic observation of two-level systems in a metallic glass model

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 twolevel 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 2 dependence of the thermal conductivity, and the plateau in ultrasonic sound attenuation [4]. Despite its success, alternative explanations for these experimental observations have been proposed [5]. In particular, the interactions between TLS [6], aspects related to the scattering of phonons by TLS [7][8][9][10], and the collective TLS dynamics [11] remain the subject of debate. Computational models can in principle be used to establish the relative merits and the ultimate validity of different hypotheses [12][13][14][15][16][17].
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 [18]. 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 [19,20].
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 f ), which is "the temperature at which the glass sample would find itself in equilibrium if suddenly brought there from a given state" [21] 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 [22], on vapor-deposited indomethacin [23], and in numerical simulations of polydisperse soft-spheres [16], where T f was decreased from above (hyperquenched glass) to below (ultrastable glass) the experimental glass transition temperature T g . However, the numerical results for continuously polydisperse systems [16] 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 [24,25] to computer models of metallic glasses [26], 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 f 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,27]. 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 de-crease in dissipation in amorphous solids at low temperatures and a concomitant improvement in the material attributes associated with a variety of practical applications.
Lastly, we revisit previous numerical protocols for TLS determination [16] and provide detailed insight on how the measured 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 [28]. 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.

A. General strategy
Our ultimate goal is to detect and characterize numerically TLS in a model metallic glass, and investigate how their number and properties evolve with glass preparation. In experiments, TLS are naturally excited at low temperature where quantum effects become relevant. However, simulating the quantum dynamics of glasses containing thousands of particles 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 connecting distinct potential energy minima) with associated quantum splitting smaller than E, counted per atom and per glass sample. The TLS model [1,2] postulates the small-E scaling The plateau value n 0 reached by n(E)/E at small E allows one to estimate n T LS , the number of active TLS per atom in a typical glass sample at temperature T Q where experiments are performed, typically around 1 K, The quantity n T LS 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 temperatures, while more detailed explanations and analysis will be presented in the following sections. There are four important 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 supercooled liquid phase at different temperatures T f ranging from the mode-coupling crossover temperature T mct down to the experimental glass transition temperature T g using the swap Monte Carlo algorithm [24]. We then form glasses by rapidly quenching the configurations to a lower "exploration" temperature T exp using conventional molecular dynamics (MD). Therefore, our glassy states at T exp would correspond to equilibrium supercooled liquid configurations if brought back to T f , which allows us to identify T f with Tool's fictive temperature [21]. Varying T f is similar to varying the cooling or deposition rate in experiments. Each of the N g swap-generated equilibrium configurations defines a "glass sample" or a "glass metabasin." 2. Explore using MD the potential energy landscape of each glass metabasin [29,30], which contains multiple potential energy minima, or inherent structures (IS) [31][32][33], 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 minimumenergy path (MEP) connecting these pairs of IS using the nudged elastic band (NEB) method [34,35]. We obtain a library of N DW double-well potentials.
4. The minimum energy path provides an effective onedimensional 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 [12,16]. We call N DW (E) the total number of double wells with a quantum splitting lower than E, with N DW (E → ∞) = N DW . The quantity N DW (E ∼ 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 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 Q over an observation time t w , the TLS need not only to have the relevant energy splitting E ∼ k B T Q , but also to have a tunneling rate Γ such that Γt w 1. Because the TLS model postulates a flat distribution of log Γ at small Γ, the effective number of active TLS grows as logt w at large t w [1,2,36]. In summary, according to the TLS model we should expect a finite n(E) ∝ k B T Q logt 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 → ∞. The TLS model then predicts that the remaining ratio N DW /(N g × 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 ∝ N g and the ratio N IS /N g then converges to a finite value when N g → ∞ 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 T LS 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 r i and r j , separated by a distance r i j = |r i − r j | interact with a Lennard-Jones (LJ) pairwise potential V (r i j ), given by if the particle pair distance r i j /σ i j = 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 trunca-tion, we consider the smoothing polynomial [37] S(r i j ) = 4ε i j C 0 +C 2 r i j The choice of coefficients: C 0 = 10/x 6 c − 28/x 12 c , C 2 = 48/x 14 c − 15/x 8 c and C 4 = 6/x 10 c − 21/x 16 c ensures the continuity of the potential V (r i j ) 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 σ AA and the unit energy is ε AA , noted σ and ε for simplicity. The interaction parameters are ε AB = 1.5ε, ε AC = 0.9ε, ε BB = 0.5ε, ε BC = 0.84ε and ε CC = 0.94ε for the energies; and σ AB = 0.8σ , σ AC = 1.25σ , σ BB = 0.88σ , σ BC = 1.0σ and σ CC = 0.75σ for the ranges. All particles have the same mass m. The unit time is τ LJ = mσ 2 /ε. This ternary mixture is an extension of the well-known Kob-Andersen model (KA) [38] 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. [26] where the model is referred to as KA 2 . We simulate N particles at number density ρ = 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 ρ = 1.2 [38]. Indeed, the pairwise attraction gives rise to a liquid-gas spinodal at low enough temperature, which may intersect the glass transition line [39], leading to a gas-glass instability [40]. The choice ρ = 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.

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 modecoupling 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 [41]. A time-step of dt = 5 × 10 −3 τ 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 [42]. Based on previous work [26], we have found that short blocks of 10 MD steps with dt = 5 × 10 −3 τ LJ interspersed with blocks of 2100 (1.75 × 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 thermodynamics of the system, provided proper decorrelation can be achieved within the simulation time window [24].
In order to monitor structural relaxation during standard and swap MC dynamics, we have estimated the structural re-FIG. 1. Angell plot of relaxation time τ α versus inverse temperature measured in physical molecular dynamics (MD, green) and equilibration particle-swap dynamics (Swap, magenta). The onset T o 1.0 (grey), and mode-coupling crossover T mct 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". laxation time τ α from the self-intermediate scattering function F s (q = 7.2,t = τ α ) = 1/e. The relaxation time τ α 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 ∼ 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, τ α ∝ (T − T mct ) −γ , in a region of moderate supercooling, which results in an empirical estimate of T mct ∼ 0.62. Finally, the laboratory glass transition temperature T g is the temperature at which τ α (T g ) = 10 12 τ 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 [43] and the parabolic law [44]. These extrapolations define a range for the laboratory glass transition temperature of T g ∼ 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×10 −3 τ 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 [45]. 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 τ 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 τ total = N cg τ 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 α and IS β in a matrix T α→β , where the labeling is ordered by energy (see Fig. 2). If a transition occurs between α and β in both directions at least once (T α→β ≥ 1 and T β →α ≥ 1), we consider that the pair (α, β ) 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. [46], which is implemented in LAMMPS as a modification of the nudged-elastic band (NEB) approach [35]. 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 κ = 0.125 εσ −2 . The optimization of the string method was done using the FIRE minimizer [47,48], while interpolation of the MEP was carried out with a cubic spline. This should result in an error scaling as O(N −4 images ) [49]. The location of the transition state is identified using a climbing image approach [50] 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 τ cg and a relatively small number of independent glass 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 τ 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. Dataset 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 = τ total .

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 [32,33]. 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 (∼ 10 5 steps) during which the system is quenched to a lower temperature We can assign to each IS α its energy and compute the Euclidean distance between a pair (α, β ) of IS as where r (α) i is the position of the ith particle in the α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 α observed during exploration. The IS in the matrices are ordered by increasing energy e α .
The most striking aspects that emerge from these figures are (i) the large heterogeneity of the number of IS (we find between n IS ∼ 6 and n IS ∼ 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 αβ matrix [51][52][53], (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 τ cg = dt, i.e. minimizing the potentialg energy at every MD step. The number of distinct IS per glass is counted, and a distribution of persistence times (τ p ) is inferred from the resulting time series: here τ 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 τ p , for different T f , is shown in Fig. 3a. The shape of the CDF indicates that the distribution of τ p is bimodal, with a first peak around τ p ∼ dt (i.e. a single MD step) and a second one around τ p ∼ 10 2 dt. A more careful analysis of the IS time series reveals that the peak at short τ 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 results.
We find a sub-linear power-law dependence of the number of minima on the exploration time, i.e.
with an exponent β ∈ [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 [51]. 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 τ cg /dt MD steps only, instead of after every step. We have chosen τ cg in our final data production runs for landscape exploration as a compromise between accuracy and efficiency. A low value of τ 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 ∼ 100 − 1000 steps. The values τ cg 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 ∼ τ LJ . The time unit τ LJ seems to be close to the characteristic time for moving outside of the basin of a single IS.

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 α→β ≥ 1 and T β →α ≥ 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) ∝ logt ∝ 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 im- plies 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 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 ∼ 1 K in experiments and can be obtained from comparing the interparticle distance with the thermal wavelength, leading to: 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 ∆V , barrier height V b = V a − ∆V /2 (black arrows), and the curvature of the potential-energy surface at the lowest minimum ω 1 , highest minimum ω 2 , and transition state ω T S (blue curved arrows). The quantum splitting E is the difference between the first two quantum energy levels ε 1 and ε 2 (dashed red).
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, ∆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 − 1 2 ∆V . As a first approximation, we reduce the tunneling problem to the effective one dimensional potential of the MEP. Considering a normalized reaction coordinate ξ ∈ [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 An interpolated MEP obtained from the NEB calculation is discretized by dξ = d/N ξ , where N ξ = 2000, using a cubic spline. The Laplacian is evaluated with a 5-point finite difference stencil. We useh = 1 to evaluate the quantum properties of the DW. This leads to a dimensionless effective mass parameterm which controls the "quantumness" of the problem, as first introduced by Vineyard [54], Them 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 [16]. In this case, for the two sets of units detailed in Sec. IV B, namely argon (Ar) and nickel-phosphorous (NiP), the value ofm is 1200 and 5250, respectively. The corresponding values of k B T Q are 0.005ε (for Ar) and 0.0012ε (for NiP).
The five smallest eigenvalues ε 1 , ..., ε 5 are then evaluated using ARPACK [55]. The splitting is given by the first two energy levels, E = ε 2 − ε 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 in the localized representation. The diagonal splitting of the TLS can be estimated by where ω 1 and ω 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 The tunneling matrix element ∆ 0 can be obtained from the effective one-dimensional potential via the WKB approximation as 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 ε 2 for the calculation of Γ, andε for ∆ 0 .
As such, the quantities ∆ and ∆ 0 are expected to approximately obey the relationship E ≈ ∆ 2 + ∆ 2 0 . The quantum splitting E is plotted against ∆ 2 + ∆ 2 0 in Fig. 6a, showing the expected correlation at all temperatures. Deviations from E ≈ ∆ 2 + ∆ 2 0 are observed at low splittings where discretization effects may become significant. Fig. 6b shows the cumulative probability distribution function F(∆ 0 ) of the tunneling matrix element ∆ 0 . The TLS model predicts a flat PDF for log ∆ 0 , which corresponds to F(∆ 0 ) ∼ log ∆ 0 at small ∆ 0 . Our data instead suggest that F(∆ 0 ) ∼ ∆ ϕ 0 with a very small exponent ϕ ≈ 0.1, corresponding to a PDF p(∆ 0 ) ∼ ∆ ϕ−1 0 . Such behavior was discussed in more detail elsewhere [2,56].
The decay rate Γ is also expected to approximately obey Γ ∝ ∆ 2 0 due to Landau-Zener tunneling [57], which is seen in the inset of Fig. 6b. Finally, the CDFs of ∆ and E are shown in  The ratio of tunneling TLS to generic double-wells, N DW (T Q )/N DW , seems to be roughly constant ∼ 0.3%, at all temperatures. Because the function N DW (E)/N DW converges to a stable finite limit when N DW → ∞, 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.

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 DW potentials in relation to their estimated quantum splitting E and the preparation temperature T f of the glass. The DW asymmetry, ∆V , is shown in Fig. 7a where one can notice a threshold value ∆V ∼ 10 −2 , below which the distribution of quantum splittings becomes much broader.
The statistics of the energy barrier, V b , 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 V b . DW sampled in lower T f glasses typically have a higher barrier than those found in glasses prepared at high T f .
The statistics of the tunneling decay rate Γ (∼ ∆ 2 0 ) are shown in Fig. 7c. The distribution of tunneling rates becomes narrower around higher values as the tunnel splitting decreases. As expected, low-splitting tunneling DW tend to have a high tunneling rate. Still, some of the low-E DW have small tunneling rate and will thus tunnel too slowly to be readily observed on a typical experimental timescale. 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 i αβ = |r defined such that 1 ≤ PR ≤ 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 T f . 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 T f T mct . 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. [58] The patterns of real-space displacements between two IS forming a DW potential are illustrated in Fig. 8. Atoms are shown with a radius proportional to their interaction range σ and colored by d 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. We highlight the first PR (to the nearest integer) particles that move the most. In general, displacements in our database of DW are highly localized and involve only several atoms. The most localized DW potentials have a PR ∼ 1 such as the one shown in Fig. 8a. A typical TLS is shown in Fig. 8b, it has a PR ∼ 7 close to the average of 8. 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 T mct 0.62 glasses.

B. Depletion of tunneling two-level systems with increasing glass stability
The standard TLS model predicts a plateau n(E)/E → n 0 as E → 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 ∼ 5 × 10 −2 . In addition to the low-E plateau behavior observed at each 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 f = 0.617. Data for two sets of units: NiP (circle,m = 5250), and Ar (square, m = 1200). PSS data from Ref. [16]. To compare the models we have fixed the exploration time t = 10 3 τ LJ to estimate n 0 . 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 polydisperse soft-spheres (PSS) [16].
In order to compare the TLJ data with the PSS model we have plotted n 0 /m versus the stability of the glass, as encoded in the relaxation time ratio log(τ α (T f )/τ o ), in Fig. 9b. In this representation, glass stability increases from left to right. The relaxation time τ α (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 [59]. The n 0 values for PSS are from Ref. [16], while we have refined our τ α -extrapolation using the latest data obtained via long MD simulations of PSS [60]. 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 τ 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. [16] 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 ofm to probe quasi-universality by emulating different material properties within the same model. This confirms that the precise value ofm does not affect our main conclusion.
The choices ofm 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 [13], and has the following parameters: σ = 3.405 × 10 −10 m, ε/k B = 125.2 K and a mass of m = 6.634 × 10 −26 kg. The unit of time is τ LJ =2.1 ps. In these units, the numerically estimated glass transition temperature is T g ≈ 59 K, T Q ≈ 0.658 K andm ≈ 1200. For the preparation temperatures T f = 0.488, 0.509, 0.558 and 0.617 one gets n sim 0 ∼ 0.04, 0.09, 0.14 and 0.674 ε −1 σ −3 which corresponds to n exp 0 ∼ 4.05 × 10 47 , 1.19 × 10 48 , 2.61 × 10 48 and 8.95 × 10 48 J m −3 .
We find a modest variation in the absolute value of the tunneling defect density withm. 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 ofm 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 [62]. 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, Once the moduli have been determined, we calculate the transverse c t = G/(mρ) and longitudinal c l = (K + 4G 3 )/(mρ) sound velocities, and from them the frequency ω t of the first phonon, and the Debye frequency, 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 ω t , when the model size is small enough to avoid significant hybridization with phonons. It is typically assumed that the density of states D(ω) scales as ω 4 in the low-frequency limit [64].
In Fig. 10a we show the integrated density of states I(ω), scaled by the expected ω 5 for the larger 12000-atom mod-els at each T f . Provided D(ω) ∝ A 4 ω 4 , we would expect a plateau in I(ω)/ω 5 just below ω 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 [65]. The low-frequency behaviour of D(ω) 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(ω) 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 ω 5 D 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 τ 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 [66]. In particular, our data seems to preclude a direct proportionality between the two quantities of the kind predicted in Ref. [28]. 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 [16]. The number of two-level systems converges relatively quickly at ∼ 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 twolevel 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 ≈ 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. [51][52][53]. The observation of these extremely strong sample-to-sample fluctuations in glass samples deserves further investigation [67,68].
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 tradeoffs might come in to play with other types of defects such as quasi-localized modes.
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 logt 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 ω 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 T LS of TLS even stronger, some of the approximations used to estimate n T LS should be removed from the simulations. The most ambitious task would be to obtain the quantum splitting for candidate TLS from a full quantummechanical calculation rather than a 1D projection along the reaction path. A second task would be to improve the computational tools used in this work to significantly increase the library of TLS found numerically. 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 understand and quantify this evolution and understand how it impacts the physical properties of glassy and viscous liquid states.

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/xx.xxxx/xxxxx.yyyyyy.
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 onedimensional local reaction coordinate x i , with an associated potential taken, for the sake of illustration, to be of the form If κ 3 i + 9h 2 i /8 < 0, this potential has two local minima, leading to a DW with energy splitting ∆V i (κ i , h i ). Let us introduce a spin variable σ 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 Such a system has N IS = 2 K local minima, corresponding to all possible combinations of the σ i , with energy E IS ({σ 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 (σ i = 0, σ j = 0) → (σ i = 1, σ j = 1) would be decomposed into elementary events, e.g. (σ i = 0, σ j = 0) → (σ i = 0, σ j = 1) → (σ i = 1, σ j = 1), leading to the appearance of an intermediate local minimum in (σ i = 0, σ j = 1) in the MEP. This very simple argument thus suggests a relationship of the kind N DW ≈ log N IS .
In summary, for K ∝ N independent excitations in each glass, and for N g glasses, we would have in total N IS ∼ N g exp(K) ∼ N g exp(N), N DW ∼ N g K ∼ N g N, and N DW /(N N 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. [67] and references therein) For weak enough coupling J i j , 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 ∼ N 2 IS . Depending on the nature of elementary excitations and their interactions, one would then observe log N IS N DW N 2 IS , 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 softsphere glasses studied in [16]. 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.
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 powerlaw. 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.