A density-functional theory study of the Al/AlO_x/Al tunnel junction

One of the key components for building a superconducting-circuit based quantum computer is the Josephson junction (JJ), which is essentially a thin insulating layer between two superconductors and ideally facilitates efficient quantum tunneling between the superconductors.¹ The JJ acts as a nonlinear inductor, allowing two states in an otherwise harmonic potential to be selected as the discrete computational basis for a qubit.² The specific parameters of the JJ, particularly the critical current, sensitively affect the resulting qubit energy levels (or, equivalently, operating frequency); thus, predictive control over the JJ fabrication process is desired. However, current fabrication processes limit the ability to a priori predict the resulting JJ parameters within an appropriate tolerance, and the ideal circuit element is limited in reality by imperfect materials, thus requiring empirical process tuning and resulting in inadequate control over critical current dispersion for scalable fabrication.^3,4 Currently, the aluminum oxide tunnel junction fabricated between aluminum superconducting leads has proven to be among the highest performance and best developed JJs for quantum information devices, even as this junction was used for the early observations of the Josephson effect almost 60 years ago.^5,6 However, the efficiency of this junction is far below 100%, and fabrication tolerances are not sufficient to produce perfect uniformity of critical current across a wafer. Better understanding of the atomic structure related to these inhomogeneities and how they may be controlled is needed.

The innovations we have seen over the last 20 years that have led to five orders of magnitude improvement in quantum coherence time of qubits fabricated from aluminum oxide JJs have been mainly driven by clever device engineering and empirical fabrication improvements,^7–9 but still little is understood about the details of junction structure and the correlation with performance. Here, we present detailed atomistic calculations of the structure-processing-property relations for Al/AlO$x$/Al JJs. The main questions we seek to address are (1) What characteristics of the atomic and electronic structure of the junction material (aluminum oxide) determine the JJ parameters and performance? (2) How can atomistic simulations be used to inform improved fabrication processes? and (3) Can JJ performance be predicted for new materials to enable accelerated design of novel fabrications? This work is complementary to experimental investigations of the JJ structure,^10–12 which are challenging and can also provide misleading data for JJ design, as we show later. In addition, the available experimental tools for direct observation of the JJ structure (e.g., electron microscopy and synchrotron radiation studies), as well as empirical experimental fabrication processes, are painstaking, expensive, and low-throughput, inhibiting rapid investigations into new material designs. However, advanced predictive computational modeling, such as based on density functional theory (DFT), can provide a means to accelerate the optimization and design of new, higher performing JJs, especially when combined with experimental studies. When DFT is properly combined with finite element analysis and/or machine learning methods, as we show here, we can reach beyond the atomic and electronic scale to impact the design of an optimal superconducting tunnel junction.

Several prior ab initio computational studies of the aluminum oxide tunnel junction appear in the literature,^13–16 upon which we aim to develop more sophisticated model building of realistic junction structures and enable predictive junction parameter extraction. Some of the cited works used coherent crystalline interface models in relatively small periodic boundaries.^13,14,16 However, because of large lattice mismatch between metal aluminum and its native oxide, experiments commonly find substantial amorphization of the oxide,^11,12 questioning if a small crystalline interface model can represent a real system. In another previous work, a more sophisticated approach was taken: an amorphous solid was separately created, sliced, and placed between crystalline electrodes,¹⁵ but this slice-and-plug way of model making misses the fact that the oxide is supposed to be a surface oxide, and the amorphous oxide is supposed to be grown under the impact of metal/oxide interface. Here, we emphasize the impact of the surface oxidation step on the Al/AlO$x$/Al tunnel junction fabrication process, which has not been addressed in the prior literature. The consequences of the spontaneous amorphization of the surface oxide and the disordered interface structures on junction performance have not been addressed in detail. Without these essential features of the structure, a simulation does not accurately predict the true effective thickness of the JJ, which is a crucial parameter that sensitively affects the critical current.

In this work, we aim to conduct a predictive ab initio study of the materials used in aluminum oxide tunnel junctions, using realistic structural models, and utilize the ab initio data to characterize the JJs at the atomistic level and to derive variations in the JJ critical current. We first conduct structural analyses determining the extent of amorphization of both the oxide and the Al/oxide interfaces, calculate electronic structure to establish a relationship between the degree of oxidation and the bandgap formation, solve the single-particle 3D Schrödinger equation to compute critical current directly from the atomic and electronic structure, and devise a neural network trained with ab initio data to determine local thickness variations of the insulator in an unbiased manner. We find that structure and tunneling properties are related, but in a complicated way. Details at the atomic scale of the structure of the junction material and the uniformity of the superconductor/junction interfaces affect tunneling and critical current sensitively, but can be evaluated predictively with the ab initio methods we present. Crystalline junctions are predicted to be $∼10×$ more transparent than amorphous junctions, and intrinsic structural inhomogeneity of amorphous junctions imparts more than an order of magnitude variation of local critical currents, impacting the effective junction area. We also show that the physical and electronic (tunneling) “thicknesses” of the junctions differ in general.

The atomistic calculations are based on the density functional theory,¹⁷ as implemented in the Vienna ab initio simulation package (VASP).^18,19 The plane wave kinetic energy cutoff for wavefunctions is set to 400 eV, and the electron–ion interactions are represented by the projector augmented wave (PAW) method.^20,21 The exchange-correlation functional by Perdew, Burke, and Ernzerhof (PBE) is used. PBE reproduces the experimental lattice parameter of bulk aluminum within 0.1% and $α$- and $γ$-Al$2$O$3$ within 1.0%.²² Total energy convergence of $1×10−5$ eV and force convergence of $1×10−2$ eV/Å  are achieved for all calculations. A $Γ$-centered $4×4×1$ grid of $k$-points is used for the $p(4×4)$ supercell of Al(111).

