Identification of structural motifs as tunneling two-level systems in amorphous alumina at low temperatures

One of the most accepted models that describe the anomalous thermal behavior of amorphous materials at temperatures below 1 K relies on the quantum mechanical tunneling of atoms between two nearly equivalent potential energy wells forming a two-level system (TLS). Indirect evidence for TLSs is widely available. However, the atomistic structure of these TLSs remains an unsolved topic in the physics of amorphous materials. Here, using classical molecular dynamics, we found several hitherto unknown bistable structural motifs that may be key to understanding the anomalous thermal properties of amorphous alumina at low temperatures. We show through free energy profiles that the complex potential energy surface can be reduced to canonical TLSs. The predicted tunnel splittings from instanton theory, the number density, dipole moment, and coupling to external strain of the discovered motifs are consistent with experiments.


I. INTRODUCTION
Since the first measurements of the thermal behavior of glasses at temperatures below 1 K, it became clear that such systems possess very peculiar properties that unite almost all disordered materials regardless their chemical composition and bonding 1-4 and distinguish them from their crystalline counterparts. Namely, the heat capacity and thermal conductivity of amorphous systems increase with temperature almost linearly and quadratically, respectively, which are in sharp contrast to the Debye T 3 behavior of crystals 5 . This anomalous behavior was attributed to microscopic tunneling twolevel systems (TLSs) that are sparsely present in disordered materials 6,7 .
The same TLSs are believed to be the main source of decoherence and noise in Josephson junction-based superconducting qubits 8,9 , various nanomechanical 10,11 and optical 12 resonators using amorphous materials. Exacerbating the problem is that each sample has its own set of TLS frequencies that cannot be controlled at will 8,[13][14][15][16] . As often happens, the phenomenon that is parasite for many applications at the same time opens the way to others. Such applications, as TLS-based qubits 17,18 , devices for non-linear optical measurements 19 or quantum memory 20 based on TLSs, have been proposed. However, inability to control TLSs and lack of understanding of their atomistic structure hinders these developments as well.
The recent breakthrough experiments by Grabovskij et al. 21 on microwave response of individual TLSs under strain leave no doubt that at least in the Josephson junctions with the amorphous alumina (Al 2 O 3 ) barrier, the low energy excitations are described by the standard TLS model. Nevertheless, no consensus has been reached yet on the nature of the tunneling particles. The fact they should carry a substantial electric charge that couples to the Josephson junction electric field has generated interpretations based on tunneling of electrons [22][23][24][25] . To fit into the experimental range of small energy splittings, addi-tional arguments had to be invoked such as polaronic dressing of tunneling electrons 22 . TLS models based on the delocalization of a single oxygen atom 26 or hydrogen impurities 27,28 have been also proposed. We think, however, that electronbased models may not be needed because atoms in an ionic material such as amorphous alumina already carry significant partial charges. 29 Neither is there the need to assume the presence of any specific impurities. The disordered alumina structure should already contain intrinsic bistable motifs. Here, we undertake the computational search and characterization for such structural motifs.
The paper is organized as follows. In Sec. II, we review the standard TLS model. The computational methods and details are described in Sec. III In Sec. IV, we present and discuss the structural motifs found in our computational simulations on amorphous alumina at low temperatures. Finally, the conclusions are drawn in Sec. V.

II. THE TLS MODEL
In 1972, Phillips 6 and Anderson et al. 7 proposed independently that in amorphous materials some atoms could tunnel quantum-mechanically between two nearly isoenergetic configurations forming TLS, just like the umbrella motion of the ammonia molecule. They assumed that this scenario could be modelled by an effective double well potential whose barrier height V 0 and energy difference between the two minima ∆ were broadly distributed across the amorphous sample. The energy splitting between the two lowest vibrational eigenstates should be small (∼ 0.1 meV) to contribute to the properties of glasses at temperatures below 1 K 1,2, 8,9,[13][14][15][16]21 .
This tunneling splitting is given by δE = ∆ 2 + ∆ 2 0 , where the coupling energy ∆ 0 = Ω exp −d √ 2m eff V 0 / is proportional to the overlap between the unperturbed localized wavefunctions in each well 1,2,4 . Here, d is the well separation, Ω arXiv:1411.0529v2 [cond-mat.mtrl-sci] 16 Nov 2014 is a typical well vibration frequency, and m eff is the effective mass of the tunneling particle.
The energy asymmetry between the two wells of a TLS can be also varied using an external strain field as δ∆ = 2γ , where is the dimensionless strain and γ is the strain-asymmetry coupling that was measured to be about 1 eV 1,2,21,30 . The coupling between the TLSs and an external electric field F is of the dipole type δ∆ = 2p · F 8,9,[13][14][15][16]21 with the TLS dipole moment p corresponding to a single effective charge moving by a distance of a single atomic bond 9,31 .
These TLSs are sparsely present in the bulk of amorphous materials (∼ 10 21 eV −1 cm −3 according to experiments 1,2,21 ) so measurements on those samples yield averages over many different TLSs. This explains why the experimental verification of the TLS model remains extremely difficult and, though supported by numerous indirect data 1,2,4 , is still hotly debated 32,33 . Modern extensions of the original TLS model, including the "soft-mode" model 34 and the "defect-interaction" model 35,36 , were recently reviewed 4 .

