Atomic delocalisation as a microscopic origin of two-level defects in Josephson junctions

Identifying the microscopic origins of decoherence sources prevalent in Josephson junction based circuits is central to their use as functional quantum devices. Focussing on so called"strongly coupled"two-level defects, we construct a theoretical model using the atomic position of the oxygen which is spatially delocalised in the oxide forming the Josephson junction barrier. Using this model, we investigate which atomic configurations give rise to two-level behaviour of the type seen in experiments. We compute experimentally observable parameters for phase qubits and examine defect response under the effects of applied electric field and strain.

Identifying the microscopic origins of decoherence sources prevalent in Josephson junction based circuits is central to their use as functional quantum devices. Focussing on so called "strongly coupled" two-level defects, we construct a theoretical model using the atomic position of the oxygen which is spatially delocalised in the oxide forming the Josephson junction barrier. Using this model, we investigate which atomic configurations give rise to two-level behaviour of the type seen in experiments. We compute experimentally observable parameters for phase qubits and examine defect response under the effects of applied electric field and strain.

I. INTRODUCTION AND CONCEPT
Superconducting qubits and Josephson junction based quantum devices in general are often limited by decoherence sources. A common and important source of decoherence is the environmental two-level system (TLS) 1,2 . Experiments have probed these defects directly and shown them to be stable, controllable and have relatively long decoherence times, [3][4][5][6][7] although little is known about their true microscopic nature. Many phenomenological theories attempting to describe them exist; including charge dipoles 8 , Andreev bound states 9 , magnetic dipoles 10 , Kondo impurities 11 and TLS state dependence of the Josephson junction (JJ) transparency. 12 Detailed fitting of experimental data can place limits on these models 13 but the scope of free parameters of each model allows them all to fit experimental results -rendering them presently indistinguishable. It is therefore important to construct microscopic models of these systems to increase our understanding of their composition. Polaron dressed electrons 14 and surface aluminium ions paramagnetically coupling to ambient molecules 15 are two recent models that may shed new light in this area.
Bistable defects in glasses and amorphous solids in general have been understood for some time 16 . Amorphous insulating barriers (either in the form of Josephson junctions or simply a native oxide) form an integral part of superconducting circuits, so it comes as no surprise that TLSs are often considered to be an important source of noise in these circuits 1,2,8 . Developments in controllable qubit architecture (charge, flux and phase) has enabled the study of so-called 'strongly coupled defects' 4,6,7 . These defects have comparable resonance frequencies to the qubit circuit and coupling strengths as well as decoherence times long enough to allow coherent oscillations between the qubit and TLS. Probing individual defects has promoted their bistable nature from hypothesis to observable fact, as well as providing clues to their microscopic origin.
As described in previous work 17 , we consider the origin of some defects to be within the amorphous oxide layer itself, specifically an oxygen atom in a spatially delocalised state. This has important ramifications for materials sci-ence based efforts to reduce the effects of TLSs. If, as alternative models suggest, TLS defects indeed stem from surface states 18 or the accidental inclusion of an alien species 19,20 ; future fabrication techniques or more robust qubit designs may suppress or diminish the response of such noise sources as has been the case historically 8,[21][22][23][24] . Nevertheless, if the oxygen atoms themselves form a noise source, a perfectly clean but amorphous dielectric may still harbour a large ensemble of TLSs.
Our approach considers only atomic positions as input parameters in an attempt to construct a microscopic model rather than a phenomenological one. The premise of our TLS model is that positional anharmonicity of oxygen atoms arises within the AlO x barrier of the Josephson junction primarily due to its non-crystalline nature. As an illustrative example, consider an interstitial oxygen defect in crystalline silicon: the harmonic approximation for atomic positions cannot be applied due to the rotational symmetry of the defect as oxygen delocalises around the Si-Si bond axis 25 . This forms an anharmonic system with a quasi-degenerate 17 ground state, even in a "perfect" crystal. This ansatz allows the existence of many different spatial configurations throughout the layer, causing unique TLS properties based solely on atomic positions and rotation in relation to the external electric field. To simplify the configuration space we initially consider the introduction of small lattice irregularities from an idealised trigonal-like 2D oxide lattice. Possible defects of this nature are depicted in Figure 1.
It has often been suggested that a bath of TLSs are responsible for the weakly coupled, ohmic 1/f noise 1 . However, it is unclear if the identified strongly coupled systems are from the same origin. Ultimately many separate microscopic suspects may be identified; although work in this area is not mature enough to speculate. The model in this work therefore only attempts to describe TLSs that are strongly coupled in nature.
The outline of this paper is as follows. Section II introduces the theory from which we build our model to investigate delocalised oxygen based two-level systems. The following three sections then start from a minimal example of this model and slowly add complexity, so interactions can be examined and understood in a system-arXiv:1408.5687v1 [cond-mat.mes-hall] 25 Aug 2014  26,27 . A: the bond length is shortened causing an oxygen to form a dipole perpendicular to the bond axis, B: the bond length is lengthened causing an oxygen to form a dipole parallel to the bond axis, C: a cluster of three oxygens are rotated about a central metal atom. Types A and B are of interest to this work. Type C requires a many body investigation that is beyond the current scope, although it has been discussed extensively in the literature 28-31 . atic way. Section III considers a defect comprising of an Al-O-Al chain in one dimension, perturbed from a crystalline lattice, simulating defects A and B in Figure 1. In reality, an oxygen atom in the amorphous layer of a Josephson junction will be surrounded by atoms in all three dimensions. Moving towards a model representation of this, Section IV extends the model to two dimensions, with four aluminium atoms confining an oxygen in a plane. The completed model is then described in Section V, which extends oxygen confinement into three dimensions (with six aluminum atoms), whilst still projecting oxygen delocalisation on a plane. Although in general an oxygen can delocalise in three dimensions, for this investigation we focus on an effective 2+1D model, minimising both computational and descriptive complexity while still modelling the relevant behaviour of the system. The following sections then apply the 2+1D model and compare delocalised oxygen responses to experimental TLS data. Section VII discuses qubit-TLS coupling and Section VIII observes the effect of mechanical strain.