A thin surface oxide of aluminum is typically amorphous.¹¹ To create the amorphous oxide junction models, we simulate the oxidation process by beginning with a $p(4×4)$ Al(111) slab model with four Al layers and the bottom two layers constrained to their bulk coordinates. Then, 25 oxygen molecules are placed on the surface with random translation and rotation. We performed NVT ab initio molecular dynamics (aiMD) at 300 K for 2.0 ps using 1.0 fs time steps and a Nosé–Hoover thermostat with the Nosé-mass parameter of 0.15, which resulted in the thermal relaxation on the order of 10 THz, coinciding with typical optical phonon frequencies of aluminum. We found the surface formed a passive saturated oxide layer, and the excess oxygen molecules hovering in the vacuum were removed. We considered the formation of the passive oxide layer complete when the mean displacement of ions reached a plateau (typically less than 2000 time steps). After the simulated surface oxidation, we placed randomly perturbed aluminum atoms on top of the oxide, followed by another NVT aiMD simulation. The first two layers of Al atoms conformed to the amorphous oxide substrate, while the rest recovered the bulk face-centered cubic ordering. After these finite temperature simulation were completed, the model was force relaxed to find a ground state structure that represents the atomic structure of the tunnel junction at low temperatures. (See Fig. 1.) Five separate models were generated in this way to better sample the structural phase space.

We use the radial distribution function, $g(r)$⁠, to analyze the amorphous structure formation,

where $ρ$ is the density of ideal amorphous phase, $4πr2dr$ is the volume of a shell in spherical coordinates, and $dn$ is the number of ions located in the shell. At the limit $dr→0$⁠, $g(r)$ shows the probability of finding an ion at distance $r$⁠. When $g(r)$ converges to 1, it means the material is in a homogeneous phase with no structural order. An amorphous material has short-range order [$g(r)$ shows peaks], but no long-range order [$g(r)$ converges to 1], making $g(r)$ a useful tool for amorphous phase identification.

We also generated a series of hypothetical nearly epitaxial crystalline junction models of different thicknesses for comparison. A $p(5×5)$ supercell of Al(111) and a $p(3×3)$ supercell of Al$2$O$3$(0001) have a small supercell mismatch, less than 1% of both in-plane lattice parameters. Thus, we used these supercells to construct coherent metal/insulator/metal sandwich structures consisting of four layers of $p(5×5)$ Al(111), variable thickness (1–5 stoichiometric layers) of $p(3×3)$ Al$2$O$3$(0001), and another four layers of $p(5×5)$ Al(111). These models were relaxed taking into account the in-plane translation of the Al$2$O$3$ layer as an additional degree of freedom, since the lattice parameter matching alone does not specify the local bonding arrangement of the lowest energy.

The interface regions between Al and Al oxide have complex bonding structure, especially in the amorphous case; in particular, the coordination number of aluminum and oxygen ions vary over extended distances across the interfaces. In order to understand the complex bonding property of the interfaces, we used electron localization function (ELF) analysis followed by electronic density-of-states analyses.

The ELF classifies the nature of interatomic bonding based on topology of the electron density.²³ Metallic bonds show continuous ELF contours, while ionic or covalent bonds show visible localization. The ELF enables us to understand the complex interface between the aluminum electrode and the oxide junction, where strong localization due to under-coordinated ions is found, and also allows us to identify the location of the interface in an unbiased manner (discussed below). Atom-projected density-of-states (PDOS) are calculated using a $Γ$-centered $4×4×1$ grid of $k$-points and the tetrahedron method with Blöchl corrections, and the computed PDOS are averaged over the number of ions.

To determine the junction thickness directly from the ab initio data, we devised a neural network (NN) model to identify the boundaries between metallic Al and the insulating junction material automatically without human bias. The NN model is constructed to identify metallic and insulating regions in space point-by-point using the electronic structure in a surrounding three-dimensional volume, allowing for automatic and consistent boundary extraction, including identification of interface roughness down to the atomic scale. Such an approach is needed to avoid bias from, for example, arbitrary assignment of a flat interface plane location; in addition, the structural symmetry-breaking associated with amorphous material generates possible ambiguity of the location of the interface that the NN resolves without bias.

The NN model is trained on the electronic structure data from DFT, as follows. The cost function of the neural network is defined as

where $m$ denotes the size of the training set, $K$ denotes the size of the solution set [in this case, $K=2$ to classify points in the model as either Al (1) or oxide (0)], and $hθ(x(i))$ denotes the hypothesis function that takes a training input vector $x(i)$⁠. We use the NN parameter $θ$ and the Sigmoid function $g(z)$⁠, where $z=θa$⁠, to propagate activation, $a$⁠, through the NN model. For example, the first activation layer takes the single line of training input, $a(1)=x$⁠, and the activation component of the next layer, $a(2)$ is determined by $a(2)=g(z(2))$⁠, where $z(2)=θ(1)a(1)$⁠. The final activation layer is taken as the $hθ(x(i))$⁠. These NN layers are created for each classification, 1 for metal or 0 for oxide in this work. The electron localization function (ELF) is used as the descriptor, since the metallic part shows dispersive, non-localized ELF, while the oxide part shows strongly localized ELF. The input to the classifier NN was taken as a $11×11×11$ grid-point average of the ELF centered on the grid point of interest. Typical grid point spacing from the DFT calculations was 0.07 Å in each direction. The NN model was trained from the domains of each model clearly identifiable as metal or oxide, far from the interfaces, and tested against a subset of the data.

