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

Superconducting qubits and Josephson junction (JJ) 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–8], although little is known about their true microscopic nature. Many phenomenological theories attempting to describe them exist; including charge dipoles [9], Andreev bound states [10], magnetic dipoles [11], Kondo impurities [12] and TLS state dependence of the JJ transparency [13]. Detailed fitting of experimental data can place limits on these models [14] 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 [15] and surface aluminium ions paramagnetically coupling to ambient molecules [16] 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 [17]. Amorphous insulating barriers (either in the form of JJs 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, 9]. 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 [18], we consider the origin of some defects to be within the amorphous oxide layer itself, specifically an oxygen atom in a spatially delocalized state. This has important ramifications for materials science based efforts to reduce the effects of TLSs. If, as alternative models suggest, TLS defects indeed stem from surface states [19] or the accidental inclusion of an alien species [20, 21]; future fabrication techniques or more robust qubit designs may suppress or diminish the response of such noise sources as has been the case historically [9, 22–25]. 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 JJ 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 delocalizes around the Si–Si bond axis [26]. This forms an anharmonic system with a quasi-degenerate [18] 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 idealized trigonal-like two-dimension (2D) oxide lattice. Possible defects of this nature are depicted in figure 1.

Zoom In Zoom Out Reset image size

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 2 introduces the theory from which we build our model to investigate delocalized oxygen based TLSs. 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 systematic way. Section 3 considers a defect comprising of an Al–O–Al chain in one-dimension (1D), perturbed from a crystalline lattice, simulating defects A and B in figure 1. In reality, an oxygen atom in the amorphous layer of a JJ will be surrounded by atoms in all three dimensions. Moving towards a model representation of this, section 4 extends the model to 2D, with four aluminium atoms confining an oxygen in a plane. The completed model is then described in section 5, which extends oxygen confinement into three dimensions (with six aluminum atoms), whilst still projecting oxygen delocalization on a plane. Although in general an oxygen can delocalize in three dimensions, for this investigation we focus on an effective $2+1$ D model, minimizing both computational and descriptive complexity while still modelling the relevant behaviour of the system. The following sections then apply the $2+1$ D model and compare delocalized oxygen responses to experimental TLS data. Section 7 discusses qubit–TLS coupling and section 8 observes the effect of mechanical strain.

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

$$\begin{eqnarray}&&H=-\frac{{{\hbar }^{2}}}{2{{m}_{{\rm oxy}}}}{{\nabla }^{2}}+V({\bf r}),\end{eqnarray} \tag{ 1 }$$