II. METHODS
We are interested in what happens to the oxygen atom as it responds to an external potential exerted via its' nearest neighbours, particularly if the bonded aluminium atoms are displaced in a manner similar to the defect types outlined in Figure 1. If we also ignore any time evolution properties of the system for now and derive an effective single particle Hamiltonian where m oxy is the mass of an oxygen atom and V (r) is the potential due to the surrounding (mostly amorphous) lattice.
The two most striking properties of a strongly coupled TLS, assuming the charge dipole framework is physically representational, are its ground to first excited state splitting E 01 and a strong electric dipole moment. Observed values of E 01 differ between qubit architectures: approximately 1-10 GHz for transmons 22 , 4-5 GHz for flux qubits 6 (although recent designs can tune this gap down to the MHz range 32 ) and nominally 7-8 GHz for phase qubits 13 . Dipole moment strengths also vary, but are usually on the order of 1 eÅ 5,13 . Whilst many other properties have been measured, this work will focus on obtaining respectable values for E 01 and dipole strength, assuming our model defect is embedded inside a fictitious phase qubit. When a representative example is required in this work, E 01 = 8 GHz and a dipole strength of 1.2 eÅ will be used, to compare directly with a TLS studied in Ref. 13.
A potential which represents the junction, requires a number of capabilities. It needs to describe interactions between atomic species (in this case, Al-Al, O-O and Al-O interactions), as well as many body interactions to be accurate as possible (which will be required to investigate quasi-degenerate states). As a complete description of many body interactions does not currently exist; potentials of this type tend to be empirically fitted to experimental data in order to obtain physico-chemical properties of a studied system as accurately as possible. The trade off here is that whilst any given potential may describe some properties of the system well because of certain fitting parameters, other properties may be well out of range as they were not included in the study that constructed the potential. Great care was taken to choose the best potential that accurately represented the junction, in this case the empirical Streitz-Mintmire potential 33 , which describes a myriad of aluminium oxide properties over quite a large range of temperatures and pressures with high accuracy when compared against similar potentials. It was also chosen over simpler fixedcharge models 34,35 due to the complex geometry of the Josephson junction. A variable charge potential such as Streitz-Mintmire can capture the variable oxygen states when present in a predominantly metallic environment through the minimisation of the electrostatic term in the potential (Eq. 2). This capability is particularly important here, as our junction has two metal-oxide interfaces, and our TLS defects may reside close to these boundaries. The Streitz-Mintmire potential is given by where q i,j is the atomic charge, . J 0 i is an empirical parameter known as "atomic hardness" or a self-Coulomb repulsion 36 , δ ij = 1 when i = j and δ ij = 0 when i = j, and all summation is calculated for N atoms of the target system. The square bracket notation represents Coulomb interaction integrals between valence charge densities and/or effective core charge densities and take the form 37 : with a = i, j; b = i, j; a = b in Eq. (2). dV a,b are integrating volume units. r av is therefore the center distance between atom a and dV b , and r vv is the center distance between dV a and dV b .
The first term in Eq. 2, E EAM , does not depend on the partial charges q i and therefore describes a chargeneutral system, represented here with a quantum mechanical based empirical embedded atom model (EAM) for the Al-Al and Al-O interactions with F i [ρ i ] as the energy required to embed atom i in a local electron density ρ i , and φ ij (r ij ) describing the residual pair-pair interactions by way of Buckingham and Rydberg potentials where r ij is the interatomic (Euclidean) distance between atoms i and j, all other constants are listed in Table I. Further formalisms and parameters can be found in Refs. 33 and 37, implementation is also discussed in Ref. 38. Using this potential, a single body time-independent Hamiltonian is constructed using a 7-point central difference method where f k = f (x 0 + kh). The step size h = 0.1Å was found to be optimal for the grid range relevant to our purposes 39 . The resulting matrix is then diagonalized, obtaining eigenvectors and eigenenergies with precision better than 10 kHz. In contrast, the energy scale for JJ defects observed in experiments is 10 GHz 4,6,7 .
The dipole element is computed using numerical integration of the ground-and first-excited states (ψ 0 , ψ 1 ), where is an example of the dipole in the x direction. Relative differences in energy levels (i.e. energy splittings) are an important measure of the model. We therefore define a convention where E ij = E j − E i , such that the ground state (E 0 ) to first excited state (E 1 ) energy splitting is defined as E 01 . For comparison with experimental results, energy is expressed in frequency units throughout this paper.