III. COMPUTATIONAL METHODS AND DETAILS
From the computational viewpoint, finding TLSs in amorphous materials is difficult due to the required large system sizes and long time scales. Firstly, the low density of TLS imposes severe limitations on the sizes of the systems that have to be involved in search for TLS. Secondly, the necessity of descending to very low temperatures to freeze all structural transformations except for flips between the configurations forming TLS makes these processes hardly activated. As a result, TLSs appear as extremely "rare events" both in space and time, which virtually excludes First-Principles methods from the available theoretical tools. For this reason, accurate and transferable force fields have been instead adopted in most computer simulations in this research field. To add further complication, the quenching simulated by computers is typically several orders of magnitude faster than the natural cooling process 37 . Moreover, there is no way to guarantee that TLSs are present after the melting-quenching computational protocol used to generate the amorphous configurations. If present, the TLSs are metastable entities that are easily destroyed by a variety of physical factors such as temperature and/or mechanical strain. A final caveat for computer simulations is that one should use a computational cell with a number of atoms much larger than the atoms comprising the TLS, which is often unknown a priori. In summary, these limitations in time scales and system sizes -and hence the need of much sampling-render the use of ab initio techniques prohibitive.
Nowadays, classical molecular dynamics (MD) represents the best compromise between computational cost and accuracy in the quest for TLS candidates. For instance, in amorphous silica (SiO 2 ), TLSs were reported to involve large cooperative re-orientations or coupled rigid rotation of nearly 30 SiO 4 tetrahedra with an estimated energy barrier of 0.06 eV [38][39][40] . The MD approach adopted in this work entails the following three consecutive steps: (1) Generating amor- where q Al = 1.244690|e| and q O = −0.829793|e| are the partial charges for aluminum and oxygen atoms, respectively. (2) Searching for any large amplitude fluctuation in atomic positions over time during MD simulations.
(3) Finally, for each large hop, we reconstruct the free energy profile to assess whether it qualifies as a TLS candidate. We calculate the tunnel splittings on the potential energy profile obtained at 0 K for the successful TLS candidates.
In all calculations, we used two force fields, hereafter called, I 42 and II 43 , to check the reproducibility and robustness of our motifs. The first of these force-fields (I) has been extensively tested against experimental properties of crystalline [44][45][46][47] , liquid 48,49 and amorphous alumina 41 . The second force-field (II) was fitted to ab initio data on crystalline and liquid alumina 43 . The structural properties of our alumina systems have been compared to previous calculations (see Fig. 5) and excellent agreement is found. The definitions of force fields I and II are presented in Table I and Table II, respectively.
Two cell sizes were considered for the MD simulations of alumina: a medium-sized cell containing 1500 atoms (300 Al 2 O 3 units) and a smaller one with 360 atoms (72 Al 2 O 3 units). Different MD codes and simulation parameters were used for each system size. In all cases, we cross-checked our findings with the different simulation protocols employed.
All samples of amorphous alumina were generated following the standard quench-and-melt protocol described in Ref. 41 . The liquid was obtained from the melting of the experimental crystalline structure of α-alumina (30 atoms per hexagonal cell with parameters a = b = 4.76020 Å, c = 12.9933 Å, α = β = 90 • , γ = 120 • and density of 3.98 g cm −3 ). Liquid alumina was equilibrated at high temperature 5000 K and low density 2.75 g cm −3 during 200 ps. The system was then cooled down to 3000 K at a cooling rate of 20 K ps −1 . The resulting configuration was subsequently equilibrated at 3000 K for over 100 ps. Finally, the systems were further quenched to low temperature (25 and 5 K) at a cooling rate of 4 K ps −1 . Our radial distribution functions (see Fig. 5) agree very well with the previous MD simulations 41,43 and demonstrate that our structures are truly amorphous.
The amorphous system with 1500 atoms equilibrated to a final cell volume of 24.1335 × 24.1335 × 26.3496 Å 3 , yielding a density of 3.24 g cm −3 . The electrostatics was handled with the standard Ewald summation technique. A cutoff radius of 11 Å was used to truncate the short-ranged interactions. During equilibration, we controlled the temperature and pressure using the Nosé-Hoover implementations of the thermostats 50 and barostat 51 , respectively. A time step of 0.4 fs was used to ensure a stable integration of the equations of motion. These calculations were performed using the freely available codes CP2K 52 and LAMMPS 53 .
The molecular dynamics simulations for the 360-atoms system were performed with the MD-kMC code 54 . The long-range Coulomb interactions were also calculated with the standard Ewald summation technique. The system was equilibrated to a final cubic cell of volume 15.61777 3 Å 3 after isotropic compression to reach a density of 3.20 g cm −3 . A cutoff of 7.81 Å was used to truncate the short-ranged interactions. The time step was 1 fs. The Berendsen thermostat was employed to control the temperature 55 .
To study the dynamics of amorphous systems at low temperatures and find possible two-level systems (TLSs), we conducted long microcanonical (constant energy) MD simulations on the equilibrated amorphous configurations for over 1 ns. We carefully monitored the time evolution of the atomic fluctuations about their average positions. We selected those atoms that deviate from their average positions more than 0.15 Å during time intervals longer than 1 ps. To search for such fluctuations, in practice we ran up to a hundred trajectories, which is statistically more efficient than performing long simulations on larger systems. While at temperature above 50 K many of such fluctuations happen simultaneously and render the analysis complicated, at temperatures below 10 K they are too rare and not accessible by standard MD simulations. In practice, we chose an intermediate temperature (∼ 25 K) to search for bistable motifs. After visual inspection of the MD trajectories, a collective variable δ (typically an asymmetric stretch) was chosen to describe the observed transformations. All visualization and snapshots were done with the VMD code 56 .
Free energy profiles were computed for the found "rare events" to check whether they qualify as TLSs, that is, to confirm their bistability. To this end, we used and compared a variety of well-established techniques, including unbiased long MD runs, the metadynamics method of Laio and Parrinello 57 , the average biasing force method 58 . These last two biased MD simulations improve drastically the sampling of rare events. We also estimated the minimum potential energy path between the wells at 0 K by running geometry optimizations with an added harmonic restraint potential V (δ) = k 2 (δ − δ 0 ) 2 on the collective variable δ. We effect the mapping of the potential energy profile by moving the center δ 0 of the parabola. We found consistent energy profiles using the different sampling techniques. In Fig. 6, we provide such comparison where we include two additional sampling methods: the blue-moon ensemble technique 59 and the nudge elastic band method 60 . Fig. 1(a) shows some of the largest fluctuations in atomic positions found during a microcanonical (NVE) MD simulation of amorphous alumina at 5 K using 1500 atoms and force field I 42 . Using the asymmetric stretch δ = | r O2−Al | − | r Al−O1 | as collective variable, we observe that for over 0.5 ns it displays the typical signature of a bistable "rare event" at that temperature, which corresponds to the transfer of a central Al  (Fig. 1). (a) The effect of fixing the position of atoms beyond a given radius (in Å) around the central Al atom on the minimum energy path between the minima ("fixed", in legend). The effect of changing the cutoff distance for short-ranged interactions is also shown ("full", in legend) where no frozen atoms are imposed. The 1500-atoms system and force field I were used. (b) Tuning the asymmetry of the energy profile at 0 K by mechanical strain (given as %) along X,Y, and Z axis. Results according to force field II 43 . (c) Estimation of the tunneling splitting by solving the Schrödinger equation numerically using an effective mass m eff = 15 a.m.u. and a fit of the profile at 0 K given in Fig. 1. Both eigenfunctions |ψ 1 | 2 and |ψ 2 | 2 were scaled down by a factor of 10 to fit in the graph. Their associated eigenvalues are indicated by blue and red horizontal lines, respectively. atom along the line between two axial O atoms ( Fig. 1(b)). At the transition state, the first coordination shell of the motif resembles a distorted octahedron with six O atoms around a central Al atom (typical coordination of corundum crystal) with elongated axial Al-O bonds, that are perpendicular to the equatorial plane ( Fig. 1(b)). 61 A cavity defect and low temperatures are essential for this motif survival. The same motif was also found with the force field II from the quenching of completely different melted alumina systems using different simulation parameters. This confirms that the structural motif in Fig. 1(b) is completely reproducible. Fig. 1(c) shows the reconstructed free energy profiles at 0, 5, and 25 K computed using both force fields I 42 and II 43 . The free energy difference between wells is about 2 meV in both cases and the largest energy barriers for the transition between minima are 4 and 12 meV for force fields I and II, respectively. Fig. 6 provides a further comparison of free energy profiles computed using several different sampling techniques and good agreement is found. All our energy profiles are asymmetric because such configurations are entropically favored due to the inherent disorder of glasses. Also, we do not have enough statistics to see these even "rarer" symmetric TLSs as they would demand extremely large system sizes and/or long simulation times. Regarding the effect of temperature, it is well-known that the TLS only survives at low enough temperatures (< 100 K). Fig. 1(c) shows that raising the temperature increases the asymmetry of the double-well profile. We note, however, that the effect is more evident for the force field II than for the I. There is a non-trivial influence of the atoms surrounding the octahedron that will be discussed below. Fig. 1(b) may give the wrong impression that the motif is rather localized and independent of its surroundings, however further tests show that this is not the case. Fig. 2(a) shows that the environment around the motif has a non-trivial influence on its stability. The motif involves tens of atoms rather than just one aluminum atom jumping between two axial oxygen atoms. Here, we froze all atoms beyond a certain radius about the central aluminum atom and recomputed the energy profile ("fixed", in legend) at 0 K along the collective variable δ. We see how the double-well profile deteriorates as we gradually shrink the sphere around central aluminum atom (that is, as more atoms are frozen). We estimate the characteristic radius of the bistable motif to be ∼ 8 Å. In Fig. 2(a), we also investigated the impact of changing the cutoff distance in the short-range interactions ("full", in legend) on the energy profiles and found little effect provided that a cutoff distance greater than 7 Å is used. We further explored the finite-size effects on the energy profiles by re-optimizing a larger cell containing the TLS motif in Fig. 1(b). A 2 × 2 × 2 replication of our 1500-atoms system (12000 atoms in total) yields a similar energy profile at 0 K (see Fig. 7) as the one shown in Fig. 1(c). This result confirms that our motifs are independent beyond a certain radius and much less affected by the farthermost surroundings, in agreement with one of the basic tenets of the TLS model.

