{'fraction_non_alphanumeric': 0.026296944621260344, 'fraction_numerical': 0.014759707192870782, 'mean_word_length': 4.84173832210086, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A 3D investigation of delocalised oxygen two-level defects in Josephson junctions

Environmental two-level systems (TLS) have been identified as significant decoherence sources in Josephson junction (JJ) based circuits. For such quantum devices to be functional, the removal or control of the TLS is a necessity. Understanding the microscopic origins of the 'strongly coupled' TLS type is one current path of investigation to that end. The delocalized oxygen model suggests the atomic position of an oxygen atom is spatially delocalized in the oxide forming the JJ barrier. In this report we extend this model from its previous 2+1D construction to a complete 3D description using a Wick-rotated time-dependent Schrodinger equation to solve for time-independent solutions in three dimensions. We compute experimentally observable parameters for phase qubits and compare the results to the 2+1D framework. We devise a Voronoi classification scheme to investigate oxygen atoms delocalizing within strained and non-strained crystalline lattices, as well as realistic atomic positions a JJ amorphous tunnel barrier constructed in previous density functional studies.


Introduction
Decoherence sources such as the environmental two-level system (TLS) 1, 2 are presently a major hurdle to the realisation of operational superconducting qubits and other Josephson junction based quantum devices. Many properties of the TLS have been probed experimentally, demonstrating relatively long decoherence times whilst showing them to be stable and controllable. [3][4][5][6][7][8][9] However, identifying their true microscopic nature remains an open problem. Controllable superconducting qubit architectures (charge, flux and phase) make it possible to study single 'strongly coupled' TLSs, 4,6,7 as opposed to weakly coupled ensembles of TLSs that may be responsible for 1/ f noise. 1,10 Although the Standard Tunneling Model (STM) 11,12 can be used to explain all existing experimental results, to improve devices we need a detailed microscopic understanding of these defects.
Many microscopic models have now been proposed: delocalised oxygens, 13,14 surface state interactions, 15,16 polaron dressed electrons 17 and alien species (primarily hydrogen); [18][19][20] all of which are possible decoherence sources and/or explain the origin of strongly coupled TLSs. Better fabrication methods, higher vacuums and more robust qubit designs have historically suppressed or diminished the response of such noise sources. [21][22][23][24][25] However, continuing down the path of a perfectly clean but amorphous dielectric may no longer be the optimal direction. Whilst crystalline layers are difficult to construct in this environment, they are possible and show a dramatic decrease (of up to 80%) in visible TLSs over their amorphous counterparts. 26 Such a decrease could be explained via our previous work on the delocalised oxygen model, 13,14 where we suggest the origin of some TLS defects are oxygen atoms in a spatially delocalised state due to the non-crystalline nature of the dielectric layer. As such, this is a concrete microscopic ensemble of the classic atomic tunneling model: typically cited as moti-vation for the STM (which is more general, as the tunneling object need not be a single atom).
Our previous model discussed a simplified configuration space using low dimensional arguments, culminating in a 2+1D framework where a cage of six aluminium atoms surrounded a central oxygen atom in three dimensions and calculated observations based on the oxygen delocalising on a two dimensional plane. The dimensional simplicity of this model allowed a direct matrix diagonalisation approach to solve a time independent Schrödinger equation, which due to memory limits could not be extended into a complete 3D representation of a TLS defect.
Here we introduce an alternative method for this model which can obtain a complete three dimensional picture of the delocalised oxygen ground and excited states. We use a Wick-rotated time-dependent Schrödinger equation to solve for time-independent solutions in three dimensions (see Methods for more detail). As the work in this report is a higher dimensional extension of the 2+1D model, its terminology has been retained herein. It is suggested the reader be familiar with the discussions in Ref 14, although a short summary of the nomenclature is outlined below. Using this more powerful technique, we are able to move beyond a simplistic cubic representation and investigate true lattice configurations and realistic amorphous positions from our DFT studies. 13,27 Two defect 'types' are investigated, both of which stem from deforming the crystalline lattice of corundum: the low temperature and pressure phase of Al 2 O 3 (the closest Al-O bond distance in this phase is ∼1.85Å). Defect type A; where the aluminium-oxygen bond distance is shortened, forcing the oxygen to occupy one of two off-axis positions, and defect type B; where the opposite occurs: the aluminium-oxygen bond distance is lengthened, allowing two preferred oxygen positions, on-axis.
The trigonal lattice structure of this crystalline material and the amorphous nature of a JJ barrier give rise to many complications for a microscopic model, due to potential contributions from many neighbouring atoms to a possible defect site. Therefore, to fully understand the response to these interactions, we start with a very simplistic cluster of atoms and add complexity once an understanding of rudimentary behaviour is achieved. Defining the oxygen position to be at an origin, two introduced aluminium atoms can be considered as pairs (x = −X, +X = {±1.85Å}) lying on a cardinal axes; which are identified as |X| or simply referred to as the 'defect pair'. Displacing |X| equidistantly from this origin (i.e. moving away from optimal crystalline configuration) will yield either an A or B type defect, depending on the direction of displacement. |Y | and |Z| describe aluminum pairs in the y and z directions respectively. Eigenenergies of a system are frequently discussed in a relative fashion using a convention where E i j = E j − E i , such that the ground state (E 0 ) to first excited state (E 1 ) energy splitting is defined as E 01 .

Results
The capability to investigate oxygen delocalisation in 3D, using realistic atomic positions within an amorphous layer of Al 2 O 3 generates an extensive state space. Hence we start the investigation using a simple three atom system comprising of an oxygen atom and a confining aluminium pair: |X|. An eigenspectrum of the ground and five lowest excited energy levels of this system is generated by increasing or decreasing the pair separation distance from the corrundum lattice distance 1.85Å which is depicted in Figure 1 (for comparison with the 2D model, see Figure 2 in Ref 14).
In this figure an (an)harmonic division about |X| ∼ 1.85 A separates two regions with dissimilar properties. In the region where |X| > 1.95Å, a degenerate E 01 is observed and as the separation distance has increased, can be identified as a possible B type defect. The second region (|X| < 1.55Å) has a tightly bound eigenspectrum compared to the rest of the map although this yields no degeneracy.
In three dimensions, the potential minima manifests as a sphere around the location of an atom. As a consequence the three atom Al-O-Al chain produces a unique, rotationally symmetric ground state which is reminiscent of an oxygen interstitial defect in crystalline germanium 28 and depicted in Figures 2a&b. Comparatively, this region in the 2D case indicated the presence of an A type defect, as the 2D potential minima ring would be projected onto a plane, collapsing a degree of freedom.
Two other observations concerning the 2D to 3D transition is the dramatically reduced total energy values as wavefunctions are no longer artificially confined, and the higher excited states in the B type region now exhibit a quad degeneracy rather than two split doubly degenerate pairs -again because states are not artificially confined in the extra dimension.
Eigenspectrum of a three atom system Al-O-Al, over a varying distance separation |X|. Each excited state has been normalised with the ground state, which shows a clear E 01 degeneracy at large separation distance -indicative of a B type defect. An intermediate region exists as the separation distance decreases, which is approximately centered about the optimal Al-O bond distance of corundum (1.85Å). At small separation distances, the eigenspectrum is distributed over comparatively small energy levels, but is not degenerate (see Figure 2). This simplistic three atom description can now be built upon to understand the interactions of a delocalised oxygen in a three dimensional space. A confining pair is introduced in an additional dimension to the configuration in Figures 2a&b, and the unique ground state is lifted to a degenerate E 01 pair; observing an A type TLS defect. Figures 2c&d illustrate this effect with a confining pair |Y | = 3.5Å. This distance is an arbitrary choice, as the shift occurs even at the cut-off limits imposed by the potential portion of the model 29 (see Methods). In other words, any additional potential contribution which does not share the axial symmetry of the Al-O-Al chain will lift the degeneracy and result in a TLS.

Oxygenic Orbitals
A complete picture of the ground and excited state wavefunctions in 3D is now possible, which allows us to further investigate the properties of the two-level system and illustrate the importance of crystalline dielectrics in future Josephson junction devices.
A confined harmonic state can exist with many atomic configurations, and as with the low dimensional model, setting variables symmetrically is one of these cases. Figure   approximation holds and the system is not considered as a TLS candidate. From this localised case, we extend the |Y | confinement out to 2.1Å. Using the Figure 1 eigenspectrum we can predict this configuration to be a B type defect; where the |Y | separation distance is increased, and an oxygen dipole forms parallel to the y axis. This is indeed the case as shown in Figure 4.
A particularly complex phenomenon emerged from the low dimensional model: quad degenerate ground states, where the energy difference between two degenerate pairs E 01 and E 23 approach the difference of the pairs themselves, i.e. E 01 = E 23 ≈ E 12 . Figure 5 shows a configuration expressing this behaviour.
The discussion in Ref 14 indicated that this form of degeneracy has a low probability of occurrence compared to A or B type defects. Additionally, if E 12 ≥ 100 GHz the higher states (E 2 and E 3 ) can be ignored completely as energies of this magnitude dissociate Cooper-pairs, hence can be viewed as an operational upper bound for Josephson junction devices. Splitting energies in 3D are much smaller than the 2+1D case, which suggests that the probability of these defects being experimentally visible may be higher than predicted with the low dimensional model based solely on this measure. However, a comprehensive analysis of the three dimensional configuration space would need to be undertaken before real estimates of this behaviour can be given.
Another interpretation of this phenomena has to do with the energy resolution of the coherent single-pulse resonant driving spectroscopy 7 used to obtain qubit-TLS coupling (S max ) values. The minimum resolvable energy using this process is approximately 1 MHz, which may indicate that S max anticrossing measurements are manifestations of qubit-TLS couplings at the E 12 level, if a configuration such as Figure 5 is scrutinised. E 01 and E 23 splitting levels would be hidden in their degenerate subspace and manifest as an effective twolevel defect with E 01 = 31.26 GHz.

A type defects and dipole considerations
Whilst adding additional confinement pairs causes the unique, rotationally symmetric ground state of |X| < 1.55Å (Figures 2a&b) to become degenerate; the third dimension yields a complication in the dipole measurement for the A type region.
Consider a system with parameters |X| = 1.53, |Y | = 1.97, |Z| = 1.95Å, illustrated in Figure 6. This system exhibits TLS behaviour, with E 01 = 40.63 MHz, and a dipole strength in y of 0.26 eÅ: perpendicular to the confining x axis as expected. The leftmost axis of Figure 6 shows a 2D projection of the first excited state, illustrating no major differences in the response of the 3D and 2+1D models.
However, a small change in the y confinement alters the system in a non-trivial manner. Moving |Y | from 1.97 to 1.9 A (for example) crosses a bifurcation in state space. As |Z| is now the least confining pair, the dominant dipole direction flips to z as shown in Figure 7.
Ultimately this does not add any complexity to the model:  minimal confinement in y generates an A type defect with a dipole in y, perpendicular to x; minimal confinement in z generates an A type defect with a dipole in z, perpendicular to x. As strongly coupled defects are experimentally identified via avoided level crossings in qubit spectroscopic diagrams, 7 the dipole moment of the defect couples to the electric field across the junction. 22 Therefore, to be identified as a TLS, the defect must be aligned to this field. It follows then, that the dipole directions for an A type defect in this model can be considered equivalent. Additionally, as with B type and any other possible TLS defect, both A type alignments will not appear on a spectroscopic scan as an avoided level crossing if their dipole is perpendicular to the external field.

Defects in corundum
Results such as Figures 1 and 3 identify configurations which exhibit harmonic eigenspectrums and localised oxygenic positions. This behaviour is expected when an oxygen atom is confined appropriately and has no reason to tunnel. One way to check this assumption (and consequently further validate the TLS model) is to calculate the potential landscape seen by a single oxygen atom from bulk corundum (α-Al 2 O 3 ). The TLS model, delocalising an oxygen across all potential space about its minimum energy position, should yield a localised, harmonic wavefunction positioned at the location where the oxygen was removed (i.e. the lattice position associated with the vacancy).
An oxygen atom in a corundum lattice is identified by a Voronoi polytope (see Methods) and a nearest neighbour atom cluster is used to calculate the wavefunctions and energies of an oxygen atom located inside the potential landscape generated by the corundum lattice. As predicted, the result is a localised oxygen atom locked into the lattice position of the vacancy, and the oxygenic orbitals behave in a similar fashion to the harmonic system calculated previously in Figure 3 (data |1 |2 x y x y z E 01 = 40.63 MHz E 0 E 12 = 1.14 GHz not shown).

Straining the Local Lattice
In an attempt to put bounds on the amount of amorphousness required for a near-crystalline structure to exhibit TLS properties, we apply a localised strain tensor on the defect cluster and observe its effects. Two strain ratios ∆L L = ∓0.019Å representing a compressive force along the x axis, and a tensile force also along x respectively, attempt to simulate a local strain from the localised Al-O separation distance of 1.85Å in the direction of the A type and B type regions respectively (see Figure 1). The original crystalline system has a ground to first excited state splitting of E 01 = 13.48 THz. Comparing the E 01 values of the strained results (compressive: 35.42 THz, tensile: 34.84 THz) to the crystal result, we see an increase in splitting energy (data not shown). This response is opposite to the simple straining of Figure 1 where the E 01 energy was a maximum at the crystalline Al-O separation distance. This suggests that local Poissonian strain does not give rise to immediate TLS behaviour-at least along the sampled axis x.
One possible explanation of this result comes from our previous study concerning strain on the 2+1D model: 14 there exists a preferred strain axis which responds rapidly in comparison to other directions, although this axis is difficult to acquire directly from a completely harmonic case. It is therefore pertinent to investigate the amorphous configurations in more detail.

Defects in Josephson Junctions
Applying the Voronoi classification scheme the myriad atomic configurations within the Josephson junction models constructed in Ref. 27 allows us to investigate realistic amorphous clusters residing within a JJ tunnel barrier.

|1
|0 |1 |2 x z x y z E 01 = 4.14 GHz E 0 E 12 = 1.49 THz Using the structure with an AlO 1.25 amorphous barrier with a density 0.75 that of corundum we classify clusters for all oxygen atoms extant in the tunnel barrier. It is posited that a large void space around an oxygen may allow it to delocalise over a large region and manifest as a B type defect (however, this is pure speculation). The volume of each Voronoi cell can be calculated (as the cell is essentially a convex hull of some vertices describing a Delaunay tessellation) and we choose the five cells with the largest volume to investigate. Another conjecture is that an asymmetric cell (obtained by some form of sheer in the local cluster configuration) yields TLS behaviour, thus a final cluster is investigated which has the most asymmetry -identified using the length parameter of the spherocylinder formalism outlined in 30. The five cells with the most asymmetry were also identified; these clusters and their properties are displayed in Figure 8 along with the large volume candidates.
As with the cluster results of the crystalline and Poisson strain investigations, the few amorphous clusters studied failed to identify any TLS behaviour. Whilst each oxygen displays a different spectra, none of the E 01 splitting are degenerate and dipole strengths are no larger than 0.034 eÅ (from the red cell, identified as the oxygen with the most asymmetry in its cluster) in any cardinal direction.

Discussion
Using the direct diagonalisation method in three dimensions is seemingly intractable with current memory limitations, although this does not exclude investigations in the 3D domain. The Wick-rotated time-dependent Schrödinger equation method has proved to be a valid tool to study this space, albeit with computational time disadvantages. Whilst the 2+1D model allows us to investigate phase space in a much more flexible way, one of the major variations in the models is the final energy values, as the 3D model is no longer collapsed into fewer dimensions. This result suggests that if one is attempting to identify exact values for a delocalised oxygen TLS, they should be utilizing the 3D method. Contrarily, if trends are required, the 2+1D model provides the more efficient solution. Whilst splitting energies and dipole results may not completely reflect the atomic positions of the surrounding clusters in this case, it is evident from discussions like those concerning Figures 6 and 7 that there are 2D to 3D equivalencies. Using a Voronoi classification scheme to identify atom clusters which describes the local potential environment of an oxygen atom is shown to be a powerful method, describing a lattice position in crystalline corundum as well as the properties of any selected oxygen atom in a model junction.
The free parameter concern stated in Ref. 14 is somewhat mitigated through this process, as clusters with many atoms can be described with one Voronoi cell (clusters of up to 24 atoms are described in Figure 8 for example).
The affect of local strain on a crystalline cluster observed behaviour which does not align with the simple strain investigations undertaken with the idealised 6 atom model, suggesting that the deformation axis chosen is not the preferred strain axis identified in Ref. 14 and that localised strain may not lead to TLS behaviour directly.
With such a small sample size it is difficult to draw any direct conclusions from the results in Figure 8. Other classification metrics or methods of identification such as extending the Euclidian distance Voronoi classification discussed here to the Voronoi S-networks considered in 30 may yield TLS behaviour with the present model.
However, the more probable explanation is that the surrounding atom clusters exert enough confinement energy to the oxygen atoms to keep them spatially localised, unlike the simpler 6 atom confinements. One phenomena which may assist in successful TLS classification, is the inclusion of many body interactions (the C type defect in Figure 1 of Ref. 14). As oxygen atoms tunnel through potential barriers, the surrounding atomic neighbours may deform and create a multiatom TLS candidate. Such behaviour has been studied previously, [31][32][33][34][35] and multi-body identification approaches have recently been investigated. 36 Generalising the delocalised oxygen atom to a multi-body description is a topic for future work.

Methods
To identify A and B defect types in either amorphous or strained crystalline systems, an effective single particle Hamiltonian can be used. This simplifies the system by ignoring any time evolution properties or many-body interactions and yields the equation where m oxy is the mass of an oxygen atom and V (r) is the potential due to the surrounding lattice. This potential is modelled by the empirical Streitz-Mintmire potential, 29,37 which adequately captures the variable oxygen charge states occurring when oxygen is present in a predominantly metallic environment like a JJ metal-oxide interface. A complete description of the parameters used can be found in Ref 14 as the implementation remains consistent with this model.
To obtain full 3D eigenstates of this Hamiltonian we implement a Wick-rotated time-dependent Schrödinger equation to imaginary time t = iτ, 38 which yields a general solution to the wavefunctions Ψ(r, τ) = ∞ ∑ n=0 a n ψ n (r)e −E n τ .
Here, a n are coefficients based on initial conditions of the system where n = 0 indicates the ground state, n = 1 the first excited state etc. and E n is the corresponding eigenenergy. As E 0 < E 1 < E 2 < . . ., evolving Eq. 4 to large imaginary time will provide a good approximation to the ground state influenced by the time-independent potential V (r).
This solution is obtained by numerically approximating the spatial derivatives with finite differences. For stability and precision, we descritise the space using a 7-point central difference method where f k = f (r 0 + kh), calculated with a step size h = 0.01 A.
Excited states are obtained by Gram-Schmidt orthogonalisation: 39,40 choosing eigenfunctions within a degenerate subspace to be orthogonal to the system's basis. For example, consider a converged, orthonormal ground state |ψ 0 and an initial, non-orthonormal guess for the first excited state |ψ 1 . Subtracting the projection of the excited state along the ground state from the initial guess yields a vector orthogonal to the ground state This resultant vector can then be normalised and evolved through the same process as the ground state until convergence is achieved. Similarly, the second excited state requires the same treatment, although projections along both |ψ 0 and |ψ 1 are needed. An active repository of the software which implements these methods can be found via Ref 41. When dealing with the corundum crystal and amorphous JJ barrier clusters, we apply a Voronoi tessellation 42 to the lattice coordinates of each system using a Euclidean distance metric. A polytope representing the convex hull encompassing an origin position symbolises the oxygen vacancy. Atoms situated within polytopes sharing an edge with the origin polytope are selected as the model cluster. These atoms constitute a first order screening of the surrounding lattice potential, representing an acceptable approximation to higher order screenings (e.g. also including second nearest neighbours, which add appreciable computational intensity). Clusters identified by this method use the Streitz Mintmire potentials generated by GULP, 43 which also account for monopole-monopole and self energy interactions to counteract any net charge generated by selecting an arbitrary cluster of atoms from a lattice with periodic boundaries.
Arbitrarily straining a corundum crystal, which has trigonal symmetry and six elastic constants, 44 will exert force on the structure in a non-trivial manner and act on many different crystal planes. The materials science complications of this phenomena are well beyond the scope of this investigation, and instead we approximate a local strain by the Poisson effect 45 using the values of the adiabatic bulk and shear moduli to obtain a Poisson's ratio ν = 0.234. This value is calculated as an equivalent isotropic polycrystalline aggregate to Sapphire (a mineral variety of corundum). 46 The relationship can be used to strain our atom cluster for example, with a length increase of ∆L in the x direction, and a length decrease of ∆L in the y and z directions.