III. TLS DEFECTS AS PERTURBED BOND ANGLES IN A LATTICE
Consider the two simplest cases in Figure 1 -defect type A; where the aluminium-oxygen bond distance is shortened, forcing the oxygen to occupy two off axis positions, and defect type B; where the opposite occurs: aluminium-oxygen bond distance is lengthened, allowing two preferred oxygen positions on axis.
These two defect types can be modelled by solving our system hamiltonian (Eq. 1) for three atoms: an oxygen with two aluminium atoms at a lattice coordinate apart, displaced toward and away from the oxygen. For example, corundum: the low pressure and temperature phase of aluminium oxide (also known as α-Al 2 O 3 ) has an Al-O bond distance of ∼ 1.85Å. If we define the oxygen position to be at an origin, the aluminium atoms can be considered as pairs (x = −X, +X = {±1.85Å}) lying on a cardinal axes; which are identified throughout this paper as |X|. 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.
An eigenspectrum of the six lowest energy levels of this system over a continuum of values in |X| is depicted in Figure 2. Each energy is measured relative

FIG. 2. (Color online) Eigenspectrum of a three atom system
Al-O-Al, over a varying distance separation. Each excited state has been normalised with the ground state, which clearly shows two regions with a degeneracy at the lowest level. Section A is indicative of the A type defect (Figure 3), section B of the B type ( Figure 4). An (an)harmonic crossover point is also extant, labelled as section C, which is approximately centered about the optimal Al-O bond distance of corundum (1.85Å).
to the ground state, which shows two particular regions