IV. RESULTS AND DISCUSSION
As mentioned, the energy profile of the identified motif is typically asymmetric due to the intrinsic disorder of the amorphous state, but it can be symmetrized by applying strain to the alumina sample. Fig. 2(b) demonstrates that the energy profile obtained by force field II on the 360 atoms system (Fig. 1(c)) becomes symmetric upon stretching by 0.1-2.5 % depending on the strain direction. The coupling coefficient between strain and energy asymmetry is found to be γ ∼ 0.2-1.0 eV for different strain directions, in good agreement with experiments 1,2,21,30 .
Another property of TLSs that can be measured is its direct coupling to external fields. Using the partial charges of force field II, we estimated the change in electric dipole moment δp for motif in Fig. 1 from where i and f refer to our initial and final local minimum configurations, respectively, in the energy profile of Fig. 1(c). The transition between the two energy minima is associated with a change in electric dipole moment of δp = 4.2 D. This value corresponds to a single electron charge moving by 0.9 Å, in good agreement with experiments 9,31 . Nevertheless, the most important feature that determines whether a proposed structural motif can influence the lowtemperature thermal properties of the amorphous material and its microwave spectroscopy is the tunnel splitting. To this end, we estimated the effective mass of the collective variable for the 0 K profile of Fig. 1(c) using Eq. (E2). The calculated effective mass using both force fields I and II varies from 8 to 16 a.m.u. (Fig. 8) along the minimum energy path. It is noted that the reduced mass for the asymmetric stretch of a system of 3 collinear atoms (O-Al-O) is m O m Al / (2m O + m Al ) = 7.3 a.m.u.. This is another manifestation that not only one atom participates in the transition between the two energy minima. Using the full coordinate dependence of the effective mass (Fig. 8) we estimated a tunnel splitting of δE = 0.3−0.7 meV using instanton theory. Taking an effective mass equal to m eff ∼ 15 a.m.u., we also solved the one-dimensional Schrödinger equation numerically for a symmetrized version of potential energy profile in Fig. 1(c). The estimated value of 0.65 meV (Fig. 2(c)) is somewhat higher than the experimental data 1,2, 8,9,[13][14][15][16]18,21 . However, smaller values of the tunnel splitting can be achieved for the motifs with higher barriers and larger well separations. We note that reproducing such small experimental splittings is beyond the precision of current ordinary ab initio calculations.
After inspecting Fig. 1(b), it is reasonable to wonder if a similar motif is feasible but with only three oxygen atoms in the equatorial plane. Fig. 3(b) shows that this is indeed possible. This new motif has a trigonal bipyramidal shape and was found using the force field I 42 . Fig. 3(a) shows the typical time evolution of the relevant collective variables and displays a clear bistable behavior. Fig. 3(c) features a more symmetric profile and a higher energy barrier (8 meV) than its counterpart in Fig. 1 (4

meV).
Are these all the possible structural TLS motifs in amor- phous alumina? We think that the answer to this question is negative. There exists an ensemble of different structures responsible for the low-temperature behavior and that we just found a few of them. For instance, in Fig. 4 we show another plausible TLS motif found for the 360-atoms system using the force field II 43 . This motif is characterized by a central tetracoordinated oxygen atom (typical coordination in the crystal) jumping back and forth between two aluminum atoms, labelled as Al1 and Al2 in the Fig. 4 (Top). Although a mobile oxygen atom was also proposed by DuBois and coworkers 26 , our motif geometry is very different to theirs. We remark that their motif was an ad-hoc model and might not be stable in the real amorphous environment, while ours is the direct result of MD simulations that consider explicitly the surrounding lattice. Although we have not found their motif in our calculations, it would be interesting to explore its existence in further atomistic MD simulations of amorphous alumina. In Fig. 4 (Bottom) we show the corresponding free energy profile at 0 and 25 K along the collective variable The plot features profiles with an energy barrier of the order of 15 meV. Interestingly, the temperature seems to have a symmetrizing effect on the profiles. The relative abundance of found TLS motifs in our amorphous alumina systems was rather small, of the order of 1 motif per 10000 atoms. This value was a rough estimate from the (limited) statistics of our MD simulations, that is, the number of TLS found (7) per the total number of atoms simulated (>50000). Assuming a uniform energy distribution (a basic assumption in the standard TLS model) we estimate the density of experimentally relevant TLSs (those with energy asymmetry below 0.1 meV) as 1 TLS per 200000 atoms, approximately. A more accurate value of the TLS density could be calculated with the systematic search algorithm proposed by Reinisch and Heuer 62,63 that improves upon standard methods that may miss several TLS candidates.