A regularization function, $R$⁠, is employed to avoid overfitting by penalizing any unusually high NN parameters, via

where $L$ denotes the number of layers of NN, $al$ denotes the number of activation components of the $l$th layer (excluding bias component, $a0$⁠), and $Θl$ denotes the NN matrix of the $l$th layer (we used $θ$ for individual component of the matrix). The strength of this regularization is adjusted by the $λ$ parameter, which has to be optimized.

Since there are two parameters to optimize (⁠$Θ$ and $λ$⁠), we divided the input dataset into three classes: training (⁠$N=12000$⁠), cross-validation (⁠$N=4000$⁠), and final test (⁠$N=4000$⁠), where $N$ denotes the number of assigned data points. From the learning curves, we determined a training set size of 12 000 and regularization parameter $λ=0.5$ gave optimal performance with minimal bias or variance behavior. The learning curves of this neural network are shown in Fig. S6 in the supplementary material. The accuracy of the trained neural network is 96.9% as evaluated from 10 000 unseen random samples taken from all junction models used in this study.

To analyze tunneling (and characterize the critical current) through the junction models, we solved the single-particle Schrödinger equation in a finite element basis in three dimensions using the COMSOL Multiphysics® software v. 5.5 with the Semiconductor module.²⁴ The analysis is based on the formulation for the JJ critical current, computed from the normal electron tunneling, as derived by Ambegaokar and Baratoff²⁵ 

where $Δ(T)$ is the superconducting gap parameter in volts, $T$ is the temperature, $T¯$ is a transparency factor equal to the tunneling probability (⁠$T$⁠) times the density of states at the Fermi level (DOS), $e$ is the electron charge, $h$ is Planck’s constant, and $kB$ is Boltzmann’s constant. Thus, the Cooper pair tunneling rate and JJ critical current may be extracted from a calculation of the electron tunneling probability $T$ through the junction potential barrier.

We obtained the DOS factor by comparing with experimental results for Al/AlO$x$/Al Josephson junctions. Specifically, Dorneles et al.²⁶ showed that the effective area, barrier potential, and thickness can be fit to the Simmons model.^27,28 From analyzing several devices, they fitted the model $RAeff=(10−9Ωcm2)/T$⁠. We may set this equal to the expression²⁵ for $R=2πh/(e2T×DOS)$ to estimate the DOS as $1.62×1014$ cm$−2$⁠. For Al, we take the superconducting gap at $∼10$ mK to be $2Δ=340$ $μ$V.^29,30

To compute $T$ from our atomistic models, we used the stationary study mode in COMSOL with the energy, $E$⁠, of the tunneling electron set equal to the computed Fermi level energy in Al from DFT. The electron potential energy field inside the insulator, $V(r)$⁠, was taken from the self-consistent Kohn–Sham effective potential computed by DFT that includes both Hartree potential as well as the exchange-correlation contribution. The region of non-metal was determined using the ELF and the neural-network model described in Sec. II C. The potential in the metal regions was set so that the electron energy corresponded to the Fermi wavevector of aluminum (11.684 eV, determined from a separate DFT calculation of bulk fcc Al) to describe the Cooper pairs in the superconducting regions. Specifically, the potential for the superconducting metal regions (at the grid point level) was set to an energy corresponding to $Eflat=EF−ℏ2kFmetal2/2m$⁠, where $EF$ is the Fermi energy in the Al, $kFmetal$ is the Fermi wavevector in Al, and $m$ is the mass of an electron. This procedure ensures that the incoming wave is appropriately set simultaneously to the Fermi energy of the junction and the Fermi wavevector of the metal for calculation of the tunneling probability $T$⁠. The effective mass, $meff$⁠, was set equal to that of a free electron. Open boundary conditions in $z$ (the junction normal direction) are used, while periodic boundary conditions are imposed in the $x$ and $y$ directions. Transmission through the junction was computed by impinging an incident plane wave with unit amplitude at the bottom edge of the Al layer and measuring the amplitude of the exiting wave at the top open boundary.

We performed aiMD simulations to approximate the Al/AlO$x$/Al Josephson junction fabrication process, as described in Sec. II and depicted in Fig. 1. This method resulted in a realistic Al/AlO$x$/Al tunneling junction structure with an oxide thickness of about 10 Å, in agreement with experimental measurements derived from cross-sectional scanning transmission electron microscopy (TEM) images.¹¹ 

We characterize the structure of the generated junction models into four regions, denoted (a) through (d) in Fig. 2, based on the analysis of the local structural order using the radial distribution function (RDF), $g(r)$⁠. The RDF approaches unity at large distances for amorphous structures, while distinct peaks are observed even at long distances in crystalline materials. To describe the structure of the obtained JJ models, we begin with the substrate and continue upwards, as this is both how the junction is fabricated experimentally and also how it was constructed with simulation. The Al substrate (d) is crystalline, while the RDF indicates the buried oxide (c) is amorphous. Some number of aluminum layers (b) directly on top of the amorphous oxide retain amorphous character even after simulated annealing. After three or four layers of aluminum deposition, section (a) recovers the Al face centered cubic crystalline order, although the stacking does not necessarily conform to the bottom substrate, (d). Similar results were obtained for the other four generated amorphous junction models, whose RDF profiles are shown in Fig. S1 in the supplementary material. The crystalline junction models are epitaxial by construction and maintain their crystallinity after structure relaxation; the RDF profiles for the crystalline junctions (Fig. S2 in the supplementary material) are characteristic of ordered Al$2$O$3$⁠, even down to 1 monolayer thickness.