FIG. 3. (Color online)
A top down view of the potential seen by an oxygen atom with two aluminium atoms closer than an appropriate lattice distance (truncated at 0 eV, gray). The ground state eigenmode of the oxygen can be seen in the center (green/yellow). This configuration has a separation distance of |X| = 1.5Å and is representative of an A type defect (see Figure 1). 1D potential values and projected wavefunctions, which are anchored to the ground state energy (thin gray line) are plotted in the outsets to better indicate the depth of the potential well.
where E 01 (green, solid) is quasi-degenerate (labelled sections A and B: both associated with the respective defect type). There exists a third (an)harmonic region (section C), which reaches a harmonic state at a separation distance of |X| ∼ 1.85Å: the optimal corundum Al-O bond distance. At this distance the spatial harmonic approximation holds and the oxygen can be considered to be localised.
To investigate the potential landscape and the resultant oxygen wavefunction in sections A and B of Figure 2, we choose two separation distances: |X| = 1.5Å (Figure 3), which lies in the A type defect section, and |X| = 2.2Å (Figure 4), which exists in the B type defect section. Although the Al-O-Al chain is arranged in a line, we consider delocalisation of the oxygen in two dimensions so that both defect types can be identified on a continuum separation in |X| (as Figure 2 depicts) rather than using two separate coordinate systems.
Both figures show a top down view of the potential exerted on the oxygen by the two aluminium atoms (gray), whose positions lie at the centre of the circles. We truncate potential energy values above 0 eV for clarity, where energy is measured relative to the zero point set by the Streitz-Mintmire potential's electronegativity correction. The ground state wavefunction (green/yellow) indicates the spatial probability of the oxygen atom in these configurations. It can be clearly seen in Figure 4, where the pair is far apart, two local minima exist in the form of rings around the base of each aluminium position: due to the fact that aluminium oxide is only partially covalent and mostly ionic.
Equivalences between the A and B type defects illus- with two aluminium atoms further apart than an appropriate lattice distance (truncated at 0 eV, gray). The ground state eigenmode of the oxygen can be seen in the center (green/yellow). This configuration has a separation distance of |X| = 2.2Å and is representative of a B type defect (see Figure 1). 1D potential values and projected wavefunctions, which are anchored to the ground state energy (thin gray line) are plotted in the outsets to better indicate the depth of the potential well. trated in Figure 1 are clearly visible: Figure 3 showing a shortened bond length, causing an oxygen dipole perpendicular to the bond axis, and Figure 4 depicting a lengthened bond and an oxygen dipole parallel to the bond axis.
The functional form of the 1D potential (outset axes of Figures 3 and 4, black; either V (x, 0) or V (0, y)) in the direction of the defect is a double well. Also depicted in these outsets are projected wavefunctions, which have been scaled for visual purposes, but anchored at the ground state energy (thin gray line). This representation is equivalent to the standard two level physics depiction of the TLS, and as such, potential offsets and tunneling matrix elements of our model may be estimated in this limit.

IV. TLS DEFECT CONFINED IN TWO DIMENSIONS
The ideal case discussed in Section III ignores many real world complications, in particular any potential constraints from nearest neighbour atoms that undoubtedly surround the TLS. In an attempt to add complexity to the model gradually, we start with two additional aluminium atoms on the same plane as Figures 3 and 4, confining the defect in the |Y | direction (i.e. y = −Y, +Y ).

A. Classifying Eigenspectrum Dynamics
As the values of |X| or |Y | are altered, a complex interplay between the excited states of the model can result. Simple two-level degeneracy and harmonic states are no longer the only possibilities. To interpret what is occurring in a certain domain, we define a metric using the ground and four lowest excited state energies This metric ranges from 0 to 2 and can give a qualitative understanding of the eigenspectrum of the defect.
To begin we plot a phase space diagram akin to those introduced in Ref. 17, where ξ is plotted as a function of the distance to the confining aluminium atoms (|X|, |Y |). Each phase diagram is split into at least four domains, where the properties of these domains can be explained through the interplay of potential configuration and dipole alignment (discussed in Section VI). Focusing for now on the influence of potential shape, the 2D potential can be approximated as two 1D potentials: one projected in the x direction and the other along y. There are two relevant configurations, a set of two double wells (tetra-well) or a set of a double and harmonic well (hemi-tetra-well); which are both illustrated in Figure 5.  The ξ metric is capable of identifying the tetra-(ξ = 0) and hemi-tetra-(ξ = 1/2) domains, as well as harmonic (ξ = 5/4) and unique ground state (rotationally symmetric, Mexican hat-like) (ξ = 2) regimes and finally the location of bifurcations or transitions (ξ = 1). Figure 6 shows the corresponding layouts of each interplay. It is worth noting that ξ = 3/2 can also be considered harmonic for the lowest three levels. 6. Energy level representation of the lowest five eigenenergies of a candidate defect and their associated ξ value. ξ = 0: the tetra-well domain, ξ = 1/2: the hemi-tetra-well domain, ξ = 5/4: the harmonic region and ξ = 2: the unique ground state region. The bifurcation/transition region, ξ = 1 is not clearly defined thus an example is not shown here.

