The influence of pinholes and weak-points in aluminium-oxide Josephson junctions

Josephson junctions are one of the key components used in superconducting qubits for quantum computers [1, 2, 3]. Whilst there have been significant developments in the fabrication processes of Josephson junctions which now enable us to have working examples of small scale quantum computing [4, 5, 6], we still see significant problems with device stability and reproducibility. In particular, large variabilities have been found in the coherence times with time [7, 8], and across different qubits on a quantum computing chip [9, 10, 11]. Variations have also been found in the critical currents between Josephson junctions fabricated in the same way [12]. These variations in the qubit parameters means that significant tuning is required to ensure qubits across the chip function similarly. Although qubits can be tuned [13, 14], this typically requires more gates to be patterned onto quantum computing chips, as well as recalibration during the computation [6, 15]. This significantly limits the reliability of quantum computing chips that are fabricated today. It is not yet clear what causes stability and reproducibility problems in these devices, however a large amount of research discusses the importance of the tunnel barrier in this context [16]. Often comprised of an amorphous metal-oxide material, the tunnel barrier is inherently disordered, inviting the possibility of structures that could act as two-level systems that couple to the qubit current causing instabilities in its parameters. A particular type of imperfection in the tunnel barrier that has been proposed are microscopic metallic links, coined “pinholes”. These pinholes may increase the current that couples to two-level systems in the barrier, exacerbating instabilities [17].

In the literature there is currently conflicting evidence on the role and prevalence of pinholes. Zhou et al. [18] finds that pinholes may exist in particular circumstances, such as when the oxide barrier is too thin, whereas Greibe et al. [19] develops a model which rules out pinholes as the source of excess current. On the other hand, Tolpygo et al. [20] shows that pinholes may not exist in a newly fabricated device, however they may form through degradation of the oxide when the junction is used in a circuit. Recent work has studied the effect of Josephson harmonics on the current-phase relationship (CPR), highlighting that the disordered nature of the oxide barrier could contribute to effects in the CPR, specifically if there are high transmission channels present [21].

This paper uses computational techniques to study how pinholes form and aims to investigate what effects they might have on transport properties in Josephson junction devices that are fabricated today. Commonly, Josephson junctions are fabricated as Al$-$AlO${}_{x}-$Al tunnel junctions [22, 23, 24]. This study investigates pinhole formation during AlO_x growth to discuss under which conditions they might be stable, and what effect stable pinholes might have on electronic transport through the device. We begin by discussing the stability of pinholes following Zhou et al. [18] which models a cylindrical pinhole in an otherwise uniform oxide barrier. Using the interfacial surface energies between materials in the system the stability of pinholes can be determined. Using a similar method we determine the stability of cylindrical pinholes with increasing radii for different barrier thicknesses. Using the same geometries we then calculate the transport properties for each of these devices.

To then investigate the formation and closing of pinholes more realistically, we use our previously developed molecular dynamics approach [25] to grow Al$-$AlO${}_{x}-$Al tunnel junctions with different amounts of oxide in the barrier.

Using a tight-binding description of the system we employ a non-equilibrium Green’s function (NEGF) approach to investigate normal (single-electron) transport in the structures studied herein. NEGF is a numerical approach which has been widely used for calculating transport properties in nanoscale devices [26], and has also been used explicitly for the case of Al$-$AlO${}_{x}-$Al tunnel junctions [27]. In this work, we determine the transmission, resistance-area, and the current density through full three-dimensional models of Al$-$AlO${}_{x}-$Al tunnel junctions with different amounts of oxide. Understanding variations in the normal transport based on structural irregularities can suggest causes for variation in the superconducting properties.

To understand how the normal transport translates to the superconducting (correlated-electron) transport, the Ambegaokaar-Baratoff (AB) relation is commonly used [28, 29]. It specifically relates normal-state resistance to the superconducting critical current. The validity of this relationship is unclear when there are pinholes present in the oxide barrier. As an alternative the model in this work is used to calculate the transmission eigenvalues for each of our structures, which in turn allows us to determine the contributions from different transmission channels and the effects they have on the CPR.