V. CONCLUSIONS
In conclusion, using extensive computer simulations, we show that bistable structures are naturally present in amorphous alumina at temperatures below 100 K. We found several structural motifs that form TLSs in alumina. We also found that the most mobile atoms of the TLSs could be either oxygen or aluminum atoms. We crudely estimate the density of experimentally relevant TLSs (those with energy asymmetry below 0.1 meV) to be 1 TLS per 200000 atoms. This value was estimated from the number of TLS found (7) per the total number of atoms simulated (> 50000), and from the standard assumption of the TLS model that the energy distribution of TLSs is uniform. In particular, we identified a motif resembling a distorted octahedron where the transferring aluminum atom performs an umbrella-like motion between the two axial oxygens. The robustness of this motif was confirmed by using two different force fields, different MD simulation parameters, and by cooling different alumina melts. The properties of the corresponding TLSs are consistent with experimental observations in Refs. 21 . Combination of our results with the microwave spectroscopy of TLSs 21 opens a venue for the eventual understanding, and, possibly, utilizing TLSs in amorphous alumina. Our results suggest that similar motifs may exist in other amorphous materials with a local corundum structure (see Fig. 9). The extension to such systems as well as the effect of pressure on TLSs are interesting problems for the future.

Appendix A: Pair distribution functions and static structure factors
We have validated our liquid and amorphous alumina configurations by comparing the calculated pair distribution functions and partial static structure factors to previous calculations 41,48 . We provide such comparison in Fig. 5(a) and Fig. 5(b), where we show the radial distributions functions for liquid alumina at 2200 K and its corresponding amorphous phase at 650 K, respectively. These temperatures were chosen to facilitate the comparison with previous computational data 41,48 . We report our results for the 1500-atoms system with force field I and for the 360-atoms system with force field II. Fig. 5(c) also shows a comparison between our computed partial static structure factors for the different pairs with previous calculations using force field I 48 . Our calculations using force field I are in excellent agreement with previous similar calculations, which validates our alumina configurations.

