Effect of atomic structure on the electrical response of aluminum oxide tunnel junctions

Many nanoelectronic devices rely on thin dielectric barriers through which electrons tunnel. For instance, aluminium oxide barriers are used as Josephson junctions in superconducting electronics. The reproducibility and drift of circuit parameters in these junctions are affected by the uniformity, morphology, and composition of the oxide barriers. To improve these circuits the effect of the atomic structure on the electrical response of aluminium oxide barriers must be understood. We create three-dimensional atomistic models of aluminium oxide tunnel junctions and simulate their electronic transport properties with the nonequilibrium Green’s function formalism. With this approach we are better able to understand the role that ﬂuctuations in the density and stoichiometry of the oxide play in the electrical response of the junctions. For instance, increasing the oxide density produces an exponential increase in the junction resistance. In addition we observe the formation of metallic channels in highly oxygen-deﬁcient junctions. We ﬁnd that local variations in density or stoichiometry can lead to localized conduction channels, even for a junction of uniform thickness. The atomistic approach we have taken provides a better understanding of these transport processes and guides the design of junctions for nanoelectronics applications.


I. INTRODUCTION
Superconducting qubits are one of the most promising architectures for quantum computers and are currently the favored technology for many quantum computing groups around the world [1][2][3][4][5]. These qubits rely on the nonlinear response of Josephson junctions, which are often fabricated as Al−AlO x −Al trilayer junctions [6][7][8]. As we move to largescale quantum computer engineering, it becomes critical to understand what limits junction performance and variability.
The conductance of Al−AlO x −Al junctions is commonly understood in a simplified one-dimensional picture. In the simplest case, the tunnel junction is considered to be a rectangular barrier where the transmission probability of an incident electron can be calculated using the Wentzel-Kramers-Brillouin (WKB) equations [9]. More detailed analytic models of the tunneling barrier include corrections for temperature, applied voltage, image forces, and asymmetries [10].
Two of these models-the Simmons model [11] and the Brinkman, Dynes, and Rowell model [12]-are often used to estimate parameters such as the barrier height and the oxide thickness by fitting to experimental measurements [13][14][15][16][17][18]. Barrier heights calculated by fitting to the Simmons model [15][16][17] range from 0.8 to 3.0 eV while a "typical" height of 2 eV is often quoted [19][20][21][22]. Estimates of the barrier height and oxide thickness given by such models are effective values which include contributions from oxide properties such as the density and stoichiometry implicitly. While useful, one-dimensional descriptions of the barrier system are unable to fully represent the amorphous oxide layer.
To include the full three-dimensional structure of the junction we turn to a numerical approach. There is a growing body of literature in which the nonequilibrium Green's function (NEGF) formalism is used to calculate the electronic properties of nanoscale devices. This is a numerical method which allows us to calculate properties such as the transmission probability, current, and charge density. A range of systems have been studied with this approach including graphene, silicon, and phosphorus-in-silicon nanowires and carbon nanotubes [23][24][25][26][27][28]. These systems all consist of regular repeating units which allow for the calculation to be performed in reciprocal space, potentially yielding improvements in computational efficiency. However, the Al−AlO x −Al junctions which are the subject of this study are inherently disordered; this removes any symmetries we might exploit to reduce the complexity of the problem.
The computational challenges that arise when dealing with disordered systems may explain the small number of first-principles calculations in the literature with a focus on Al−AlO x −Al junctions. One study by Zemanová Diešková et al. [29] presents ab initio transport calculations for small atomistic junction models. The conductance was calculated using a transfer matrix method and compared to the conductance of rectangular and trapezoidal barriers as well as an sp-like tight-binding model. A ground-state ab initio simulation is used to determine the parameters of the tight binding calculation. Relatively poor agreement with experimentally FIG. 1. An atomistic model of a Josephson junction (d = 14 Å, ρ = 0.7, γ = 1.1) created using a simulated annealing method. Aluminium and oxygen atoms are shown as gray and orange spheres, respectively. The yellow and blue regions correspond to the central and interfacial parts of the oxide barrier, respectively, which are referenced in Sec. IV B and Fig. 4 reported conductances is observed. Inaccurate estimation of the barrier thickness with the Simmons model is raised as a possible cause for this discrepancy.
In this paper we use molecular dynamics techniques to create three-dimensional models of Al−AlO x −Al junctions that include the detail of the atomic structure. The shape of the potential barrier-used as an input to our electronic transport model-is calculated in three dimensions from the atomic positions and charges. The electronic properties of the junction models are calculated with the NEGF formalism [30]. Due to the native disorder in the oxide noted above we calculate solutions to the NEGF equations for a three-dimensional real space representation of the system. By starting with a model of the atomic structure of the Al−AlO x −Al junction and retaining a full three-dimensional description of the structure through each part of our calculations, we probe the effect of structural changes on commonly measured quantities such as the junction resistance.
The structure of the present work is as follows. In Sec. II we describe our approach for creating atomistic models of Al−AlO x −Al junctions. The nonequilibrium Green's function formalism used to calculate the electronic properties of the junction models is laid out in Sec. III. The results presented in Sec. IV demonstrate the way in which changes in the material properties of the oxide layer such as thickness, stoichiometry, and density affect the junction resistance. Variation in the local structure at the Al/AlO x interfaces is shown to affect the uniformity of current flow through the junction.