From a toy model we determine that pinholes may be present in thinner barriers and a source of single-electron current. We then use the analysis from the toy model to understand the transport properties of a three-dimensional atomistic model describing Al$-$AlO${}_{x}-$Al tunnel junctions derived from molecular dynamics simulations. Our results show there is a significant amount of variability in the transport properties. These variabilities are still present for devices without pinholes, indicating that the variation is a product of the inherently disordered AlO_x barrier. We also find that the amorphous structure of the oxide results in a non-uniform electric potential, which makes the metal-insulator transition of the barrier less clear compared to analysis from our toy model. Specifically, some results highlight that thick barriers with no pinholes may remain conductive. These devices do not have clear pinholes in the barrier, but have weak-points which can facilitate single-electron current flow. This could indicate a situation in which a seemingly complete barrier could result in a dysfunctional device.

Given the complex structure of disordered AlO_x, we first develop an idealised model to understand how pinholes of a simple structure remain stable in the barrier and also how systematic changes in pinhole size may affect the transport properties of tunnel junctions. For this section we use a geometric model with pinholes described as cylinders in an otherwise uniform oxide layer in Al$-$AlO${}_{x}-$Al tunnel junctions. Fig. 1 schematically depicts the geometry of our model.

Having developed an atomistic model of Al$-$AlO${}_{x}-$Al tunnel junctions using molecular dynamics we now use this model to obtain information about the transport properties of these 3D simulated junctions. The transport calculations in this section do no include superconducting parameters or considerations of Cooper pairs of electrons. They describe the device in its normal state and are purely single-electron transport calculations.

We use a NEGF method to calculate the transmission and the normal resistance, $R_{\textrm{N}}$, in the same way as we did for the toy model in Sec. II.2. Fig. 9 shows the resultant $R_{\textrm{N}}A$ relationship for junctions grown on each substrate. As we expect $R_{\textrm{N}}A$ to be exponential with film thickness [57], results are plotted on a logarithmic scale.

On each of the data sets in Fig. 9 the dashed line indicates the resistance quantum that corresponds to $T=1$ as determined in Fig. 3. It corresponds to the expected transition from metallic to insulating behaviour. For every data set we see significant variation in the resistance values, and no clear transition to insulating behaviour until  400 oxygen atoms are deposited. An exception to this is in the ‘(111)1’ data set where we find that the last data point is below the metal-to-insulator transition even with approximately 400 oxygen atoms in the barrier, like the other structures. This implies there may be more metallic regions present in the barrier in the form of pinholes or weak-points.

In this study we aim to create a distinction between pinholes and weak-points. We define pinholes to be present when the Al substrate is $<$100% covered in oxide, which results in a direct metallic link once metal is deposited as a top contact (as in Sec. II.2). We define weak-points to be present when our surface coverage calculations show no presence of pinholes, however the transport calculations indicate a localised higher current density compared to a uniform barrier. Such localised variation in the current density results from variation in the oxide density or stoichiometry but may not correspond to the conventional image of a pinhole, as illustrated in Fig. 1.

Comparing these results to Fig. 8, we find that an increase in surface coverage (especially close to 100%) isn’t directly correlated to insulating behaviour in the barrier. We reach close to 100% coverage when approximately 200 oxygen atoms are deposited, however for all data sets in Fig. 9 we see that the $R_{\textrm{N}}A$ data points around 200 oxygen atoms are close to or below the metal-to-insulator transition, indicating metallic behaviour. This highlights that surface coverage alone cannot give us a good indication of a functioning barrier.

The variation we see in the resistance data in Fig. 9 is consistent with measurements in manufactured Al$-$AlO${}_{x}-$Al tunnel junctions [9, 10, 11]. For example, it has been shown that when some device barriers are grown to thicknesses which are typically expected to display insulating behaviour, a few may fail due to an excess current [20]. Dataset ‘(111)1’ shows a theoretical example which could explain these situations. As well as this, experiments have shown significant variability in the transport between devices fabricated in the same way. We believe these variabilities may arise from the disordered nature of the AlO_x region as results in Fig. 3 showed significantly less variability, likely due to the uniform potential barrier used and systematic changes to the pinhole size. It has also previously been shown that crystal grain boundaries in the aluminium bottom contact could cause thickness variations in the AlO_x layer that is grown above [52]. These thickness variations could then in turn cause irregularities in the transport measurements. In this study we have not included the effects of grain boundaries, we grow devices on single crystal aluminium substrates. Despite omitting grain boundaries, we still see significant variability, highlighting that disorder alone may be a key contributor to the variation in transport we see in junctions of this sort.