The interfaces of Al/AlO$x$/Al show unique bonding characteristics as visualized by the electron localization function (ELF) (Fig. 3). Interfacial ions directly at the transition from metallic Al to insulating AlO$x$ show both localized (and therefore insulating) and delocalized (metallic) characteristics at the same junction depth. Therefore, the transition between metal and insulator is not an abrupt one and varies laterally in space in the direction of the junction. Similar characteristics are observed for the other four amorphous junction models (see Fig. S5 in the supplementary material). This observation questions precisely where the insulating barrier begins in the atomic structure, specifically if the interface should be associated with the presence of oxygen, be associated with proximal aluminum layers, or alternatively begins deeper in the oxide layer. In fact, under-coordinated aluminum ions (e.g., aluminum ions that have less than four neighboring oxygen ions) can retain metallic character according to the ELF.

In order to better understand the insulating barrier, we examine the electronic structure of the junction with the projected density-of-states. This analysis shows that only the fully coordinated aluminum ions in the AlO$x$ layer (coordination number $≥$4) develop finite bandgaps (Fig. 4), while under-coordinated aluminum ions at the interfaces remain metallic. Therefore, coordination number, and not proximity to the amorphous aluminum oxide layer, determines metallicity and the effective thickness of the insulating portion of the junction. We confirmed this finite bandgap formation with the other four amorphous junction models (shown in Fig. S3 in the supplementary material), as well as the crystalline junction models (shown in Fig. S4 in the supplementary material). This result explains the measurement discrepancy between cross-sectional TEM¹¹ vs quantitative XPS,³¹ as we will discuss later in the article.

We deployed the neural network trained to detect the local electronic character to identify the boundary between metal and oxide. Combined with the other analyses, the neural network gives thickness profiles determined from the predicted local interface locations, as shown in Fig. 5, where results for typical amorphous and crystalline tunnel junctions are compared. (Similar analyses for all other models considered in this work are shown in Figs. S7 and S8 in the supplementary material.) These results are based on the unbiased analysis of the differentiation of metallic (superconducting) and insulting regions of the junction models in three dimensions, from which we extract interface profiles and characterize their (electronic) roughness. Table I summarizes the extracted thicknesses, as well as both out-of-plane and in-plane roughness metrics for the bottom and top interfaces of each model. (Distributions of the local junction thicknesses are shown in Fig. S9 in the supplementary material.) The roughness metrics shown in Table I are defined from standard analyses used, e.g., for atomic force microscopy and are defined as follows.³² The roughness average, $Ra$⁠, measures height variations of the interface via

where $z(x)$ is the interface height for a given spatial position $x$ and $z¯$ denotes the average of $z(x)$ over the sample area $A$⁠. The average wavelength, $λa$⁠, measures the characteristic in-plane spacing of peaks and valleys via

where $Δa=1N∑|∂z/∂x|$ denotes mean absolute slope.

The results in Table I clearly show differences between the amorphous and crystalline structures, with the crystalline structures having very regular and uniform interfaces irrespective of junction thickness, while the amorphous structures show significant variations between models and also larger roughness metrics in general. Specifically, all the crystalline models (except the thinnest 1L) show standard deviation of thickness $∼0.4$ Å, while the amorphous models range over $∼0.6$–1.2 Å or 10%–13%. Similarly, the interface roughness metrics for the crystalline models (except 1L) are nearly constant for all thicknesses (⁠$Rabot∼0.2$ Å, $Ratop∼0.25$ Å, $λabot∼2.2$ Å, and $λatop∼3.2$ Å), while greater variation is seen among the amorphous models along with greater roughness (⁠$Ra∼0.35$–1.0 Å) and longer in-plane characteristic wavelengths (⁠$λa∼4$–6 Å). The roughness and in-plane wavelength values for the crystalline models are associated with atomic/lattice dimensions, as expected, while the amorphous interfaces are chemically heterogeneous and disordered. The one-layer crystalline model shows greater thickness fluctuations and interface roughness metrics associated with increased disorder from under-coordinated atoms through much of the thickness, discussed further below. For all models, the top interfaces generally show larger roughness than the bottom interfaces, which is expected since the top interface is formed from an overlayer growth of Al. However, the data show a more consistently higher degree of roughness of the top interface compared to the bottom interface for the crystalline models than for the amorphous, since the amorphous structures already exhibit significant roughness and lateral undulations in the bottom interfaces. Below, we further analyze the effects of the complexity of these interfaces on electron and Cooper pair tunneling behavior.

We performed electron tunneling simulations using the Kohn–Sham effective potential computed with DFT to connect the junction structures with their transport properties. In addition to relating structure to JJ critical current, we observed spatially non-uniform tunneling—a channeling effect—even at very small length scales, associated with quantum mechanical interference associated with different paths through the atomic-scale potential energy variations. We note that the atomically resolved potential energy landscape of the barrier includes portions both above and below the Fermi energy, even though the spatially averaged potential energy difference is positive in the barrier region. By channeling effect, we refer to non-homogeneous transport through the junction with respect to the plane of the insulating layer, that is, certain paths through the junction carry more current than others. The channeling effect is apparent, for example, in the probability density $|Ψ|2$ depicted in Fig. 6, where the incoming electron plane wave is incident from the left and tunnels through the insulating region, approximately demarcated by the transparent blue planes. As an illustration of the atomic-structure dependence of the tunneling, the inset of Fig. 6 shows contours of the probability density inside the junction and enhanced amplitude coinciding with the position of an aluminum ion that is severely under-coordinated inside AlO$x$⁠.