B. Visualising Phase Space and Identifying TLSs
Using this ξ metric and varying the values of both |X| and |Y | generates Figure 7, a phase map of the interplay of low energy states of the oxygen atom confined in two dimensions by aluminium atoms. The x and y axes show the |X| and |Y | pair separation distances respectively over a range of 1.85−4Å and the phase colour indicates which ξ region a particular configuration exists in.
Regions that are deemed to be unimportant when searching for TLS behaviour are those where ξ ≥ 1, as these correspond to eigenspectra which do not possess a (quasi-)doubly degenerate ground state.
As discussed in Section III, in the one dimensional confinement case, both A and B type defects exist in a hemitetra-well (ξ = 1/2). This is also the case for the two dimensional landscape -contour lines corresponding to E 01 = {0.5, 2, 4, 6, 8, 10} GHz (red-yellow) on Figure 7 show the configuration properties that result in TLS E 01 splittings in various qubit architectures described in Section II. For example: anywhere an orange, 8 GHz line exists on this phase map, the |X|, |Y | coordinates of the model generate a cluster configuration that yields the same E 01 observed in phase qubit experiments 13 .
Many of these lines lie completely within a ξ = 1/2 (green) region. Consider the contour set on the left of Figure 7 where |X| 1.5Å: Section III states this distance is indicative of an A type defect. The addition of the |Y | aluminium pair does little to perturb the potential at larger distances (|Y | 3Å) and causes only slight deformations at smaller distances (2 |Y | 3Å). As the phase space is symmetric about |X| = |Y |, A type defects exist at the bottom of the plot where |Y | 1.5Å as well.
B type defects also exist in a ξ = 1/2 region, when |X| or |Y | 2.4Å. The orthogonal pair separation distance is at least 1Å larger than the defect pair in these configurations and therefore have no bearing on the potential seen by the oxygen. Complications arise when the orthogonal pair is closer and begins to confine the defect. To explain this response, a better understanding of the domains of the map is required.