Appendix B: Comparison of different free energy sampling techniques
Different sampling methods used to compute the energy profile in Fig. 1(c) are compared in Fig. 6. The alumina system investigated was the one with 1500 atoms described by the force field I. Good agreement is found for all methods.

Appendix C: Finite-size effects on energy profiles
We compared the energy profile at 0 K for the motif in Fig. 1(b) (1500 atoms and force field I) with a 2 × 2 × 2  Fig. 1(b) (1500 atoms and force field I). At 5 K, we compare metadynamics (MTD, in legend) 57 with the blue-moon ensemble (BME) 59 method. At 0 K, the minimum energy path was mapped with static restrained optimizations (OPT) and with the nudge elastic band (NEB) method 60 . replicated version of that system (a total of 12000 atoms) using different cutoff radius for the pairwise interactions. Good agreement is found for the different profiles.

Appendix D: Maximum size of avoided level crossing
It is also possible to estimate the maximum size of the avoided level crossings, S max , expected in microwave spectra due to the qubit interaction with the TLS 9,26 where h is the Planck's constant, w is the width of the amorphous oxide barrier, C is the capacitance, and E 01 = δE is the energy difference between the two lowest qubit levels. Using typical parameters (w = 2 nm, C = 1 pF, Comparison of energy profiles at 0 K for the motif of Fig. 1(b) with a 2 × 2 × 2 replication of such system using different cutoff radius for the interactions. E 01 = 0.1 meV), Eq. (D1) gives a TLS-qubit coupling strength S max ≈ 60 MHz, in good agreement with the experimental data (∼ 100 MHz) 8,9,[13][14][15][16]18,21 .