There also appears to be a larger spread in the data for junctions grown on Al (100) substrates compared to Al (111). This could indicate that the spatial variation of the junction resistance for the (111) crystal orientation is lower, however further statistics are required to be definitive. The computational expense of growing large atomic structures and performing atomic-resolution transport calculations limits our ability to perform detailed replication studies at this time.

In prior work by Cyster et al. [27] a reliable exponential trend was observed in $R_{\textrm{N}}A$ calculations for junctions with increasing thickness. That work used a melt and quench method to create devices which gave control over the oxide thickness, however it did not consider the experimental fabrication processes used. In this study we mimic the fabrication process more closely [25] aiming to provide a more accurate representation of what would be obtained from measurements of real devices. However, this does produce greater variability in the transport results. Cyster et al. [27] also uses a S-M potential, whereas this work uses a ReaxFF potential which has since been shown to more accurately describe the self-limiting behaviour we see experimentally [25].

To investigate the variation we see in the transport calculations, Fig. 10 shows the transmission function for each of the thickest structures, corresponding to more than 350 oxygen atoms in the simulation cell.

In particular, Fig. 10 shows the sensitivity of transport measurements to changes in the Fermi energy, noting that our calculations all use $E_{\textrm{F}}=0.7$ eV, marked on Fig. 10 with a dashed line. From the transmission we also see that the ‘(111)1’ structure is more transparent than the other three structures, supporting results from Fig. 9 where we found that the thickest structure was metallic compared to the junctions in the other data sets.

Our calculations discussed up to this point give us an indication of the overall transport in the junction. However to understand the specific effects that pinholes have on transport within the barrier we need to study how current flows through the junction spatially. To investigate this we calculate the three-dimensional current density throughout the junction.

The current flowing between two points in a junction is given by the element-wise product of $H$ and $G^{n}$,

where the electron Green’s function is given by

and

Using $G^{n}(E)$, the charge density can be calculated in three-dimensions with

The current in one direction is calculated from the difference between pairs of points, for the $z$ direction this looks like

Fig. 11 shows $yz$ slices of the barrier and the charge and current densities for positions in those slices.

Dark regions in Fig. 11 correspond to areas in the junction which have a higher charge density, typically indicative of a metallic region. The orange arrows show the single-electron current density at positions in the junction. In both figures it is apparent that the largest proportion of current flows through localised regions of the junction. Fig. 11 (a) shows the lowest surface coverage junction in the ‘(111)1’ data set (26 oxygen atoms), and Fig. 11 (b) shows the same substrate but for the thickest oxide for full surface coverage (397 oxygen atoms). The current in Fig. 11 (a) is on average approximately four times higher than in Fig. 11 (b) as is expected for a more transparent barrier. Of particular interest, despite the device in Fig. 11 (b) having a thick enough barrier to expect insulating behaviour, there are still weak-points which facilitate peaks in the single-electron current flow.

The ‘(111)1’ structure shown in Fig. 11 (b) is the same device that has displayed a higher transparency in Fig. 10 and Fig. 9 (c). Fig. 11 allows us to probe how the microscopic structure of the device may contribute to the higher transparency across the junction, emphasising that the disordered barrier may be a crucial consideration when trying to predict the parameters of fabricated Josephson junction devices.

In general, defects in a superconducting qubit can couple via charge, flux, or critical current in the circuit [58]. In the case that defects in the junction modulate the critical current [17, 16], the qubit-defect coupling strength is proportional to the magnitude of the (normal or super) current flow at the site of the defect. Therefore, increased variance in the current density through the junction can result in increased coupling to defects.

Fig. 12 shows $xy$ slices of the magnitude of the current density perpendicular to the transport direction for the thickest barriers in ‘(100)1’ and ‘(111)2’ data sets. For figures that follow we compare the ‘(100)1’ and ‘(111)2’ data sets where the thickest two barriers are both insulating.