where m_oxy is the mass of an oxygen atom and $V({\bf 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₀₁ and a strong electric dipole moment. Observed values of E₀₁ differ between qubit architectures: approximately 1–10 GHz for transmons [23], 4–5 GHz for flux qubits [6] (although recent designs can tune this gap down to the MHz range [33]) and nominally 7–8 GHz for phase qubits [14]. Dipole moment strengths also vary, but are usually on the order of 1 eÅ [5, 14]. Whilst many other properties have been measured, this work will focus on obtaining respectable values for E₀₁ 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 will be used, to compare directly with TLSs studied in the phase qubits of [14].

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 as 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 [34], 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 fixed-charge models [35, 36] due to the complex geometry of the JJ. A variable charge potential such as Streitz–Mintmire can capture the variable oxygen states when present in a predominantly metallic environment through the minimization of the electrostatic term in the potential (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

$$\begin{eqnarray}&&V({\bf r})={{E}_{{\rm EAM}}}+\mathop{\sum }\limits_{i}^{N}{{q}_{i}}{{\chi }_{i}}+\frac{1}{2}\mathop{\sum }\limits_{i,j}^{N}{{q}_{i}}{{q}_{j}}{{V}_{ij}},\end{eqnarray} \tag{ 2 }$$

where ${{q}_{\{i,j\}}}$ is the atomic charge, ${{\chi }_{i}}=\chi _{i}^{0}+{{\sum }_{j}}{{Z}_{j}}\left( \left[ j\mid {{f}_{i}} \right]-\left[ {{f}_{i}}\mid {{f}_{j}} \right] \right)$ and ${{V}_{ij}}=J_{i}^{0}{{\delta }_{ij}}+\left[ {{f}_{i}}\mid {{f}_{j}} \right]$. $\chi _{i}^{0}$ is the electronegativity and J_i⁰ is an empirical parameter known as 'atomic hardness' or a self-Coulomb repulsion [37] of atom i. Z_j is an effective core charge of atom j, ${{\delta }_{ij}}=1$ when i = j and ${{\delta }_{ij}}=0$ when $i\ne 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 (e.g. j) and effective core charge densities (e.g. f_i) and takes the form [38]

$$\begin{eqnarray}&&\left[ a\mid {{f}_{b}} \right]=\int \frac{{{f}_{b}}({{r}_{b}},{{q}_{b}})}{{{r}_{av}}}{\rm d}{{V}_{b}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\left[ {{f}_{a}}\mid {{f}_{b}} \right]=\int \int \frac{{{f}_{a}}({{r}_{a}},{{q}_{a}}){{f}_{b}}({{r}_{b}},{{q}_{b}})}{{{r}_{vv}}}{\rm d}{{V}_{a}}{\rm d}{{V}_{b}}\end{eqnarray} \tag{ 4 }$$

with $a=i,j;b=i,j;a\ne b$ in (2) to prohibit self interactions. ${\rm d}{{V}_{a,b}}$ are integrating volume units. r_av is therefore the center distance between atom a and ${\rm d}{{V}_{b}}$, and r_vv is the center distance between ${\rm d}{{V}_{a}}$ and ${\rm d}{{V}_{b}}$.

The first term in (2), E_EAM, does not depend on the partial charges q_i and therefore describes a charge-neutral system, represented here with a quantum mechanical based empirical embedded atom model (EAM) for the Al–Al and Al–O interactions

$$\begin{eqnarray}&&{{E}_{{\rm EAM}}}=\mathop{\sum }\limits_{i}^{N}{{F}_{i}}\left[ {{\rho }_{i}} \right]+\mathop{\sum }\limits_{i\lt j}^{N}{{\phi }_{ij}}({{r}_{ij}}),\end{eqnarray} \tag{ 5 }$$

with ${{F}_{i}}\left[ {{\rho }_{i}} \right]$ as the energy required to embed atom i in a local electron density ${{\rho }_{i}}$, and ${{\phi }_{ij}}({{r}_{ij}})$ describing the residual pair–pair interactions by way of Buckingham and Rydberg potentials

$$\begin{eqnarray}&&{{\phi }_{ij}}({{r}_{ij}})=A{\rm exp} \left( -\frac{{{r}_{ij}}}{\rho } \right)-B\left[ 1+C\left( \frac{{{r}_{ij}}}{{{r}_{0}}}-1 \right) \right]{\rm exp} \left[ -C\left( \frac{{{r}_{ij}}}{{{r}_{0}}}-1 \right) \right],\end{eqnarray} \tag{ 6 }$$

where r_ij is the interatomic (Euclidean) distance between atoms i and j, all other constants are listed in table 1. Further formalisms and parameters can be found in [34, 38], implementation is also discussed in [39].

Using this potential, a single body time-independent Hamiltonian is constructed using a 7-point central difference method

$$\begin{eqnarray}&&{{f}^{\prime\prime }}({{x}_{0}})=\frac{2{{f}_{-3}}-27{{f}_{-2}}+270{{f}_{-1}}-490{{f}_{0}}+270{{f}_{1}}-27{{f}_{2}}+2{{f}_{3}}}{180{{h}^{2}}}+\mathcal{O}({{h}^{6}}),\end{eqnarray} \tag{ 7 }$$

where ${{f}_{k}}=f\left( {{x}_{0}}+kh \right)$. The step size h = 0.1 Å was found to be optimal for the grid range relevant to our purposes [40]. A sparse matrix representation of H (1) is then solved via direct diagonalization, obtaining eigenvectors and eigenenergies with precision better than 10 kHz. In contrast, the energy scale for JJ defects observed in experiments is $\lesssim 10$ GHz [4, 6, 7].

The dipole element is computed using numerical integration of the ground- and first-excited states (${{\psi }_{0}}$, ${{\psi }_{1}}$), where

$$\begin{eqnarray}&&{{\wp }_{x}}=\int \int \psi _{0}^{*}(x,y)x{{\psi }_{1}}(x,y)\;{\rm d}x{\rm d}y\end{eqnarray} \tag{ 8 }$$

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₀) to first excited state (E₁) energy splitting is defined as E₀₁. For comparison with experimental results, energy is expressed in frequency units throughout this paper.

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 (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 $\alpha$-${\rm A}{{{\rm l}}_{2}}{{{\rm O}}_{3}}$) has an Al–O bond distance of $\sim 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 $\mid X\mid$ or simply the 'defect pair'. Displacing $\mid X\mid$ 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 seven lowest energy levels of this system over a continuum of values in $\mid X\mid$ is depicted in figure 2. Each energy is measured relative to the ground state, which shows two particular regions where E₀₁ (green, solid line) 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 $\mid X\mid \sim 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 localized.

Zoom In Zoom Out Reset image size

To investigate the potential landscape and the resultant oxygen wavefunction in sections A and B of figure 2, we choose two separation distances: $\mid X\mid =1.5$ Å (figure 3), which lies in the A type defect section, and $\mid X\mid =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 delocalization of the oxygen in 2D so that both defect types can be identified on a continuum separation in $\mid X\mid$ (as figure 2 depicts) rather than using two separate coordinate systems.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 illustrated 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.

The ideal case discussed in section 3 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 $\mid Y\mid$ direction (i.e. $y=-Y,+Y$).

4.1. Classifying eigenspectrum dynamics

As the values of $\mid X\mid$ or $\mid Y\mid$ 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

$$\begin{eqnarray}&&\xi =\frac{{{E}_{1}}-{{E}_{0}}}{{{E}_{2}}-{{E}_{0}}}+\frac{{{E}_{3}}-{{E}_{0}}}{{{E}_{4}}-{{E}_{0}}}.\end{eqnarray} \tag{ 9 }$$

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 [18], where ξ is plotted as a function of the distance to the confining aluminium atoms ($\mid X\mid ,\mid Y\mid$). 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 6). 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. It is clear from the outset potential projections of figures 3 and 4 that both A and B type defects reside in hemi-tetra-wells.

Zoom In Zoom Out Reset image size

The ξ metric is capable of identifying the tetra- ($\xi =0$) and hemi-tetra- ($\xi =1/2$) domains, as well as harmonic ($\xi =5/4$) and unique ground state (rotationally symmetric, Mexican hat-like) ($\xi =2$) regimes and finally the location of bifurcations or transitions ($\xi =1$). Figure 6 shows the corresponding layouts of each interplay. It is worth noting that $\xi =3/2$ can also be considered harmonic for the lowest three levels.

Zoom In Zoom Out Reset image size

4.2. Visualizing phase space and identifying TLSs

Using this ξ metric and varying the values of both $\mid X\mid$ and $\mid Y\mid$ generates figure 7, a phase map of the interplay of low energy states of the oxygen atom confined in 2D by aluminium atoms. The x and y axes show the $\mid X\mid$ and $\mid Y\mid$ pair separation distances respectively over a range of 1.85–4 Å and the phase colour indicates which ξ region a particular configuration exists in.

Zoom In Zoom Out Reset image size

Regions that are deemed to be unimportant when searching for TLS behaviour are those where $\xi \geqslant 1$, as these correspond to eigenspectra which do not possess a (quasi-)doubly degenerate ground state.

As discussed in section 3, in the 1D confinement case, both A and B type defects exist in a hemi-tetra-well ($\xi =1/2$). This is also the case for the 2D 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₀₁ splittings in various qubit architectures described in section 2. For example: anywhere an orange, 8 GHz line exists on this phase map, the $\mid X\mid ,\mid Y\mid$ coordinates of the model generate a cluster configuration that yields the same E₀₁ observed in phase qubit experiments [14].

Many of these lines lie completely within a $\xi =1/2$ (green) region. Consider the contour set on the left of figure 7 where $\mid X\mid \simeq 1.5$ Å: section 3 states this distance is indicative of an A type defect. The addition of the $\mid Y\mid$ aluminium pair does little to perturb the potential at larger distances ($\mid Y\mid \gtrsim 3$ Å) and causes only slight deformations at smaller distances ($2\lesssim \mid Y\mid \lesssim 3$ Å). As the phase space is symmetric about $\mid X\mid =\mid Y\mid$, A type defects exist at the bottom of the plot where $\mid Y\mid \simeq 1.5$ Å as well.

B type defects also exist in a $\xi =1/2$ region, when $\mid X\mid$ or $\mid Y\mid \simeq 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.

4.3. Analysis of phase space domains

The case of $\xi =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 $\xi =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 $\mid X\mid$ and a confining pair at a constant value of $\mid Y\mid =2.5$ Å. If  $\mid X\mid$ is varied from an initial value of 1.6–2 Å, this extension spans phase space from one point where ${{E}_{01}}=8$ GHz to another. This separation is traced on figure 7 (white dashed–dotted line labelled A), where the defect visibly moves from a hemi-tetra- regime ($\xi =1/2$) and goes through a transition region before reaching the tetra-well regime ($\xi =0$). Figure 8 depicts the eigenspectrum of this trace. This is in contrast to larger confining pair values such as $\mid Y\mid =3.6$ Å (also traced on figure 7: white dashed–dotted line labelled B) where the eigenspectrum is largely unchanged from the 1D case presented in figure 2.

Zoom In Zoom Out Reset image size

As the $\mid X\mid$ separation distance is increased, the effect on E₀₁ is negligible on the scale of the higher energy levels (although differs greatly on the TLS splitting energy scale). However, the degenerate pair E₂₃ (a degenerate second and third excited state pair) 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 $\mid X\mid$ and $\mid Y\mid$ in these domains usually have two quasi-degenerate pairs which are still observable, akin to the hemi-tetra- domains, with the difference between the pairs approaching the difference of the pairs: ${{E}_{01}}={{E}_{23}}\approx {{E}_{12}}$ ($1.95\lesssim \mid X\mid \leqslant 2$ Å in figure 8 for example). In other words, when the ground to first excited and second to third excited state differences are equal, the first to second excited state difference trends from a much larger value to one effectively equivalent to the aforementioned pairs in this region.

When considering the influence of higher lying excited states, it is 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 JJ devices (typical operating frequencies are however device specific, and are much lower in practice). For much of the tetra-well domain ${{E}_{12}}\geqslant 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 $\mid X\mid \to \mid Y\mid$ is also caused by this response (see figure 7). Confinement pairs start interacting with the defect, causing E₂₃ to approach the value of E₀₁ (again, $1.95\lesssim \mid X\mid \leqslant 2$ Å in figure 8 depicts this phenomena). The ξ metric does not clearly differentiate between two quasi-degenerate pairs that are marginally separated and two pairs that are actually degenerate.

If however, ${{E}_{12}}\lt 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 $\mid X\mid$ and $\mid Y\mid$ 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 localized about the four confining aluminium atoms and as such should not be considered as TLS candidates.

To understand the properties of a delocalized 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 $\mid Z\mid$. Interactions with these atoms in the third dimension are now considered, although the model is still 2D (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 figure 10 for A (left) and B (right) type defects.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The selection of a fixed $\mid Z\mid$ distance changes the phase landscape in a manner that can be qualitatively extrapolated between two arbitrary values even a few angstroms apart. Values of $\mid Z\mid =2.75$ Å (figure 11) and $\mid Z\mid =2.25$ Å (figure 12) have been chosen to analyze in detail.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

TLS behaviour can be observed on both maps and each value of $\mid Z\mid$ has been selected based on model parameters. Oxygen confinement occurs for $\mid Z\mid$ values lower than 2.25 Å (i.e. ${{E}_{01}}\gg 100$ GHz). $\mid Z\mid$ values larger than 2.75 Å show similar phase behaviour to that of figure 7, which in completely unbound in z. Large $\mid Z\mid$ separation distances also decrease the validity of the 2+1D model, in addition: the radial distribution analysis in [18] suggests large separation distances for nearest neighbour atoms have a low probability of occurrence.

As the pair separation distance $\mid Z\mid$ decreases, the tetra-well ($\xi =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 idealized system, figure 12 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- ($\xi =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 $\xi =5/4$ (harmonic) and $\xi =2$ (unique ground state) domains where the oxygen atom can be considered under the spatial harmonic approximation (i.e. localized). 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.

As stated in section 2, another important experimentally measurable property of the TLS is its strong electric dipole moment. Using the same $\mid X\mid ,\mid Y\mid$ and $\mid Z\mid$ 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 13 shows the same phase space as figure 11, where $\mid Z\mid =2.75$ Å, and figure 14 matches figure 12, where $\mid Z\mid =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 $\mathcal{O}\left( {\rm nm} \right)$ junction widths [9, 14].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The dipole moments are presented as $\left( \mid {{\wp }_{y}}\mid -\mid {{\wp }_{x}}\mid \right)/e$ rather than separate plots because the dipole elements are discontinuous at the bifurcation points (i.e. when $\mid {{\wp }_{x}}\mid \gt 0,\mid {{\wp }_{y}}\mid =0$ and vice versa). Comparing these plots to the phase maps in section 5, it is apparent that only the tetra- and hemi-tetra- domains ($\xi \lt 1$) exhibit a dipole response—which is appropriate for our model as E₀₁ splittings representative of a TLS only appear in these regions. Localized oxygen atoms ($\xi \gt 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 $\mid Z\mid =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 $\mid X\mid =\mid Y\mid$) separates the dipole domains into four regions which can therefore be clearly identified in terms of dipole moment ($\mid {{\wp }_{x}}\mid$ or $\mid {{\wp }_{y}}\mid$) and defect type (A or B).

Clusters without as much z confinement (such as $\mid Z\mid =2.75$ Å) exhibit tetra-well behaviour: generating two additional regions. Tetra-well domains, as discussed in section 4.3, 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, $\mid {{\wp }_{y}}\mid$ domain in figure 13, it is clear that the dominant dipole direction remains constant as $\mid X\mid$ separation is in increased and the tetra-well domain is entered. The same $\mid X\mid ,\mid Y\mid$ parameters on figure 14 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 $\mid Y\mid =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₀₁ and dipole strength) tetra-well, B type, $\mid {{\wp }_{y}}\mid$ TLSs are identical to B type, $\mid {{\wp }_{y}}\mid$ TLSs that reside in a hemi-tetra-well domain.

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₀₁, with the qubit–TLS coupling. For the phase qubit [9], the qubit–TLS coupling S_max is a function of E₀₁ and ℘ [41], the effective dipole moment due to an electric field applied in the direction of delocalization

$$\begin{eqnarray}&&{{S}_{{\rm max} }}=2\frac{\wp }{w}\sqrt{\frac{{{e}^{2}}}{2C}{{E}_{01}}}.\end{eqnarray} \tag{ 10 }$$

Throughout this discussion we assume a junction width w = 2 nm and capacitance C = 850 fF.

The dipole magnitudes in figures 13 and 14 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 JJ DFT models [18] we partition the charge density associated with atoms across the lattice into Bader volumes [42]. 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\pm 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 ${\rm A}{{{\rm l}}_{2}}{{{\rm O}}_{3}}$). We can use this value such that

$$\begin{eqnarray}&&\tilde{{{\wp }_{x}}}=1.4e{{\wp }_{x}}\end{eqnarray} \tag{ 11 }$$

(for a dipole in the x direction) to gain a better estimate of S_max.

In figure 15 we plot contour lines representing constant values of E₀₁ which correspond to the purview of experimentally observed qubit resonant frequencies for constant values of $\mid Z\mid =2.25$ Å. The S_max (10) response to these frequencies is plotted as a function of $\mid X\mid$, in which we see maximum coupling strengths which correspond exceptionally well with experimental observations [5, 6, 14]. The most comprehensive of these studies [5] measures S_max values of 3–45 MHz, assuming a 1 Å dipole moment.

Zoom In Zoom Out Reset image size

Whilst the S_max response is neither smooth nor singular over the phase space investigated, the value range is surprisingly small. Figure 16 shows the range of S_max couplings for all ${{E}_{01}}=8$ GHz configurations calculated with confinements in $\mid Z\mid$ from 2.25 Å to unbound and $\mid X\mid$ (hence $\mid Y\mid$ 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.

Zoom In Zoom Out Reset image size

Unusually long coherence times of strongly coupled defects are a key observation of TLS–qubit experiments [4, 43]. 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 [44], which we can mimic in our 2D model by introducing a series of deformations onto the cluster, and subsequently measure the E₀₁ response.

Whilst many deformation types were tested [18], only two types (depicted in figure 17(a) show an active response over a length change ($\Delta L$) of 1 pm. Four clusters are chosen that lie on the 8 GHz contour when $\mid Z\mid =2.25$ Å, indicated as (b), (c), (d) and (e) in figure 17(a). As can be seen in figure 14, each of these cases are within a $\mid {{\wp }_{x}}\mid$ dominant domain.

Zoom In Zoom Out Reset image size

The investigation in [18] focused on $\mid X\mid$ translation due to its large, hyperbolic response—which is also seen in experiments [44]. The cluster configuration considered had a $\mid Z\mid$ confinement of 2.5 Å and the same $\mid X\mid$, $\mid Y\mid$ separations as published here; thus it is not surprising that $\mid X\mid$ translation yields an E₀₁ response is of $\sim {{10}^{4}}$ GHz for $\mid Z\mid =2.25$ Å also. Therefore, translation of $\mid X\mid$ has not been included in figure 17 or the following discussion.

Clusters (b) and (c) are both identified as A type defects and are relatively insensitive to $\mid X\mid$ dilation and $\mid Y\mid$ translation. Dilation of $\mid Y\mid$ 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 ($\mid X\mid$) is increased.

B type defects, labelled (d) and (e), respond discordantly depending on their configuration. Dilation in $\mid Y\mid$, 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 $\mid X\mid$, and point (e) specifically is confined with $\mid Y\mid =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 $\mid X\mid$ 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 $\mid Y\mid$ translation. As point (e) is essentially unbound in $\mid Y\mid$ 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 $\mid X\mid$ direction. Translation of $\mid Y\mid$ forces local minima 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.

Even with the extreme idealization 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 delocalization in 2D, a range of different behaviour can be seen. The existence of effective two-state 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 TLSs 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 [14]. Whilst this delocalized oxygen model only uses one type of parameter as input (atomic positions), it still requires an $x,y,z$ coordinate set for up to six 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 [18]. 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.

This research was supported under the Australian Research Council's Discovery Projects funding scheme (project number DP140100375). Computational resources were provided at the NCI National Facility systems at the Australian National University through the National Computational Merit Allocation Scheme supported by the Australian Government.