Appendix E: Calculation of the effective mass and tunneling splitting
To estimate the effective mass associated to the motif, we used three methods: (1) We used the equipartition theorem on the collective variable δ Here, the angular brackets denote a thermal average over the canonical ensemble. The variables T δ andδ are the temperature and the velocity of the collective variable, respectively. The latter was evaluated analytically. We are aware that application of equipartition theorem to non-ergodic systems such as amorphous alumina is questionable.
(2) We performed a  (Fig. 1(b)). The effective mass of the collective variable δ = | r Al−O 2 | − | r Al−O 1 | was estimated using several techniques: The equipartition method given by Eq. (E1), the spectral method, and Eq. (E2). The two missing points at about δ = −0.3 Å were artifacts due to a "kink" in the minimum energy path and were eliminated from this graph.
Fourier transform of the time evolution of the collective variable δ (t) from a microcanonical MD trajectory to get characteristic frequencies, ω, in the lower part in the spectrum. For the motif in Fig. 1, we found a large vibrational peak at ν = ω/2π ∼ 4.12 THz. Then, from the known curvature of one of the wells (κ = 3.62 J m −2 ), we estimated the effective mass according to the harmonic approximation: m eff = κ ω 2 . (3) The third and last method is based on classical mechanics. The effective mass associated to a collective variable δ is given by where the sum is taken over all atoms N, each with Cartesian coordinates x i , y i and z i . The partial derivatives in Eq. (E2) were computed numerically by finite differences. With these three methods, we estimated the effective mass for the motif shown in Fig. 1(b). This motif was found using the force field II on the 360-atoms system. The results are in Fig. 8. Using only the coordinates of main three atoms (labelled as O1, Al, O2 in Fig. 1(b)) in Eq. (E2) yields m eff ≈ 4 a.m.u., which is consistent with the other two local methods (equipartition and spectrum, in legend). However, this value seems to underestimate the truly effective mass (m eff > 9 a.m.u.) because is missing the contribution from the rest, which is obtained using a summation over all atoms in the system (360, in legend). We note that for the local methods, m eff is weakly dependent on the collective variable δ, whereas if we consider all 360 atoms in Eq. (E2), the dependency becomes more pronounced. On the left well (negative δ), the effective mass for the full system was m eff =9 a.m.u., which is consistent with the expected reduced mass for a Al-O pair (10.04 a.m.u.), and increases as δ increases. This is due to the fact that as we move along the reaction coordinate δ, the differences between adjacent snapshots progressively increases.
Based on the obtained potential energy profiles at zero temperature, the characteristic barrier height V 0 of the TLS can be calculated by instanton theory 64 . Due to the heavy atoms involved, we can safely assume that for qualitative estimates, the instanton path is very similar to the calculated classical minimum energy path. Therefore, the collective variable corresponding to motion along this instanton path can be chosen the same as the generalized coordinate q describing the activated transition between the two potential energy wells. Then, the optimized action integral reads where the integration limits correspond to the turning points of the classical trajectory and m eff (q) is the effective mass. The tunnel splitting is given by δE = F exp(−S 0 / ), where the "fluctuation" factor F is of the order of the vibrational energy Ω in the potential energy wells. In general, neither the potential V(q) nor the effective mass m eff (q) in Eq. (E3) are constant along the trajectory. However, we can approximate S 0 by using V(q) ∼ V 0 and taking the effective mass equal to the mass of an aluminum atom, m eff (q) ≈ 27 a.m.u.. Using the interwell distance d = 0.5 Å and a typical vibration frequency Ω = 10 13 s −1 , we estimate that for an experimental splitting δE = 0.1 meV, the characteristic barrier should be about V 0 ≈ 8 meV for alumina. Thus, a temperature below 100 K should be enough to discern between thermally activated classical transitions between minima of TLSs and tunneling.
Once the effective mass of the collective variable was estimated, we fitted a symmetrized version of the energy profile at 0 K in Fig. 1(b)to a double-well like polynomial V (δ) = D δ 2 − δ 2 0 2 . The fitted values were D = 3.28 kcal/(mol Å 4 ) and δ 0 = 0.45 Å. We then solved the one-dimensional Schrödinger equation numerically using the code OCTOPUS 65 to get the tunneling splitting between the ground and first excited states for a tunneling particle of mass m eff ∼ 15 a.m.u.. In the OCTOPUS calculation, we used an ultrafine grid spacing (0.01 a.u.) on a one-dimensional box of length 16 a.u..