Our tunneling simulations allow us to predict the critical current density for each model using Eq. (3). The results are summarized in Table II. The transmission coefficient $T$ [Eq. (3)] for the five amorphous models was computed between 0.000 64 and 0.016, showing significant variation that is not correlated with mean junction thickness; the corresponding JJ critical current densities range over 170–4100 A$/$cm$2$⁠. These values are higher compared to experimental results for typical JJs used in qubit devices (1–100 A$/$cm$2$⁠, Ref. 33); however, since our models are quite thin and very uniform in thickness over large area (from the small areas of our models), we expect an overestimation of the critical current. We note the high sensitivity of the value of critical current on junction thickness, since the transmission coefficient is exponentially dependent on thickness, and values only a few Å thicker would produce current densities in the experimentally observed range.

The corresponding tunneling and critical current results for the different thickness crystalline junction models are also given in Table II. These are found to have tunneling probabilities and current densities of 0.068–0.27 and $∼20000$–70 000 A$/$cm$2$⁠, respectively. These high values (⁠$∼10$–$50×$ larger than for the amorphous of comparable thickness) are consistent with the regularity of the crystalline junction material, producing greatly enhanced transmission.

In the literature, Josephson junction tunneling is commonly approximated by a rectangular barrier model using the WKB approximation,²⁷ but direct interpretation of such a model is problematic, since the current depends sensitively (exponentially) on the junction thickness and real junctions have spatially non-uniform thickness. Only the thinnest portions meaningfully contribute to the current, and Dorneles et al.²⁶ used experimental measurements to fit effective model parameters of barrier height, junction thickness, and effective area to real junctions, finding discrepancy with physical dimensions and a naïve value of barrier height, e.g., taken from electron affinity differences of Al and AlO$x$⁠. Their analysis found an effective barrier height of 0.81 eV for Al/AlO$x$ JJs.

From our tunneling calculations and by setting our computed junction current equal to that predicted from a rectangular barrier model using the WKB approximation with the width set to the mean thickness derived from our ELF analysis (Table I), we compute equivalent effective barrier heights corresponding to the detailed transmission coefficients determined for each of our models from the direct solution of the Schrödinger equation. The resulting equivalent barrier heights are given in Table II. For the amorphous junction models, effective rectangular WKB barrier heights of 0.20–1.0 eV are obtained, with the variations associated with the local atomic structure in each model. The experimentally derived effective barrier of 0.81 eV falls within this range. The spread in the calculated effective barriers for these different junctions is due to the heterogeneity in their structures, representing different local regions in space. In contrast, for the crystalline models, the effective rectangular barriers are significantly smaller, between 0.04 and 0.11 eV, further demonstrating that the regular structure of the crystalline barriers is much more transparent for tunneling. Interestingly, we note that the values of $T$ for the crystalline cases do not strictly decrease monotonically with increasing barrier thickness; rather, the effective barrier heights calculated fall into two values: $∼0.10$ eV for the three thinnest models and $∼0.04$ eV for the two thicker models. These two different values of effective barrier height for the crystalline cases, depending on the thickness, suggest subtle effects of constructive and destructive interference of the quantum mechanical tunneling that seemingly divides the models into “thick” and ”thin” regimes.

To understand the detailed correspondence between atomic structure and JJ performance, we first analyze the true effective thickness of the insulator part of the junction, according to the electronic and tunneling characteristics. Aluminum ions in the middle of the junction with four or more oxygen neighbors [Al(4) and Al(5) in Fig. 4] develop a finite bandgap, while under-coordinated ions exhibit metallicity. Thus, a thin region of the barrier near the Al interface that contains under-coordinated Al ions provides direct access to the Fermi sea in the metal (superconductor) as well as smoothing of the barrier potential at the interface, enhancing tunneling. Interestingly, this situation may be contrasted with 2D transition metal dichalcogenide (TMDC)-based tunneling junctions:³⁴ while TMDCs are among the thinnest crystalline materials, they are known to suffer high contact resistance. Even when the insulating part of the AlO$x$ junction is only monolayers-thick, the transition from metal (superconductor) to insulator is gradual over atomic length scales. This unique configuration of the electronic structure may contribute to the historically high performance of near-atomically thin Al/AlO$x$/Al tunneling junctions.

Previous work in the literature has discussed how the amorphous nature of the tunnel junction may lead to formation of two-level-system (TLS) fluctuators, which can be detrimental to performance of qubits. The TLS formation is considered an intrinsic property of an amorphous material.³⁵ When TLS is coupled to a tunneling quantum object such as a Cooper pair or resonant with microwave energies in the superconducting circuit, they can cause energy loss (decoherence) from the qubits³⁶ and/or dephase them.³⁷ DuBois et al. used the Streitz–Mintmire interatomic potential to find quantum instability of oxygen ions in aluminum oxide and point out that the quantum noise may be arising from “delocalization of the atomic position of the oxygen.”³⁸ Analysis of the energy landscapes in our structures reveals that oxygen atoms in the amorphous materials, or some atoms at the interface of crystalline materials, most probably contribute to these types of configurational TLS. Our present results contribute an ab initio view on this issue, giving a more detailed and accurate picture into the nature of the amorphous tunnel junction.