II. ATOMISTIC JUNCTION MODEL
To study the effect of atomic scale structure on the electronic properties of the junction we create atomistic models, starting with a large supercell of crystalline Al 2 O 3 (corundum). We adopt a convention where the thickness of the barrier d is measured along the z axis while x and y are the lateral directions. The variable ρ is used to describe the density of the junction oxide as a multiple of the density of crystalline Al 2 O 3 (3.97 g cm −3 ) [31]. The variable γ represents the stoichiometric ratio of oxygen to aluminium in the center of the oxide (the yellow region in Fig. 1). A crystalline Al 2 O 3 structure would therefore be described by values of ρ = 1.0 and γ = 1.5. By modifying the Al 2 O 3 crystal we produce oxide structures with a range of densities and stoichiometries. Bulk amorphous AlO x is experimentally reported to have lower density and stoichiometry than the crystalline structure [32]. To create an oxide barrier of a given thickness d but a reduced density, a volume is cut from the corundum supercell of size x × y × ρd after which the structure is expanded in the z direction by a factor of ρ −1 . The desired stoichiometry is then obtained by randomly removing oxygen atoms from the structure. Following this, a geometry optimization is performed to find the lowest energy configuration of the atoms during which the atoms are free to move, but the size of the simulation box is fixed to ensure that the density remains constant. We use the General Utility Lattice Program for both this optimization and the subsequent molecular dynamics calculations [33]. Interations between the aluminium and oxygen atoms are described with an empirical potential parameterized by Streitz and Mintmire [34].
To introduce disorder in the structure we run a molecular dynamics calculation at 3300 K (which is 1000 K above the melting point of corundum) for 4 ps with a time step of 1 fs. Following this the simulation temperature is linearly reduced to 300 K over 6 ps to quench the oxide in a specific disorder configuration. Crystalline aluminium regions are then placed adjacent to the oxide (in the positive and negative z directions) and a second geometry optimization is performed to reconstruct the interfacial regions between the oxide and the aluminium contacts [35]. During this optimization the box can expand or contract along the z axis, and atoms in the aluminium contacts and up to 4 Å into the oxide on each side are free to move. By fixing the atoms in the central region we are able to retain the desired density and stoichiometry even in cases where the final structure may not be energetically optimal. An example of an atomistic junction model produced in this way is shown in Fig. 1.
Experimental studies of aluminium oxide structure [32,36] show that the amorphous phase has a density approximately 0.8 times that of the crystal phase and a stoichiometry of γ = 1.10. Though epitaxial stoichiometric Al 2 O 3 junctions have previously been realized [38], we focus on the more common amorphous barrier systems. The oxide layer in a single junction [37] varies from 10 to 20 Å. On the basis of these values we create three sets of junction models summarized in Table I. In each data set we vary one parameter while keeping the other two fixed at realistic values based on experimental data.
The lateral dimensions of each junction model are x = y = 24 Å, whereas those of junctions in real circuits usually exceed 100 nm [7,18]. For this reason periodic boundary conditions are applied during the development of the junction model.