C. Analysis of Phase Space Domains
The case of ξ = 0 is particularly interesting in this scenario: a tetra-well region, causing a quad degeneracy in the ground state. In this region, the characteristics of a B type defect are strongly modified. To understand why, we first explain the large rounded ξ = 0 domains in the centre of Figure 7's phase space.
Tetra-well domains exist when the confinement potential on the oxygen consists of two double wells (see left plot in Figure 5. This phenomena emerges when both the defect pair and the confining pair of atoms are close enough to interact. Consider an A type defect with a defect pair in |X| and a confining pair at a constant value Overlayed contour lines corresponding to E01 = 0.5−10 GHz (red to yellow) are comparable to existing experimental qubit results. Cases where quad-degeneracy exists are denoted as dotted rather than solid contours. Two traces (white, dashdotted lines) are also depicted. The first, at |Y | = 2.5Å (labelled A) is plotted in Figure 8. In contrast, the trace at |Y | = 3.6Å (labelled B) yields an equivalent eigenspectrum to the 1D case in Figure 2.
of |Y | = 2.5Å. If |X| is varied from an initial value of 1.6Å to 2Å, this extension spans phase space from one point where E 01 = 8 GHz to another. This separation is traced on Figure 7 (white dash-dotted line labelled A), where the defect visibly moves from a hemi-tetra-regime (ξ = 1/2) and goes through a transition region before reaching the tetra-well regime (ξ = 0). Figure 8 depicts the eigenspectrum of this trace. This is in contrast to larger confining pair values such as |Y | = 3.6Å (also traced on Figure 7: white dash-dotted line labelled B) where the eiginspectrum is largely unchanged from the one dimensional case presented in Figure 2.
As the |X| separation distance is increased, the effect on E 01 is negligible on the scale of the higher energy levels (although differs greatly on the TLS splitting energy scale). However, the degenerate pair E 23 rapidly shifts from a level much higher than ground to quasi-degenerate at ground. Complete quad-degeneracy is not always extant in the tera-well regime, configurations of |X| and |Y | in these domains usually have two quasi-degenerate pairs which are still observable, akin to the hemi-tetradomains, with the difference between the pairs approaching the difference of the pairs: E 01 = E 23 ≈ E 12 (1.95 |X| ≤ 2Å in Figure 8 for example).
FIG. 8. (Color online) Eigenspectrum of a 2D TLS: an oxygen atom caged with four aluminum atoms. Each excited state has been normalised with the ground state. The confinment pair |Y | is held at a constant distance separation of 2.5Å and the |X| range shows how higher eigenvalues behave as the confinement atoms force the defect into a tetra-well regime. Figure 7 depicts this data in terms of ξ (white dash-dotted trace labelled A).
When considering the influence of higher lying excited states, its important to keep the fundamental energy scales of the problem in mind. In typical qubit experiments, the superconducting properties of the device put a rigorous upper limit on TLS frequencies of interest. At frequencies greater than approximately 100 GHz, there is enough energy to dissociate Cooper-pairs and therefore it can be viewed as an operational upper bound for Josephson junction devices (typical operating frequencies are however device specific, and are much lower in practice). For much of the tetra-well domain E 12 ≥ 100 GHz and consequently can effectively be ignored, the system can be considered as a two-level defect even in this quad degeneracy domain.
The B type defects' sudden behaviour change as |X| → |Y | is also caused by this response (see Figure 7). Confinement pairs start interacting with the defect, causing E 23 to approach the value of E 01 (again, 1.95 |X| ≤ 2Å in Figure 8 depicts this phenomena). The ξ metric does not clearly differentiate between two quasidegenerate pairs that are marginally separated and two pairs that are actually degenerate.
If however, E 12 < 100 GHz, higher lying energy states must still be considered and the model exhibits true quad degenerate behaviour. Regions in which this occurs are denoted in Figure 7 as dotted contours, which over the entirety of phase space are extremely rare -which suggests a reason why quad-level systems are yet to be experimentally observed.
The final domain yet to be discussed on Figure 7 is the upper right hand corner where both |X| and |Y | are large. This region is tetra-well dominated but can be considered as a region where the TLS model breaks down. Each of the four potential minima exist localised about the four confining aluminium atoms and as such should not be considered as TLS candidates.

V. TLS DEFECT CONFINED IN THREE DIMENSIONS
To understand the properties of a delocalised oxygen, we have considered confining aluminium atoms in both a line and in a plane. In reality, aluminium atoms will surround the oxygen in all three dimensions. Two more confining atoms are therefore added into the system in the z direction labelled as |Z|. Interactions with these atoms in the third dimension are now considered, although the model is still two dimensional (i.e. 2+1D); thus oxygen continues to be confined to the xy plane. An illustration of this cluster configuration is displayed in Figure 9 and representations of the cage potential and oxygen wavefunctions are shown in Figures 10 and 11 for A and B type defects respectively.
The selection of a fixed |Z| distance changes the phase landscape in a manner that can be qualitatively extrapolated between two arbitrary values even a few angstroms apart. Values of |Z| = 2.75Å ( Figure 12) and |Z| = 2.25Å ( Figure 13) have been chosen to analyse in detail.
TLS behaviour can be observed on both maps and each value of |Z| has been selected based on model parameters. Oxygen confinement occurs for |Z| values lower than 2.25Å (i.e. E 01 100 GHz). |Z| values larger than 2.75Å show similar phase behaviour to that of Figure 7, which in completely unbound in z. Large |Z| separation distances also decrease the validity of the 2+1D model,  in addition: the radial distribution analysis in Ref. 17 suggests large separation distances for nearest neighbour atoms have a low probability of occurrence. As the pair separation distance |Z| decreases, the tetrawell (ξ = 0) regimes diminish in size and no longer exhibit TLS behaviour. This suggests that quad-degenerate defects, whilst quite rare in phase space already, are extremely rare in reality. For one to exist in a junction, the amorphous layer would have to be disordered in such a way that an oxygen atom's nearest neighbour atom pair exists at a distance greater that ∼ 3Å along one orthogonal basis vector. Whilst this configuration of six cage aluminium atoms on cardinal axes around a central oxygen is still an idealised system, Figure 13 is confined in all three spatial dimensions with pragmatic distances and is therefore considered to be the most 'realistic' representation of the TLS phase space for this model published here. TLS candidates lie well within hemi-tetra-(ξ = 1/2) domains, are not mired by higher energy level complexities (i.e. quad-quasi-degeneracies) and can be clearly identified as A type and B type defects separated by an (an)harmonic boundary.
Phase space is also dominated on this map with ξ = 5/4 (harmonic) and ξ = 2 (unique ground state) domains where the oxygen atom can be considered under the spatial harmonic approximation (i.e. not delocalised). This is significant, as TLS observations in experiments are not statistically dense. We expect few defects in the junction compared to the number of atoms extant.

VI. CHARGE DIPOLES
As stated in Section II, another important experimentally measurable property of the TLS is its strong electric dipole moment. Using the same |X|, |Y | and |Z| parameters from the phase maps in the previous section, the dipole elements in the x direction ℘ x , and the y direction ℘ y can be calculated via Equation 8. Figure 14 shows the same phase space as Figure 12, where |Z| = 2.75Å, and Figure 15 matches Figure 13, where |Z| = 2.25Å. The colourmap for both figures now represent dipole strengths of each cluster configuration rather than energy level splittings (represented through the ξ metric). These computed dipole moments correspond well to observed values, assuming O (nm) junction widths 8,13 .
The dipole moments are presented as (|℘ y | − |℘ x |) /e rather than separate plots because the dipole elements are discontinuous at the bifurcation points (i.e. when |℘ x | > 0, |℘ y | = 0 and vice versa). Comparing these plots to the phase maps in Section V, it is apparent that only the tetra-and hemi-tetra-domains (ξ < 1) exhibit a dipole response -which is appropriate for our model as E 01 splittings representative of a TLS only appear in these regions. Localised oxygen atoms (ξ > 1) are also expected to not elicit dipole behaviour.
With this information, the domain boundaries and bifurcations on each phase map can now be fully explained. Two variables alter the landscape: dipole and potential. Clusters with tight z confinement (those without tetra-well regions such as |Z| = 2.25Å) have four unique regions where a TLS may reside. The dipole direction dominates two of these domains: an A type region when the confining pair is collinear to the dipole, and a B type region when the confining pair is perpendicular. A symmetry bifurcation (at |X| = |Y |) separates the dipole domains into four regions which can therefore be clearly identified in terms of dipole moment (|℘ x | or |℘ y |) and defect type (A or B).
Clusters without as much z confinement (such as |Z| = 2.75Å) exhibit tetra-well behaviour: generating two additional regions. Tetra-well domains, as discussed in Subsection IV C, are caused when the confining pair of aluminium atoms start interacting with the defect pair (and hence the oxygen as well). If we consider the A type, |℘ y | domain in Figure 14, it is clear that the dominant dipole direction remains constant as |X| separation is in increased and the tetra-well domain is entered. The same |X|, |Y | parameters on Figure 15 cross a bifurcation line, changing dipole direction and the model indicates B type defect properties. Increased confinement in z induces a deeper potential well in the xy plane and removes any major landscape changing contributions from the confining atom pair -effectively reducing the model back to a 1D description. A cluster with a comparatively shallow potential that generates tetra-well domains does so when conditions are advantageous for an A type confinement pair to become a B type defect pair (a dipole direction change is not required for this to occur). A trace like the |Y | = 2.5Å (labelled A) on Figure 7 (and an associated eigenspectrum response similar to Figure 8) is an example of this behaviour. If we do not consider the small portion of quad-degenerate defects in this domain; the measurable properties (i.e. E 01 and dipole strength) tetra-well, B type, |℘ y | TLSs are identical to B type, |℘ y | TLSs that reside in a hemi-tetra-well domain.

VII. QUBIT COUPLING
To compare our TLS model directly to experiments, we assume that our JJ lies within a phase qubit, although the model applies equally for any device comprised of amorphous junctions. The measurable signal of a TLS in a phase qubit is the resonance of the TLS and qubit splitting energy, E 01 , with the qubit-TLS coupling. For the phase qubit 8 , the qubit-TLS coupling S max is a function of E 01 and ℘ 40 , the effective dipole moment due to an electric field applied in the direction of delocalization, Throughout this discussion we assume a junction width w = 2 nm and capacitance C = 850 fF. The dipole magnitudes in Figures 14 and 15 are calculated against the electron charge for simplicity, although as we are discussing an oxygen atom, the dipole elements may in fact be larger. Using our Josephson junction DFT models 17 we partition the charge density associated with atoms across the lattice into Bader volumes 41 . The charge enclosed within each Bader volume is a good approximation to the total electronic charge of an atom. An average value of 1.395 ± 0.006 e is found for oxygen atoms in a junction comprised of AlO 0.5 at a density of 0.8 times that of corundum (a common, low temperature and pressure phase of Al 2 O 3 ). We can use this value such that (for a dipole in the x direction) to gain a better estimate of S max . In Figure 16 we plot contour lines representing constant values of E 01 which correspond to the purview of experimentally observed qubit resonant frequencies for constant values of |Z| = 2.25Å. The S max (Equation 10) response to these frequencies is plotted as a function of |X|, in which we see maximum coupling strengths which correspond exceptionally well with experimental observations 5,6,13 . The most comprehensive of these studies Whilst the S max response is neither smooth nor singular over the phase space investigated, the value range is surprisingly small. Figure 17 shows the range of S max couplings for all E 01 = 8 GHz configurations calculated with confinements in |Z| from 2.25Å to unbound and |X| (hence |Y | as well from symmetry arguments) from 1.2 to 4Å. The entire range of S max values is only 60 MHz wide, which suggests an explanation as to why large couplings (of order 500 MHz for example) have not been seen experimentally.

VIII. STRAIN RESPONSE
Unusually long coherence times of strongly coupled defects are a key observation of TLS-qubit experiments 4,42 . As our model assumes a charge-neutral defect, coherence time is linked to the dipole element (for charge noise) and the strain response (for phonons). Mechanical deformation of a phase-qubit has recently been observed directly 43 , which we can mimic in our 2D model by introducing a series of deformations onto the cluster, and subsequently measure the E 01 response.
Whilst many deformation types were tested 17 , only two types (depicted in Figure 18a) show an active response over a length change (∆L) of 1 pm. Four clusters are chosen that lie on the 8 GHz contour when |Z| = 2.25Å, indicated as (b), (c), (d) and (e) in Figure 18a. As can be seen in Figure 15 due to its large, hyperbolic response -which is also seen in experiments 43 . The cluster configuration considered had a |Z| confinement of 2.5Å and the same |X|, |Y | separations as published here; thus it is not surprising that |X| translation yields an E 01 response is of O(10 4 ) GHz for |Z| = 2.25Å also. Therefore, translation of |X| has not been included in Figure 18 or the following discussion. Clusters (b) and (c) are both identified as A type defects and are relatively insensitive to |X| dilation and |Y | translation. Dilation of |Y | for these defects is a different way of saying "extending the defect pair separation" -which has been described in the above sections. The intensity of the response however is noticeably larger as the confinement distance (|X|) is increased.
B type defects, labelled (d) and (e), respond discordantly depending on their configuration. Dilation in |Y |, whilst a sizable strain source for A type configurations, initiates no response from the (e) configuration at all. B type defects located at this position in phase space have their defect pair in |X|, and point (e) specifically is confined with |Y | = 3.1Å. As discussed in the sections above, this separation distance is close to being practically unbound. Point (d) on the other hand has a tighter separation distance and begins to confine the system. Dilation in |X| responds in the opposite manner effectively: expanding defect separation distance whilst keeping the confinement distance static ultimately confines the wavefunction. The final response that the B type defects respond to is |Y | translation. As point (e) is essentially unbound in |Y | we do not expect any response from a strain of this magnitude. Point (d) on the other hand is interacting with its tighter confinement distance, with the TLS dipole collinear to the |X| direction. Translation of |Y | forces local minima  Figure 15). Translation of |X| is not shown on any graph as its' E01 response is of O(10 4 ) over (∆L) for all cases.
of the potential landscape from points on the x axis to locations off axis, yet still within the minima rings about the aluminium atoms. In other words, the wavefunction is slightly rotated around the defect axis.