An engineering workaround to the TLS problem has been to reduce the area of the junction.^2,7–9 In this way, the number of TLS in the amorphous oxide is minimized simply by reducing the oxide’s volume and possibly also by utilizing a “self annealing effect” associated with nearby free surfaces to relax the defective structures.³⁹ However, this solution does not address the root of the issue (the material and its intrinsic defects), which becomes problematic when scaling up to many devices/qubits or when large junction areas are desired for device design. In particular, the inherent randomness of the material from junction to junction—even when fabricated with an identical process—can severely limit the predictability of individual JJ critical current (qubit operating frequency) and predictable system design for fabrication.

We note that the amorphization is a result of the fabrication process chosen and it is not an innate property of the aluminum oxide material, which could be made crystalline or even epitaxial. We compared our thin crystalline Al$2$O$3$ junction models to the amorphous AlO$x$ models to assess the hypothetical performance limits associated with removing the amorphous nature of the AlO$x$ junctions.

As shown in Fig. 5, the amorphous tunnel junction tends to exhibit significant disordering compared to the crystalline model systems, and evidence of large variations in local properties are evident from comparison of our five amorphous models: in local thickness, interface roughness, and tunneling currents. These intrinsic variations highlight the need to optimize junction uniformity, perhaps ultimately toward crystallinity, which may be challenging with fabrication processes by oxidation. The crystalline tunnel junctions are not only much more uniform with smooth interfaces, but also exhibit $∼10×$ higher tunneling current for comparable thickness.

Regarding the prediction of critical currents and effective barrier heights given in Table II, while the computed values of critical current and effective barrier heights are reasonable, their absolute values should not be taken to directly characterize actual JJ performance metrics, for the reasons given in Sec. III. Rather, these predictions give relative measures of the junction performance for the different atomic structures explored, which might represent real material or portions of those materials in space. The methodology here may be used to evaluate new materials or guide material and process design, particularly if combined with material structure characterization and/or benchmark JJ performance data from experiments. This approach is a predictive, parameter-free ab initio method to compare proposed junction materials, without bias or arbitrary approximation about the tunneling behavior.

We note that the common use of the WKB approximation, for typical junction parameters, is further problematic in being a poor approximation to the exact solution of the Schrödinger equation even for a rectangular barrier model. We demonstrate the problem of using such a model in Table III, which lists the computed transmission coefficient using the WKB approximation for a rectangular barrier with a height of 0.81 eV (taken from Dorneles et al.²⁶) and width taken from each of our models, compared to the “true” values we computed directly for the atomically resolved barrier potentials; the predicted transmission coefficients are orders of magnitude lower than from our solution of the 3D Schrödinger equation. If, for instance, one were to use for the barrier height an ab initio derived electron affinity difference from DFT (we computed 3.052 eV for Al/AlO$x$⁠), the results are even worse. A better model in 1D may be Simmons’ image-charge corrected potential²⁷ solved numerically (rather than the typical approach using the WKB approximation, as done in Simmons’ original work and often found in the literature), but even this (shown in the last column of Table III) does not properly account for the subtleties of the quantum mechanical tunneling through the complicated potential defined by the material at the atomistic level; clearly, there is no proportionality between the true and approximated transmission values. Thus, material evaluation using simulations without experimental benchmarks should account for the true atomistic potential for tunneling through such thin junction materials, as we have done here.

At this point, it is clear that, along with disorder, the thickness is a key parameter for materials design of the tunnel junction. But accurate assessment of thickness is more difficult when the interfaces feature a considerable degree of roughness and disordering. Even with complete knowledge of the atomic scale details of the interface, as we have in our atomistic models, defining a consistent criterion to determine thickness from the structure is not trivial. For example, if we choose to define the interface by discontinuities in the field of Al–O coordination numbers, it is still ambiguous to define the location of the interface where no ions are present. Another issue is that the electronic property of the barrier is not precisely described solely by a simple structural parameter, such as coordination number. Hence, we must use the full information of the self-consistent electronic structure to define the interface as relevant for the tunneling property of the JJ. To utilize the electronic information to determine the location of interfaces, we use the electron localization function (ELF). ELF can be an ideal feature for machine learning in this case, because both the absolute value and the derivative of ELF contain useful information to determine the chemical property of local domains based on the metallic vs non-metallic character: the metal domain features delocalized and continuous ELF, while the oxide features strongly localized ELF rapidly varying between 0 and 1. While thickness in the context of JJs is usually considered a one-dimensional quantity, we point out that the detailed nature and statistics of the spatial variations and roughness of the thickness is critical to JJ performance and that mean thickness provides only an incomplete characterization of the junction material, as apparent in Fig. 5.

In addition to the AlO$x$ regions of the junction, the interface Al metal has interesting structural features to be pointed out as well. Figures 1–3 illustrate that the top Al grown above the AlO$x$ layer includes a thin interface layer that is highly disordered/amorphous, compared to the bottom substrate/AlO$x$ layer which maintains the crystalline order up to the interface. The visible structural disorder in the top Al seen in Fig. 1 is quantified by the RDF analysis in Fig. 2; the amorphous nature of the oxide (c) and top Al interface region (b) is apparent by the flat RDF tending to 1.0 at large $r$⁠. The first few layers of metal on top of the AlO$x$ naturally become amorphous because its substrate (AlO$x$ part) is amorphous. After the first 2–3 layers of metal deposition, Al already mostly recovers the face centered cubic (fcc) ordering, as evidenced by comparison of the peaked RDF for the top Al layers (a) with the perfectly crystalline bottom substrate (d). In other words, the AlO$x$ tunnel junction experiences the fundamental asymmetric fabrication problem as pointed out by Copel et al.⁴⁰ 