The slices shown in Fig. 12 are taken at equal intervals through the oxygen extent of the barrier region. The distribution of the current density between slices varies greatly, specifically we see that the current peaks in different positions between slices. This highlights that for thicker structures without clear pinholes in the oxide, we find weaker, indirect paths of single-electron current (weak-points). We suspect this is due to aluminium rich regions in the oxide enabling a current flow with an indirect path. It is important to note that the edges of the oxide are not abrupt as there is significant oxygen diffusion into the contacts [25, 22]. This stoichiometric gradient results in more current hotspots at the starts and ends of the oxide compared to the middle of a junction. The central slices from Fig. 12 are then the most indicative of the behaviour of the insulating barrier in the junction.

Although these weak-points have the most significant effect in ‘(111)1’ data set, they are nonetheless present in all junctions studied in this work.

Higher transparency channels within the tunnel barrier in Josephson junctions can cause variations in the transmission dependent CPR [21]. The results in Sec. IV indicate higher transparency regions in the normal current density. Following the same methods used in Sec. II.3 for the toy model, we determine the transmission eigenvalues and CPR for each structure.

In this section we also follow recent work which shows that the contributions from different transmission channels can be written as a Fourier series [21]:

where $c_{m}$ are the transmission dependent Fourier coefficients of order $m$.

The ratio of the fundamental sinusoidal signal to the higher order contributions is given by

which allows us to quantify how sinusoidal the CPR is.

Fig. 13 shows the CPR and Fourier components calculated for three structures within the ‘(100)1’ data set.

The magnitudes of the CPRs in Fig. 13 (a) differ by orders of magnitude. For this reason we choose to scale them by the critical current, $I_{\textrm{C}}$, which is equivalent to the maximum for each curve. These values are given in the caption and vary between structure. We find the structures with the least amount of oxygen to have a higher $I_{\textrm{C}}$ compared to the thickest barriers. Scaling the results in Fig. 13 (a) allows us to more clearly study the shape of the CPR. For lower surface coverages we see a deviation from the expected sinusoidal Josephson relationship, $I_{\textrm{S}}(\phi)=\sin\phi$ [21, 41]. When we look at the Fourier transform of each $I_{\textrm{S}}$ curve we find that for junctions with 27.81% and 59.24% coverage there are extra contributions in the Fourier series expansion in Eq. 16. These arise from the higher transmission probability eigenvalues for the thinner barrier structures, and result in a deviation from the sinusoidal relationship. To show this clearly, Fig. 13 (b), (c), and (d) give a breakdown of different sinusoidal contributions to the CPR. Most notably, the thickest junction shown in Fig. 13 (d) has only one significant contribution, the main sinusoidal signal, and is the closest to an ideal Josephson junction. As discussed in Sec. II.3, although we see a systematic change toward sinusoidal behaviour in Fig. 13, this is not necessarily indicative of every structure in the data set. The shape of the CPR is dependent on the number of channels with transmission probabilities close to 1, which in turn depends on the structure of the oxide and the presence of pinholes. Fig. 14 summarises data from Fig. 13 for all the structures in the ‘(100)1’ and ‘(111)2’ data sets, including the transmission eigenvalues for the ‘(100)1’ structure.

Fig. 14 (a) shows the ‘(100)1’ set of data for the eigenvalues of the transmission matrix calculated from the Green’s function. We see that for thinner barriers the eigenvalues are closer to 1, indicating a more transparent junction. As the barriers become thicker the eigenvalues begin to drop to zero. The Fano factor, shown in Fig. 14 (b) is given by

which allows us to differentiate between different types of nanoscale conduction [59]. We find that as we grow thicker barriers the Fano factor approaches 1, indicating barrier tunnelling, and for a pure aluminium we obtain a Fano factor of 0 corresponding to ballistic transport as indicated on the graph with dashed lines. A Fano factor of 1/3 corresponds to a diffusive wire, and 1/2 a 2D layer of randomly distributed scatterers [29, 59]. We note that many of the points in the thinner barriers are situated around the diffusive wire regime, and we do also see a few points near the 2D layer of scatterers regime. However, given the variability in our data we cannot distinguish between these two regimes definitively.