III. ELECTRONIC TRANSPORT MODEL
In order to calculate the electronic transport properties of the atomistic Al−AlO x −Al junction models we implement a nonequilibrium Green's function model in one, two, and three dimensions. This provides a single-particle description of a free electron moving through a disordered potential obtained from the atomistic model. A finite difference approximation of the kinetic operator is used to describe the channel through which transmission occurs and the source and drain contacts which are connected to it. In the three-dimensional transport calculation the boundary conditions are periodic in x and y, with open boundary conditions in z.
To second order, the finite difference representation of the kinetic energy for a one-, two-, or three-dimensional system can be written as where k ∈ {x, y, z}. The magnitude of the on-site energy ε changes with dimensionality: . ( The magnitudes of the hopping energies t k =h 2 /2m * a 2 k are determined by the spacing between points a k 1/3 Å and the effective mass m * . We choose m * to be the free electron mass m e as the model is designed to describe electrons tunneling between two contacts composed of bulk aluminium in which m * m e [39]. The electrostatic potential V (x, y, z) in the junction structure is calculated on a Cartesian grid and added to the kinetic energy T to form the complete channel Hamiltonian H C = T + V . Details of the numerical approximations made when computing the electrostatic potential are given in Appendix A 1. To obtain the transmission function T (E ) we calculate the retarded Green's function, where I is the identity matrix, iη is a positive imaginary infinitesimal number, and S (E ) and D (E ) are the selfenergies for each contact where the subscripts S and D denote the source and drain, respectively. The matrix inversion in Eq. (4) is computationally expensive, and we take advantage of a recursive algorithm to speed up our calculations [40]. The trace over the product of the broadening matrices S,D = i( S,D − † S,D ) and the retarded Green's function yields the probability of transmission through the channel as a function of the energy of the incoming electron: In the Landauer-Büttiker formalism [30] we can use the value of T (E ) to evaluate the current in the channel as a function of applied bias: where e is the charge of an electron, h is Planck's constant, and f i (E ) is the Fermi-Dirac distribution for contact i: where μ 0 is the Fermi level, k B is Boltzmann's constant, and T is the temperature. The junction is symmetrically biased such that V S = −V D . Equation (6) could also be used to determine the junction resistance from the gradient of the linear I-V response at low bias. The computational cost can be reduced by working in the limit V → 0 where we can use the zero-bias conductance formula A further optimization is obtained by taking the zerotemperature limit where Eq. (9) becomes a δ function centered at μ 0 . This gives us an expression for the resistance which requires only the evaluation of the transmission function at a single energy: where G 0 = 2e 2 /h is the conductance quantum. When the transmission function is smoothly varying (see Appendix A 2), this becomes a sufficiently good first approximation at higher temperatures. We report the resistance area R N A given by the product of the normal resistance [calculated with Eq. (10)] and the area of the simulation cell transverse to the conduction direction. We choose to calculate the resistance area because it is commonly measured in experiment and can be calculated assuming normal state conduction. With the Ambegaokar-Baratoff relation we are then able to link the resistance in the normal state with the critical current of the device when it is superconducting [41].
Calculations of the current and resistance [with Eqs. (6) and (10), respectively] depend on the value of the Fermi level μ 0 . To estimate μ 0 we fit our simulation to an experimental value of the resistance area. For a reference junction (with typical thickness, density, and stoichiometry) we calculate the resistance area for a range of energies. The Fermi level is then found by matching the calculated resistance area with a representative experimental [18] resistance area of 600 μm 2 . For our data set this gives a value of μ 0 = 1.35 eV. In Sec. IV we are limited to a discussion of qualitative trends only as variation in μ 0 , which would occur if a different junction or 013110-3 experimental value was chosen as a reference point, leads to an offset in the calculated resistances for the junction models.
We can also calculate electronic properties which vary spatially within the junction structure. The charge density in three dimensions is given by with the electron Green's function and The current flowing between two points can be determined from the element-wise product of H and G n : The net current in a particular direction is then calculated from the difference between pairs of points and normalized by the area of the discretization in the other two directions. For example, where r = (x, y, z) and r = (x, y, z + a z ). The expressions for J x and J y are constructed similarly. It is worth noting that the equations presented here are entirely general to any one-, two-, or three-dimensional transport system that is well described by a nearest-neighbor finite difference model.

A. Current-voltage response
The response of a tunnel junction to an applied bias is expected to be linear when the bias is close to zero and to become nonlinear as the bias is increased. When a sufficiently large bias is applied the Fowler-Nordheim tunneling theory can be used to describe the response [42]. The current-voltage relationship for an atomistic junction model, calculated in three dimensions with Eq. (6), is shown in Fig. 2. We observe a linear response at low bias and find that the behavior is well described by the Fowler-Nordheim tunneling model above an applied voltage of approximately 2 V. It should be noted that this is well above the typical experimental breakdown voltage for junctions and is used here simply to benchmark the technique.
The current in the Fowler-Nordheim model is given by [9,43,44] where α is a scaling factor related to the proportion of the barrier which participates in tunneling via field emission, A is the cross-sectional area of the device, φ is the work function, and β is the inverse of the barrier thickness d. The quantities a and b are the Fowler-Nordheim constants, which are  [32,36,37]. At low bias a linear response is observed, while agreement with the Fowler-Nordheim model is seen at high bias.
given by where m is the effective mass of the electron in the oxide. An estimate of the work function φ can be found by fitting the calculated current-voltage data with Eq. (17). The effective mass used during the fitting process affects the calculated value of the work function. An effective mass of m = 0.4 m e estimated from band structure calculations [45,46] for crystalline Al 2 O 3 yields a value for the work function of φ = 2.4 eV. Alternatively, direct measurement of aluminium oxide barriers [47] gives an estimate of m = 0.75 m e leading to a value of φ = 2.0 eV. Both values are close to the commonly quoted barrier height [19][20][21][22] of 2 eV.

B. Effect of oxide morphology on resistance area
The resistance area is calculated with Eq. (10) for each junction in the three data sets described in Table I. The resistance area as a function of oxide thickness is shown in Fig. 3(a). A linear fit to the log of the resistance area data is calculated using MATLAB. This data set consists of 18 junctions with approximate thicknesses of 10-30 Å and densities in the narrow range ρ = 0.77−0.87. The exponential increase in the resistance area with barrier thickness is in agreement with experimental observations [15] and an exponential reduction in the tunneling probability. Figure 3(b) shows the relationship between the density of the barrier oxide and the resistance area. Here we observe that the resistance of the junction is also exponentially related to the oxide density. We note that each junction in this second data set has a similar thickness (d = 16 ± 1 Å). The height of the potential barrier increases with increasing density [see inset in Fig. 3(b)]. To the authors' knowledge, no systematic studies exist investigating the relationship between the junction resistance and the oxide density. Sullivan et al. [32] report that oxides manufactured with an O 2 plasma deposition process are of higher density (ρ = 0.8) when Al is evaporated simultaneously and lower density (ρ = 0.6−0.7) when the substrate is exposed only to the plasma. From Fig. 3(b) we can estimate that this variation in the density would correspond to change in the resistance area of 1-2 orders of magnitude.
Resistance-area data are presented in Fig. 4(a) for a range of oxide stoichiometries. Between γ = 0.9 and γ 1.2 the resistance area is approximately constant. This range is comparable to reported experimental values for oxide stoichiometry of γ = 0.8−1.2 (depending on fabrication conditions) [36]. A significant drop in the resistance is seen for values of γ outside this region. At low stoichiometries (γ < 0.9) this is due to oxygen deficiency in the junctions creating metallic channels which dominate conduction and lead to a decreased resistance. We can therefore identify the metal-insulator transition at approximately γ = 0.8−0.9.  4. (a) The calculated resistance area varies over several orders of magnitude as the structure moves from an oxygen deficient to an oxygen-rich configuration. Charge and current densities are presented in Fig. 6 for the structures corresponding to the three data points marked with red arrows. The shaded data points correspond to the current density histograms presented in Fig. 7. The upper and lower bounds of experimentally reported stoichiometries [36]  At high stoichiometries (γ > 1.2) the resistance drops again, which is the result of fixing the density at ρ = 0.8 to systematically compare between different stoichiometries. The reduction in resistance at high stoichiometries for this set of junctions can be attributed to changes in the distribution of oxygen within the barrier. Two additional data points are included at ρ = 0.9 and ρ = 1.0. As expected from the exponential dependence shown in Fig. 3(b), increasing the density increases the resistance. Figure 4(b) helps to explain the change in the oxygen distribution at higher stoichiometries (γ > 1.2). We define the stoichiometry in the interfaces as γ interface and plot how this changes as a function of the stoichiometry in the center of the barrier γ . The central and interfacial regions are defined in Fig. 1. We observe that γ interface begins to decrease at higher values of γ . This implies that the interfacial regions between the contacts and the oxide barrier tend to reach a limiting stoichiometry of γ 0.8 irrespective of the higher oxygen concentration in the center of the barrier. Constraints used during the preparation of the structure such as requiring a density of ρ = 0.8 and limiting the motion of atoms during the optimization are forcing the system away from thermodynamic equilibrium and may drive this variation in stoichiometry across the oxide. The oxygen-deficient interfaces are more conductive and cause a decrease in the effective thickness of the tunneling barrier leading to the observed decrease in the resistance-area product. We note that while high conductance is observed at both low and high stoichiometries (when ρ = 0.8), the transport in the high stoichiometry region is still in the tunneling regime where T (E ) < 1.

C. Charge and current density
To better understand how conductance changes as a function of stoichiometry, we calculate the charge density and current density in three dimensions for junction models with stoichiometries of γ = 0.3, 0.9, and 1.5 (corresponding to the three data points in Fig. 4 marked with red arrows). These properties were computed with Eqs. (12) and (16) at an applied bias of 50 mV. In Fig. 5 we plot n(y, z) and J (y, z) for a junction in the tunneling regime (γ = 0.9) for three planes at different positions along the x axis. Lighter regions with lower charge density are associated with the presence of oxygen. The disorder in the atomic structure of the oxide can be observed in the contours of the charge density in the barrier region. The current density varies as a function of x and y with regions of higher current around the center of Fig. 5(c) and on the bottom of Fig. 5(d).
Variation in the physical thickness of the oxide layer has been observed directly in microscopy studies where it is estimated that less than 10% of the total barrier area dominates the tunneling of electrons [37]. This gives rise to the idea that there are localized regions of higher conductance or "hot spots" in the junctions. For example, Aref et al. [18] report an effective conduction channel size of 20 nm 2 . Our results demonstrate that the effective width of the tunneling barrier can be affected by small local differences in the density of aluminium and oxygen atoms at the oxide/metal interfaces. This is evident even in our junction models with minimal variation in physical thickness across the structure. A comparison of the calculated current for the various stoichiometries is shown in Fig. 6. The charge density contours in Fig. 6(a) (γ = 0.3) show a significantly weaker suppression of the current than is evident in the insulating γ = 0.9 junction. The low stoichiometry structure (γ = 0.3) contains small regions of aluminium oxide that do not span the entire lateral width of the junction model, leaving metallic channels through which the majority of the current flows.  It is important to note here that the arrows depicting the current density are 10 3 times smaller than those in Fig. 6(a). Figure 6(b) (γ = 0.9) corresponds to a junction in the fully insulating regime, while a path of higher current density can be seen at the top of Fig. 6(c) (γ = 1.5). Oxygen deficiency in the Al/AlO x interfaces creates areas where the insulating barrier is thinner and electrons can more easily tunnel through the oxide. The current densities presented in Fig. 6 allow us to understand the drop in the calculated resistance-area values (at both low and high stoichiometries) in Fig. 4.
The formation of localized conduction channels in oxides with high stoichiometries and reduced densities can be seen directly by examining the current density in the xy plane perpendicular to conduction. In Fig. 7 we illustrate this with three examples: (I) high stoichiometry and reduced density FIG. 7. Histograms of the current density in xy slices through the junction structure averaged over the central region of the oxide (i.e., the yellow region defined in Fig. 1). The three series correspond to the data points marked on Fig. 4 in the same colors. The inset figures are representative of typical xy current densities in these different systems with the given densities and stoichiometries.
(γ = 1.5, ρ = 0.8), (II) high stoichiometry and high density (γ = 1.5, ρ = 1.0), and (III) typical density and stoichiometry for an insulating barrier as observed experimentally (γ = 0.9, ρ = 0.8). The histograms in Fig. 7 are of the current density in the xy plane averaged over the central region of the oxide (i.e., the yellow region in Fig. 1). From these distributions we see that when a localized conduction channel forms there is also a corresponding reduction of current in the surrounding region [see Fig. 7 (I)]. In contrast, a more typical barrier conducts current uniformly across its entire lateral extent [ Fig. 7 (III)].
The localized conduction channel visible in Fig. 7 (I) is approximately 2-3 Å in diameter, which is substantially smaller than the usual "hot spots" discussed in the literature [18]. The distribution of current in Fig. 7 (II) is more uniform, but some areas of higher current density are observed. In Fig. 7 (III) the disorder of the oxide at the atomic scale is visible in the calculated current density.
Due to the dominance of the localized conduction channel in junction (I), we see a large weighting near zero in the histogram as many points in the xy plane have close to zero current density relative to the maximum current density. In junction (III), the typical insulating barrier, the distribution of current is comparatively uniform. The histogram for junction (II) lies between these two extremes where there are some regions of high current density but most of the cross-sectional area conducts only a small proportion of the current.

V. CONCLUSIONS
Fine control of the critical current is highly desirable for creating addressable qubits when fabricating devices containing tens or hundreds of Josephson junctions. In this work we study the interplay between the internal structure of the oxide and its electrical characteristics using a three-dimensional description of the junction. The material properties of the oxide layer in the Al−AlO x −Al junction are found to affect the calculated resistance-area product. We observe the exponential dependence between the thickness of the oxide barrier and the junction resistance as expected. Additionally we find that the junction resistance is exponentially dependent on the oxide density. The junction resistance also changes with the stoichiometry of the barrier with conduction in highly substoichiometric structures being dominated by metallic conduction channels. These local variations in current are of smaller length scales than those typically discussed in previous work.
To study how the electronic characteristics change due to local atomic structure we calculate the charge density and current density. In highly oxygen deficient structures conduction is dominated by metallic channels. However, even with more oxygen present, particular paths through the oxide contribute more to the current flow. This nonuniformity of the current distribution has important consequences for the influence of charged defects within the amorphous structure. Defects near dominant conduction paths are more likely to couple strongly to the current, contributing to the noise in the critical current I C .
Despite their widespread usage, Al−AlO x −Al junctions suffer universally from noise caused by two-level systems whose exact physical origin is an ongoing topic of interest [48]. Magnetic surface spins [49], delocalized atoms [50,51], and many other models have been proposed to explain the observed noise [52]. Understanding the physical origin of two-level defects and their impact on the electrical properties of junctions is key in achieving improvements and consistency in fabrication. The present work provides a framework for a better understanding of how the performance of a junction in a circuit relates to its atomic structure.
We have developed a computational approach for determining the electrical characteristics of Al−AlO x −Al junction models based on their atomistic structure. Using this technique allows us to study the role of junction morphology and composition in determining junction performance. An understanding of the exponential dependence of the junction resistance on barrier thickness and oxide density can be reached using relatively simple models. However, the relationship between the atomic structure and flow of current through the junction can only be fully understood with a complete three-dimensional treatment of the problem. A future application of the techniques developed in this work may be to consider transport in parallel through many junction models with varying thicknesses and densities. Developing computational modeling tools for atomistic simulation of electronic devices at the nanoscale will prove invaluable in optimizing their fabrication, leading to more reliable and reproducible nanoelectronics. shown where interactions at r < ±r c . are described by the Gaussian potential, while the Coulombic potential is used for r > r c . The markers are separated by a distance of 1/3 Å, which is representative of the discretization used.
Computational Infrastructure (NCI), which is supported by the Australian Government. The authors acknowledge the people of the Woi wurrung and Boon wurrung language groups of the eastern Kulin Nations on whose unceded lands we work. We also acknowledge the Ngunnawal people, the Traditional Custodians of the Australian Capital Territory where NCI is located. We respectfully acknowledge the Tra-ditional Custodians of the lands and waters across Australia and their Elders: past, present, and emerging.

Truncation of the Coulombic potential
In order to calculate the electrostatic potential V (x, y, z) inside the junction structures we use the Ewald summation method. This is a standard method for computing the electrostatic energy of a particular configuration of charges in a periodic system [53]. We implement a version of Ewald summation in which a modified version of the short range real space interaction is used. This is necessary because the finite difference approach used in our NEGF calculations becomes a poor approximation when confronted with the divergences arising from the Coulombic potential close to a charged particle. To account for this we replace the Coulombic potential for short-range interactions with a potential of a Gaussian form. The junction potential is calculated on an evenly spaced threedimensional grid using the coordinates and charges obtained from the molecular dynamics calculation.
We define a radius r c inside which the Gaussian-like description of the potential will be used and write down the function h(r) which combines the Coulombic and Gaussian components: The different potential profiles are shown in Fig. 8, where the red points indicate the potential used in our calculations. The energy scale is characteristic of the atomic sites in a our calculations where the magnitude of the potential of the order of tens of electron-volts. As we are interested in energies close to the Fermi energy μ 0 (∼1 eV), the application of the truncation still allows for the atomic structure inside the junction models to be reflected in the calculated electronic properties. By using the Gaussian-like potential to describe the short-range interactions we ensure that the potential varies smoothly throughout the junction structure and avoid the numerical instabilities of the Coulombic divergences.

Finite difference order
To choose an appropriate value for r c the transmission was calculated in three dimensions for a range of radii using both three-and five-point finite difference approximations. The transmission is plotted in the left-hand panels of Fig. 9 for values of r c = 1.0, 1.1, and 1.2 Å along with smoothing splines fitted with MATLAB. On the right-hand side the residuals are plotted showing the difference between the calculated transmission and the fitted spline. As the radial truncation increases the Coulombic divergences are smoothed out, which in turns affects the stability of the calculated transmission. The behavior of the residuals is more dependent on the value of r c than the order of the finite difference approximation. FIG. 10. The variance of the residuals r in Fig. 9 is plotted as a function of the truncation radius r c for three-and five-point finite difference approximations.
To obtain a single metric for the smoothness of the transmission calculation we calculate the variance of the residuals. Figure 10 shows the decrease in the variance of the calculated residuals var(r) as r c increases and also highlights that the choice of r c affects numerical accuracy more than changing the finite difference approximation. With the view to include as much of the physics around the atomic sites as possible we use a value of r c = 1.2 Å for the rest of the work.