There have been a number of previous computational modeling works on Al/AlO$x$/Al system,^13,13–16 each of them contributing to understanding the microscopic details of the aluminum oxide tunnel junction. A main aspect missing in the literature has been that the atomistic models used in the earlier works were limited in accurately representing the interface structure of the aluminum oxide tunnel junction, which may be garnered only by explicitly simulating the surface oxidation. Some earlier works relied on building a sandwich model consisting of Al fcc metal and Al$2$O$3$ attached coherently in a small unit cell, which results in unrealistic bonding and bond lengths near the interface. A significantly improved model building procedure encompassed first creating an amorphous aluminum oxide with a variable density parameter, then selecting the structure whose density agreed with experimental data; the amorphous aluminum oxide model was then sliced to arbitrary thickness and placed between crystalline metal electrodes, followed by force-relaxation.¹⁵ In our work, we found that a realistic amorphous oxide structure is spontaneously created when we explicitly simulate the surface oxidation of the aluminum electrode. In the earlier works, the resulting tunnel junction model could not predict the thickness of the junction formed by oxidation, because the thickness was chosen. Because we simulate the oxidation process in our model creation, the resulting junction thickness is a realistic representation of the rapid oxidation process used to create thin junctions, and the oxide thickness in our AlO$x$ models—as well as structure—closely matches experimentally fabricated JJs without input parameters. Our simulations also demonstrate probably a lower limit of amorphous junction thickness for conditions where oxygen can saturate the surface (typical fabrication conditions), which is consistent with the thicknesses of our amorphous junction models being on the thin end of fabricated junctions.

Of course, the simulation time scale (and system size) for the ab initio approach we used is limited and cannot match experimental scales, where oxidation may occur over minutes or longer. However, our approach is able to capture the relevant thermodynamics of the process and effectively accelerates the kinetics via our iterative model creation process so that the resulting structures are realistic and compare well to experimental structural characteristics.^11,12 For example, the RDFs computed from our ab initio-generated amorphous models generally compare well to those generated from transmission electron microscopy by reverse Monte Carlo.¹¹ In fact, structural models as we generated may be fit even more consistently to the experimental images. Surface oxidation of aluminum is known to be very rapid, and the resulting oxide is passivating once formed.^31,41 In our MD simulations, when oxygen is introduced to the Al surface, the formation of a passive oxide film is deemed complete when the ionic displacements reach a plateau as shown in Fig. 1, which typically occurs by 10 ps of aiMD under the conditions we used. The subsequent slow kinetics of the film formation are not explicitly captured by the aiMD simulations, but the final energy relaxation step mimics the longer timescale annealing of the films.

The literature includes experimental measurements of the thickness of Al/AlO$x$/Al tunnel junctions that appear inconsistent with each other. When measured by quantitative XPS, the aluminum oxide layer seems to be much thinner³¹ compared to measurement from cross-sectional transmission electron microscopy (TEM).^11,42 Our result shows that this is not a measurement error. While XPS is sensitive to chemical changes, TEM images capture mostly the structure distortions. As shown by our results (Fig. 4), full chemical modification of aluminum is only achieved near the middle of the oxide layer. For example, the apparent physical thickness of the junction in our Amorphous E model is approximately 1 nm when including the full width of the oxygen containing region, corresponding to the interpretation of TEM measurements, while the thickness of the fully developed near-stoichiometric oxide is only about half this, corresponding to the interpretation of XPS measurements. Overlapping layers and interface roughness make cross section TEM analysis difficult to determine precise interface boundaries. Possibly XPS measures the fully chemically transformed aluminum layer (about 0.5 nm into the interior of the oxidized layer), while TEM determines the thickness that includes the interface layers. How we define the thickness of the junction has critical consequences on how we model and design quantum tunneling devices, because the thickness is a key parameter of the junction as quantum tunneling decays exponentially with distance. If one determines the thickness as a distance between Al atoms that have a bulk lattice structure, then the thickness is considerably larger than the other case when we define the thickness from the electronic property of the junction. As illustrated by Figs. 3 and 5 and Table I, the location of the interface defined by its electronic properties, which is most relevant for the JJ performance and critical current, does not necessarily match obvious structural features that can be readily identified in TEM images, for example.

To assess the variance in Josephson junction physical parameters and performance, Dorneles et al.²⁶ fit junction thickness, effective area, and potential barrier height to the resistance of 17 devices. The fitted effective area of these junctions varied by two orders of magnitude for devices that were designed to all have the same area, and the effective area was found to be five orders of magnitude smaller than that defined by lithography. This behavior was attributed to “hot spots” in the junction that contribute most of the current and form due to uncontrolled thickness variations in the oxide layer. In our simulations, we find more than an order of magnitude difference in transmission (or critical current) between similarly created (amorphous) oxide barriers. Since our models are small in area and represent statistical spatial samples of a larger area amorphous material, we interpret the models exhibiting higher current as corresponding to higher current “hot spots” of the Dorneles study. Of course, thickness variations in the oxide can create large differences in current due to the exponential dependence of current on thickness; however, we find that approximately an order of magnitude variation in current can also arise from local atomic structure variations for junctions with a similar thickness. Therefore, the “hot spots” observed by Dorneles et al. may be caused by the inherent non-uniformity of the local atomic structure of the amorphous junction material, in addition to gross thickness variations.