Appendix F: Other corundum materials: amorphous hematite
Here, we address the question of whether the motifs found in amorphous alumina could also exists in the glassy state of other corundum-like materials. To this end, we performed MD calculations of amorphous hematite (Fe 2 O 3 ) at low temperatures. The amorphous configuration was generated from the common crystalline phase α-hematite at its experimental lattice parameters (a = b = 5.038 Å, c = 13.772 Å, α = β = 90 • , γ = 120 • and density=5.26 g cm −3 ) following the same melting-quenching protocol used for alumina. An amorphous configuration of 1500 atoms (300 Fe 2 O 3 units) was equilibrated in a orthorhombic cell with final volume 25.2762 × 25.2762 × 27.5973 Å 3 , giving a mass density of 4.517 g cm −3 at zero external pressure. For the MD calculations on Fe 2 O 3 , we used the reactive force field REAX 66 as implemented in the LAMMPS code 53 together with its colvars module 67 for the free energy calculations. We used a time step of 0.2 fs and performed the charge equilibration [68][69][70] at every MD step up to a tolerance of 10 −6 eV. After several nanoseconds of MD simulation, we found only the motif shown in Fig. 9 (Top) that resembles the transformation between a trigonal-pyramidal to a tetrahedral geometry. This transformation, which is a change of coordination number (from 3 to 4) of the central iron atom, can be followed by a single reaction coordinate consisting of the distance between the involved iron and oxygen atoms. Fig. 9 (Bottom) shows the computed free energy profiles of this process at several temperatures. This figure features a more complex multilevel landscape than the previous rearrangements found for amorphous alumina. Interestingly, lowering the temperature bias the profile towards a single well corresponding to the trigonalpyramidal configuration. Thus, it is not clear whether this structural motif behaves as a TLS and remains to see whether amorphous hematite actually has TLSs.