Due to the logarithmic scale in Fig. 14 (a) for the ‘(100)1’ data set, the results in Figs. 14 (b) and (c) for the same data set may show a greater discrepancy (by eye) compared to the eigenvalues in Fig. 14 (a). Fig. 14 (c) shows the ratio of the largest sinusoidal contribution to the CPR with the other contributions, highlighting results in Fig. 13 which indicate that thicker barriers approach a pure sinusoidal CPR and a $\nu\simeq 0$, indicative of an ideal Josephson junction. For the ‘(100)1’ data set, the junction with the strongest sinusoidal response in Fig. 14, with $\nu=0.0071$, has a thickness between 15.25 Å and 23.49 Å based on the FWHM and O extent in Table 1. Fig. 14 indicates that lower oxygen content in the barrier region leads to channels with higher transmission probabilities, resulting in a deviation from the expected sinusoidal response. This is particularly true for thinner barriers where pinholes are more likely to be present. In Sec. IV we found an anomaly in the thickest barrier for the ‘(111)1’ data set, where there were no pinholes present however there were weak-points in the barrier. The term weak-points is synonymous with channels that have higher transmission probabilities.

Fig. 15 directly compares the ‘(111)1’ data set to the ‘(100)1’ for the thickest barriers. Fig. 15 (b) shows the larger proportion of higher transmission eigenvalues for the ‘(111)1’ data set compared to the ‘(100)1’ data set. As we found for junctions with pinholes, these higher transmission channels result in extra contributions to the CPR given in Eq. 16. Fig. 15 (a) shows how this results in a deviation from the expected sinusoidal response, which ‘(100)1’ shows an example of based on its $\nu$ ratio given in Fig. 14 (c) and in the caption of Fig. 15. How pure the sinusoidal response is has a direct influence on the effective Hamiltonian in superconducting qubits. To accurately calibrate qubits it is essential to be able to predict the CPR [21]. Our results indicate that both pinholes and weak-points in the barrier can lead to corrections to the conventional sinusoidal response.

The results of this work show the key role that both pinholes and disorder in the amorphous oxide barrier play in Al$-$AlO${}_{x}-$Al tunnel junctions. Prior work has shown variations in the critical current of Josephson junction devices fabricated in the same way. Using molecular dynamics we have been able to probe the microscopic structure of these devices and analyse the oxidation process step-by-step. We have shown that in thinner barriers pinholes may be present and be a source of single-electron current which could in turn affect the transport properties of Josephson junction devices. Further to this we have also observed a significant amount of variability in junctions with a complete barrier which appears to be caused by the inherently disordered oxide barrier. To understand the minimum amount of oxide required to have a functioning device we studied the metal-insulator transition and found for a uniform electric potential we should see a transition to insulating behaviour between 90-100% coverage. In reality the amorphous structure of the oxide creates a non-uniform electric potential which makes this transition less clear. In fact even when barriers are grown to have thicknesses comparable to experimentally fabricated devices, we might not see a transition to insulating behaviour which could be a contributing factor to device failure. These devices do not have pinholes in the barrier, but have weak-points which can facilitate single-electron current flow.

We find that the variations in transport that arise from the non-uniform electric potential also have consequences on the CPR. More transmissible channels due to pinholes result in more contributions to the CPR, causing it to deviate from the expected pure sinusoidal signal. As barriers are grown thicker we obtain a stronger sinusoidal signal. Although this is the case, we also find that weak-points in the barrier (in addition to fully metallic pinholes) can lead to corrections to the conventional sinusoidal response. This emphasises it is not only essential that a Josephson junction be free of pinholes but also that the oxide barrier be as uniform as possible. Such precise control of junction thickness, density, and stoichiometry may ultimately only be possible with epitaxial junctions [16, 60, 61].

Previously, the specific processes of oxide growth and why thinner films may not function as expected were not well understood. Our work uses computational models to analyse the oxidation process and the formation of structures that may cause faults in junctions with oxide barriers. The analysis in this work aims to inform and support Al$-$AlO${}_{x}-$Al device fabrication.

The authors acknowledge the support of the Australian Research Council through grants CE170100039, and the National Computation Infrastructure (NCI) which provided the computational resources to undertake this work. The authors also acknowledge useful discussions with B. Muralidharan and I. Pop, and assistance of L. Moshovelis with the creation of Fig. 1. We acknowledge the people of the Woiwurrung, Wathawurrung and Boonwurrung language groups as the traditional custodians of the land on which this research was conducted, and the people of the Ngunnawal language group who are the Traditional Custodians of the land of the Australian Capital Territory where NCI is located.