Our prediction of the metallic states of under-coordinated aluminum ions is not entirely new. Rippard et al. noted that oxidized ultrathin aluminum exhibits conducting channels until the Al atoms become fully oxidized.⁴³ They used ballistic electron emission to characterize the single electron channel effect of the ultrathin aluminum oxide tunneling junction. Our Schrödinger solution of the tunneling electron (Fig. 6) corroborates this observation. The tunneling probability, $|Ψ|2$⁠, suggests sparse conducting channels with orders of magnitude higher tunneling probability compared to neighhboring areas. The locations of these high-conductance channels appear to coincide with the position of severely under-coordinated aluminum ions at the oxide interface with only one oxygen neighbor. The details of this effect require further investigation to establish a clear relationship between the structure and enhanced local Cooper pair tunneling.

Both experimental characterization and ab initio simulation of the atomic structure of the tunnel junction suffer a common source of error, which is that determining the precise location of the interface by ad hoc means can be arbitrary. We have resolved this issue in this work by employing artificial intelligence, through the training of a neural network to determine the location of the interfaces in an unbiased way from the local electronic structure determined by DFT. We used this strategy throughout the present work, and we were able to obtain more consistent and robust characterization of the junction material and interfaces than possible by manual assignment. In addition, we found this approach to be essential for accurate description of electron (Cooper pair) tunneling through the junctions; comparison to calculations that assumed flat interfaces by arbitrarily assigning an interface position averaged over the junction area produced significantly different values of junction transparency, which were not reliable metrics of the relative critical currents.

The electronic structure of Al/AlO$x$/Al computed by DFT shows that the insulating part of thin junctions can be only a few monolayers thick, while significant portions of the interface regions contain under-coordinated aluminum atoms that are conducting. This structural transition at the interface, and the associated electronic properties of the under-coordinated Al atoms, explains the discrepancy between reported XPS and TEM measurements of the thickness of Al/AlO$x$/Al tunnel junctions, since XPS captures the chemical modification in the few monolayer-thick insulating part, while obscured contrast in cross-sectional TEM images would tend to include the under-coordinated interface atoms as part of the junction oxide. Oxides formed from surface oxidation of aluminum show asymmetry of the bottom and top interface layers, demonstrated by our amorphous oxide models generated using simulated oxidation and metal overlayer growth. The amorphous nature of the grown oxide causes the first few layers of overgrown Al to be amorphous as well, although crystalline order is restored after 2–3 monolayers of Al.

We used explicit 3D numerical solutions of the Schrödinger equation to study electron (Cooper pair) tunneling across realistic atomic structures of junctions, comparing both amorphous and crystalline materials. The results revealed a channeling effect for tunneling associated with local structure of the oxide, as well as large variations of junction transparency across different areas of at least an order of magnitude for the amorphous material. Even different amorphous models with similar mean thicknesses exhibited differences in tunneling current as large as $10×$⁠, showing that the specific degree of disorder in the junction atomic structure has large effects on performance. Similarly, our analysis also shows a significantly higher transparency (up to $∼50×$⁠) for an epitaxial crystalline junction compared to amorphous, for comparable thickness. The results further demonstrate the extreme sensitivity of JJ tunneling performance on details of the atomic structure. Junction thickness is a critical parameter, since tunneling is exponentially dependent on thickness, but thickness is a difficult-to-define parameter. We used a trained neural network model based on the local electronic structure to objectively define the position of the insulating interface spatially across our models. Interface roughness parameters, both out-of-plane and in-plane, further influence the junction performance; for junction transparency, the important metrics relate to electronic roughness, which is related to but may not be identical to physical roughness. The distribution of under-coordinated atoms in the material resulting from the fabrication process also influence the properties.

The whole of our analyses in this work illustrates an approach using ab initio calculations to assist the materials and process design of tunnel junctions. The atomic structure of the junction material, including thickness variations, under-coordinated ions, and interface roughness, was taken into account in a predictive manner and connected to parameter-free models of junction performance (transparency, critical current), ultimately to understand how structure affects performance. Our ab initio models gave reasonably predictive results for structure and gave qualitatively useful critical current behavior. The overestimation of critical current likely originated from models that are slightly too thin and from the small junction areas of our models that can only capture structural variations on very small scales, although statistical averaging over the results for multiple models could account for this affect. However, comparison of relative performance of different junction materials/structures could be reliably obtained with our predictive, parameter-free methods, which can be a useful tool for evaluation of new materials proposed in the future. We also demonstrated that care must be taken in the selection of tunneling models for meaningful results that can be used for design; rectangular barrier models with the WKB approximation often are insufficient for even crude approximation of the tunneling physics.

From the analyses presented in this work, design principles may be derived for the next generation of quantum tunnel junctions, using ab initio methods to compare candidate materials, structures, and processes, beyond overly simplified metrics of thickness and apparent process uniformities. The atomic-scale details of non-uniformity in the junction material and its interfaces are critically important, and we present direct calculations to predictively evaluate the magnitude of the impacts on JJ performance. Detailed experimental structural characterization may also form an input to calculations as shown here, enabling a better understanding of the tunneling physics of newly developed JJs.

See the supplementary material for detailed data on all of the models investigated in this work not shown in the main text. Particularly, radial distribution functions and structural analyses of each model are provided, along with projected density of states and electron localization function for electronic analysis. In addition, detailed data on thickness profile analyses are provided, including profiles and thickness distributions for each model, along with learning curve data for the neural network used to identify interface locations.

This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract No. DE-AC52-07NA27344, funded by the Laboratory Directed Research and Development Program at LLNL under Project Tracking Code No. 18-ERD-039. C.-E.K. gratefully acknowledges advice from David F. Bahr and his support for additional computing resources from Purdue University Research Computing.

The data that support the findings of this study are available within the article and its supplementary material and are available from the corresponding author upon reasonable request.