IX. CONCLUSIONS AND OUTLOOK
Even with the extreme idealisation of this model, it allows the prediction of experimentally measured properties of strongly coupled TLSs with atomic positions as the only input parameters and therefore shows as a proof of concept that these defects can arise in AlO x without any alien species present.
Our model shows that even when only considering delocalisation in two dimensions, a range of different behaviour can be seen. The existence of effective twostate systems can come about through different potential shapes, which in turn arise due to various atomic position configurations. To understand these different configurations, we consider both the symmetries of the eigenspectrum and the effective charge dipole of the defect.
We find that two-level systems with equivalent properties to those seen in experiments can be formed from atomic configurations which are entirely consistent with the material properties of amorphous aluminium oxide barriers. Most interestingly we find that the expected qubit-TLS coupling strength and the TLS strain response correspond very well to that observed experimentally.
A complication that occurs with phenomenological models that attempt to describe TLS behaviour is their free parameter count is high enough to describe all facets of the observed experimental parameters without being distinguishable from other, non-complimentary models 13 . Whilst this delocalised oxygen model only uses one type of parameter as input (atomic positions), it still requires an x, y, z coordinate set for up to 6 cage atoms. More sophisticated modelling techniques such as molecular mechanics and density functional methods can be then used to obtain more realistic values for the atomic positions 17 . Microscopic models of this type will guide future fabrication and design of superconducting circuits, leading to lower levels of noise and greater control over their quantum properties.