Simulating the fabrication of aluminium oxide tunnel junctions

Aluminium oxide (AlOx) tunnel junctions are important components in a range of nanoelectric devices including superconducting qubits where they can be used as Josephson junctions. While many improvements in the reproducibility and reliability of qubits have been made possible through new circuit designs, there are still knowledge gaps in the relevant materials science. A better understanding of how fabrication conditions affect the density, uniformity, and elemental composition of the oxide barrier may lead to the development of lower noise and more reliable nanoelectronics and quantum computers. In this paper, we use molecular dynamics to develop models of Al–AlOx–Al junctions by iteratively growing the structures with sequential calculations. With this approach, we can see how the surface oxide grows and changes during the oxidation simulation. Dynamic processes such as the evolution of a charge gradient across the oxide, the formation of holes in the oxide layer, and changes between amorphous and semi-crystalline phases are observed. Our results are widely in agreement with previous work including reported oxide densities, self-limiting of the oxidation, and increased crystallinity as the simulation temperature is raised. The encapsulation of the oxide with metal evaporation is also studied atom by atom. Low density regions at the metal–oxide interfaces are a common feature in the final junction structures which persists for different oxidation parameters, empirical potentials, and crystal orientations of the aluminium substrate.


INTRODUCTION
Superconducting quantum computers often use aluminium oxide tunnel junctions as Josephson junctions to introduce the required nonlinearity [1][2][3][4][5][6][7] . The tunnel barrier in such junctions is formed by a thin dielectric film of amorphous aluminium oxide (AlO x ) which separates two metallic contacts. As interest has expanded in superconducting quantum computing architectures so too has the importance of clarifying the materials science which governs the formation and stability of thin AlO x films. Understanding the microscopic details of the oxide layer is a present focus for identifying and mitigating noise sources in a range of superconducting electronic devices 7 .
Fritz et al. have recently stated that "systematic studies of the influence of oxidation parameters on structural and nanochemical properties are rare up to now" 8 . Here we aim to address this gap in the literature from a computational perspective. The quantities reported, such as material density and stoichiometric Al:O ratio, and the conditions under which the oxide develops amorphous or crystalline features, can inform the fabrication of superconducting qubits with increased quality and reliability.
High-quality trilayer Al-AlO x -Al tunnel junctions are most commonly produced using the double-angle evaporation process pioneered by Dolan 9 . The aluminium layers are deposited through a lithographic mask at different angles to a substrate with an intervening low pressure oxidation which forms the oxide barrier 10 . Other fabrication methods which modify or even remove the standard Dolan bridge structure also employ a low-pressure oxidation step [11][12][13] .
In this study, we use molecular dynamics (MD) to explicitly model the low-pressure oxidation and aluminium evaporation processes. The structure of the oxide and junction emerges as oxygen and aluminium atoms are consecutively added to the surface. The aims of the present work are to demonstrate the method of oxidation simulation wherein events are modelled individually and elucidate trends in the oxidation process and in properties of the final junction structures. In Fig. 1a-c, three stages of oxide growth are shown for a typical simulation. The aluminium surface is partially then completely covered with oxygen as more atoms are deposited. After the oxide layer is formed aluminium is added (metallisation) to complete the tri-layer junction structure ( Fig. 1d-f).
Structural properties such as density and stoichiometry are studied over the course of the simulated oxidation. We also investigate how the charges on the atoms change as the dielectric barrier layer is formed and investigate the effect of temperature on the structure of the oxide. Our results are then discussed in comparison to computational and experimental results from the literature. We then consider aluminium deposition onto a formed oxide layer of similar thickness to experimental reports 14 (1.4-1.6 nm) to simulate the growth of Al-AlO x -Al junctions. We examine how the structure of these junctions changes as they grow and make comparisons between the two empirical potentials we have used in this work.

Oxidation of aluminium
One standard approach to studying oxidation computationally is to first create a region of aluminium surrounded by vacuum space before filling the vacuum with a nominal density of oxygen atoms or molecules. The system is then allowed to evolve until a stable oxide layer forms on the aluminium surface. In such studies, the oxygen gas density frequently corresponds to an unrealistically high pressure (~10-500 atm) in order to accelerate the dynamics and reduce the required computational resources [15][16][17][18][19] . By comparison, experimental junction fabrication is normally performed under high or ultra-high vacuum and partial oxygen pressures can vary over many orders of magnitude from 10 −12 to 10 −2 atm 6,20-22 .
Another way to develop models of Al-AlO x -Al junctions is to use a simulated annealing approach. In this case, a crystalline Al 2 O 3 structure is simulated at a temperature above the melting point to generate disorder before being cooled to lock the atoms into a particular configuration [23][24][25] .
As an alternative to artificially raising the gas pressure or creating disorder with annealing, we simulate the oxidation process directly by iteratively adding atoms to the surface. We use MD to model the approach and bonding of individual oxygen atoms to the bare aluminium surface. As we model individual atoms approaching the surface, we need not simulate the relatively long periods when no atoms are interacting with the surface. This results in a considerable reduction in the computational cost and corresponds to the case where the pressure is sufficiently low that the surface interactions can reasonably be treated as independent events. In this model, any sufficiently low pressure results in an equivalent simulation and result. A derivation is provided in Supplementary 15 . We therefore deposit individual neutral oxygen atoms in the present work for simplicity (except where otherwise stated).
From the positions of the atoms, we calculate the material density as a function of z (the direction of the oxide growth). Figure 2a shows an example of a density profile calculated in this way for the junction depicted in Fig. 2b. Density is reported in units of ρ where ρ = 1 corresponds to the density of crystalline Al 2 O 3 (3.97 g/cm 3 ) 37 . The aluminium contacts have a density of approximately ρ = 0.7 with some oscillations due to the alignment of the lattice planes perpendicular to the z-axis. There are notable drops in the density at the metal-oxide interfaces and the oxide region in the centre has a higher density than the contacts. Figure 3a-c shows the path of the newly added oxygen atom as it approaches then bonds to the surface. The evolution of the material density in the structure over the course of a 15-ps simulation is shown in Fig. 3b. The spatially varying density at each time is calculated in the same way as for Fig. 2a. We observe that the incoming oxygen atom-which is embedded 4-5 Å below the surface by the end of the simulation-seems to initiate the transition from a semi-crystalline (Fig. 3d) to an amorphous structure ( Fig. 3e) (see Supplementary Video 1). This suggests that a locally ordered region of the growing oxide may undergo this type of phase change due to a single oxygen atom disrupting the structure, though we note that the details of this effect may be modified by the finite size of the simulation cell. A more direct investigation of this transition may be possible using Monte-Carlobased techniques 38,39 .  While Fig. 3 shows how the density changes over the course of a single 15 ps simulation, we can also examine how the structure evolves as the oxide growth is simulated. In Fig. 4a, the density is calculated as a function of z at the end of each iteration, i.e. after a new oxygen has been added to the surface. We use this to understand the evolution of the density profile as oxygen atoms are consecutively deposited on 16 × 16 Å 2 Al(100) surface. There is a low density region at the lower Al/AlO x interface which is persistent and moves down in z as more aluminium is incorporated into the growing oxide layer. The points where significant structural changes occur are marked with vertical dashed lines, e.g. the second dashed line at N = 144 atoms in Fig. 4 corresponds to the crystalline-to-amorphous transition depicted in Fig. 3.
In the same way we determine the density as a function of z, we calculate the spatial variation of the stoichiometry (O:Al ratio) and coordination number for each iteration (Fig. 4b, c). The abrupt structural changes visible in the density data can also be seen in the stoichiometry and coordination. Coordination numbers of aluminium atoms are calculated by counting the number of oxygen atoms within 2.4 Å. This distance corresponds to a position in the radial distribution function g(r) between the first and second peaks.
In Fig. 4b, we note that the AlO x /vacuum interface is oxygen rich compared with the Al/AlO x interface. Figure 4c shows that more highly coordinated aluminium atoms tend to be towards the surface. This is consistent with the analysis of the stoichiometry which shows that more oxygen atoms are available for bonding closer to the surface. In crystalline α-Al 2 O 3 , the maximum coordination value is six which is reflected here. A high density of oxygen atoms accumulating at the surface precedes the marked structural changes.
Oxidation calculations were performed on an Al(111) surface with the temperature of the thermostat set to correspond to experimental values of interest: liquid nitrogen cooled (77 K), room temperature (300 K), heated to 100°C (370 K), and heated to 200°C (470 K). The temperature of the oxygen gas-used to generate the velocities from a Maxwell-Boltzmann distributionwas 300 K in all cases. This is a computational approximation to experimental conditions where it is far more likely to be able to change the temperature of the substrate than the temperature of the gas introduced to the chamber.
The evolution of the density and stoichiometry in the growing oxides is shown in Fig. 5. Here we can again see the step-like way in which the aluminium substrate is converted into surface oxide. This is most clearly visible in the low temperature calculations. The calculations proceed in general without the abrupt structural changes such as those in Fig. 4 though some such features can be seen in panels c and g. There appears to be minimal difference between the 77 K and 370 K calculations other than the thermal atomic motion limiting the clarity of the calculated density. By examining the bond angles in the oxide (Fig. 6) after 200 atoms have been deposited, we can see a structural difference which is not evident in the density or stoichiometry. The bond angle analysis shows strong peaks for temperatures of 300 K and 370 K indicating the presence of semi-crystalline structures in the oxide. The high temperature calculation (470 K) has the same crystalline peaks which have been broadened slightly by the thermal noise. By comparison, the low temperature calculation (77 K) is significantly more amorphous. Figure 7 shows how the distribution of charge in the system changes at different stages of oxide growth at 300 K. The continued oxidation of the surface with the Streitz and Mintmire (S-M) potential gives rise to a charge gradient across the oxide. This is in agreement with the empirical understanding of oxidation. Mott-Cabrera oxidation theory is predicated on the effect of such a charge gradient on incoming oxygen atoms and molecules 40 . We also examine how the charge distribution differs between the two empirical potentials we have used (compare Fig. 7a, b). In both cases, the net charge is neutral in the bulk of the oxide and tends to become negative at the metal-oxide interface, though the charge separation is smaller in magnitude by around a factor of two with the ReaxFF potential. We are unable to compare the charges at the later stages of oxidation as the ReaxFF potential qualitatively reproduces the natural termination of the process at a limiting thickness (see following section).
The iterative method by which we form the oxide layer allows for a detailed study of the dynamics. Figure 8 shows clustering of oxygen atoms on the aluminium surface forming a hole (evident in Fig. 8b) which is filled in as the oxide continues to grow (see Supplementary Video 2). Holes which form and close in this way have previously been observed in MD calculations 18 . Experimentally, Nguyen et al. observed the formation of islands at lattice shelves (terraces) on the aluminium surfaces by making a series of time-resolved observations of the growth of oxides on pristine Al (100) and Al(111) surface 41 . The islands proceed to grow laterally and merge to cover the remaining exposed aluminium.

Junction formation
Relatively few attempts have been made to construct complete ab initio junction models, and those that exist are mostly limited by the high computational cost of density functional theory (DFT) calculations. These models have been created by placing a stoichiometric layer of Al 2 O 3 between two metallic contacts of either pure aluminium or niobium and do not include any disorder in the oxide layer 42,43 . Junction models developed using a simulated annealing method provide a more accurate representation of the real oxide layer which is known to be amorphous 25 . Amorphous models have also been constructed by alternating the addition of O 2 and Al layers to an Al substrate with intervening ab initio MD calculations at 300 K to progressively form the oxide layer 44 . We have recently reported transport properties of junction models formed with simulated annealing 45 .
When working with the S-M potential, rather than continuing the oxidation indefinitely, we create junction models by beginning to deposit aluminium on the surface when the oxide reaches the desired thickness (~1.4-1.6 nm). The oxides grown with ReaxFF self-limit at a given thickness after which we start depositing aluminium. Figure 9 shows how the density, stoichiometry, and coordination evolve during the creation of a deposited junction. Oxygen atoms are consecutively added to a 16 × 16 Å 2 Al(100) surface until the oxide layer reaches a thickness of 1.4 nm. After the oxide layer is formed, aluminium is deposited to form the second electrode of the junction structure. The vertical dashed lines indicate the point at which this change from oxygen to aluminium deposition takes place.
The development of low-density regions at the AlO x /Al interfaces is visible in Fig. 9a. In Fig. 9b, we observe an oxygenrich surface at the end of the oxidation process in agreement with Fig. 4. New aluminium atoms quickly bond to this surface oxygen and the stoichiometries at both AlO x /Al interfaces become equivalent (see Supplementary Video 3).
The spatial variation in the material density for four finished junction models is shown in Fig. 10. The structure in panel b was formed by oxidising the surface with O 2 molecules rather than single O atoms and showed no discernible difference in any of our analyses. Low densities are again observed at the interfaces between the contacts and the oxide. This is a common feature in our analysis of the density as a function of position regardless of the crystal orientation, temperature of the aluminium substrate, or the empirical potential used. The interfacial and central regions of the structure are shaded in blue and orange, respectively. The  Fig. 10 by eye, this appears to be a good heuristic approach to defining the central and interfacial regions. Figure 11 shows histograms of the minimum and central density for a range of junction models formed with both the S-M and ReaxFF potentials. The minimum densities are determined from the lowest density value in the blue shaded interfacial regions in Fig. 10 and the centre density is the mean of the orange shaded region at the centre. Both potentials predict a reduced density at the interface (ρ = 0.56-0.58, i.e. 56-58% of the density of crystalline Al 2 O 3 ) which is a persistent feature across all simulations. Junctions deposited with the S-M potential have a higher density in the centre of the oxide.
Looking at the partial charges on the atoms as we add aluminium to the oxide surface (Fig. 12a-d), we can see the shape of the net charge become negative at the interfaces and neutral in the centre of the barrier. The same profile is observed for a junction created with ReaxFF (Fig. 12e) though the barrier is significantly thinner due to the self-limiting of the oxidation calculation. As in Fig. 7 we observe that the magnitude of the charge separation is reduced relative to the S-M results.
We use the bond angles in the aluminium contacts as a measure of the crystallinity of the structures. Figure 13d shows the bond angles calculated for the Al(100) and Al(111) substrates we generate at the beginning of the oxidation simulation after thermalisation at 300 K. The data shown in Fig. 13a are for the deposited aluminium contacts demonstrating that the crystal structure forms naturally as a result of the atomic interactions described by the empirical potential. Figure 13b and c shows the junction structures as grown on Al(100) and Al(111) substrates, respectively. The ordering of the atomic layers is partially evident in these images, however it is somewhat obscured as the orientation of the top contact is not aligned with the substrate direction.

DISCUSSION
The oxidation of aluminium is known to self-terminate when a thin amorphous oxide layer has been formed 22,46 and, as the magnitude of the tunnelling current in Josephson junctions is exponentially dependent on the thickness of the oxide layer, the factors which affect the self-limiting thickness are important considerations for device design 47 . In order to optimise processes for device applications, the effect on the uniformity and morphology of the barrier obtained by heating or cooling the aluminium crystal substrate, using single-crystal substrates of different orientations, or varying the oxidation pressure has been studied 48,49 .
In our calculations with ReaxFF, we observed self-limiting behaviour on both Al(100) and Al(111) surfaces (averaged over eight simulations of each crystal orientation) at thicknesses of:  In recent studies, the oxide thicknesses have been measured directly by taking images of the structure at nanometre scales with scanning transmission electron microscopy (STEM) 14 As an alternative to measuring the barrier thickness in STEM images, an estimate of the average thickness over a large area can be made by comparing the relative intensities of the aluminium metal and aluminium oxide signals obtained from x-ray photoemission spectroscopy (XPS) 50 . Measurements of the limiting thickness made in this way find values in the range 5-10 Å for Al (100) and Al(111) substrates 51,52 . For Al(111) surfaces oxidised at room temperature over a wide range of pressures between 1 × 10 −6 Pa and 650 Pa, the self-limiting thickness was found to increase monotonically from 0.2 to 1.  54 . A self-limiting oxide thickness of 1.6 nm is reported as well as low densities at the AlO x /Al interfaces. Due to the manner in which we approximate the deposition process, our work is most reasonably compared to the lowest pressure experimental reports which also produce the thinnest oxide layers. The thicknesses we report here are of the same order as existing experimental and computational reports although, based on the most recent studies, they are likely to be lower than the true values.
Many studies which report thickness also investigate the composition of the oxide layer (i.e. ratio of oxygen to aluminium) which can also be estimated from XPS measurements 50 . On Al(100) and Al(111) substrates, oxides have been reported to be superstoichiometric with final O:Al ratios of 1. 6-1.7 (refs. 22,53 ). The stoichiometry at the surface has been found to be lower than that in the centre of the oxide 41 . The overall composition of the oxide on Al(431) substrates is reported to be stoichiometric (O:Al = 1.5) whereas the surface is highly substoichiometric (O:Al = 0.3-0.7) 46 .
Fritz et al. report stoichiometries for oxides grown using four different techniques: thermal oxidation with and without UV illumination, plasma oxidation, and physical vapour deposition achieved by heating Al 2 O 3 -pellets with an electron beam 55 . The stoichiometries were determined using STEM electron energy loss spectroscopy (EELS) and are in the range 1.1-1.3 in the amorphous oxide regions except for the thermally oxidised sample without UV illumination which has a reported stoichimetry of 0.5. In some cases, nanocrystals of either Al or stoichiometric γ-Al 2 O 3 are formed.
High oxygen concentrations at the surface, such as those we report, are in agreement with other computational work on the topic. Zeng et al. investigated the microstructure of an oxide barrier with STEM imaging 38 . Based on these measurements, an atomistic model of a possible tunnel barrier structure is then reconstructed which predicts oxygen deficiency at the Al/AlO x interfaces. Sankaranarayanan et al. also find higher oxygen concentrations at the AlO x /gas interface than the AlO x /Al interface in their simulations 54 .
Jeurgens et al. observe a change from amorphous to crystalline morphology in oxide layers at the temperature was raised from 573 to 773 K 56 . A change from amorphous to semi-crystalline "γ (-like)-Al 2 O 3 " structures was observed at between 400 and 550 K by Reichel et al. depending on the crystallographic orientation of the lower Al substrate 51 . A similar trend of temperature dependence has been observed in recent work which reported features in EELS spectra corresponding to amorphous structures at 343 K and crystallinity at 523 K 8 .
We observe crystalline features in the bond angle distribution even at room temperature (300 K). These features are not present at 77 K which suggests that the temperature at which the amorphous-crystalline transition takes place may be reduced by the periodic boundary conditions in the simulation. A more detailed future study focusing solely on temperature effects and using a range of substrate sizes would be ideal to further understand this effect.
We can also find that the apparent density in the centre of the oxide increases as it grows, which we can see in Figs. 4a, 5, and 9a. From Fig. 11, we see that the two empirical potentials give similar densities at the metal-oxide interface while S-M predicts a higher density at the centre of the oxide than ReaxFF, although part of this discrepancy may be caused by the reduced oxide thickness in the ReaxFF junctions. The density of AlO x barriers formed with thermal oxidation is not widely reported in the literature (possibly due to the difficulty of measuring the nm-thick layer). Studies which use different experimental methods 57,58 to deposit thicker layers (1 μm) report densities in a wide range from ρ = 0.58 to 0.95. Oxide densities reported in simulations of thin film oxides 15,23,24,38,59 lie within a narrower range of ρ = 0.73-0.88. Spatial variation of the density is evident in many of our results with a pronounced reduction at the metal-oxide interface. This is in agreement with Auger analysis by Evangelisti et al. which "suggests density variations across the oxide layer, with lower densities near the surface and the metal-oxide interface" 60 . The authors also note that they measured minimal variation in the stoichiometry across the thickness of the oxide which is in agreement with our results in Fig. 9b. We also observe oxygen  49 . In this case, thickness fluctuations in the AlO x are minimised, and matching of the crystallographic orientation between the lower and upper aluminium layers is observed. In the present work, we observe crystallinity in both aluminium layers but no alignment between the top and bottom contacts. It would be an interesting extension to perform junction formation calculations as a function of temperature and the thickness of the oxide layer to increase our understanding of how this information is transferred across the oxide layer.
Using our iterative approach to oxide growth, we have created Al-AlO x -Al junction models with both the S-M and the ReaxFF potentials. A key difference in the behaviour of the potentials is that ReaxFF qualitatively reproduces the self-limiting behaviour which is observed experimentally. The final densities of the oxides formed with ReaxFF are closer to the mean of the experimental reports though the densities in the S-M models are still within the experimental range. Without more accurate reports of the oxide density for direct comparison, it is difficult to comment on the reliability of the empirical potential in faithfully reproducing the physics of the oxide formation. Making a comparison in relative terms, ReaxFF is a more modern potential which qualitatively reproduces results closer to experimental reports. It is possible that a reparameterisation of the force field for the oxidation of aluminium surfaces rather than nanoclusters may further improve the accuracy of the results.
In general, ab initio models of Al-AlO x -Al junctions are difficult to develop due to the inherently amorphous oxide layer. The iterative approach we adopt in building the oxide layer atom by atom allows us to see dynamic changes in the structure that would be missed when creating oxide models with simulated annealing. The formation and closing of holes in the oxide, the transition of surface oxide between amorphous and semicrystalline configurations, and the development of a charge gradient are all examples of these observations. We believe this type of simulation to be a promising approach as many results in the present work-such as self-limiting oxidation, the trend of temperature dependence of the oxide crystallinity, the reduced density at Al-AlO x interfaces, and the crystallisation of the deposited aluminium contacts-are in line with experimental reports.
We also note that the iterative deposition approach is easily adaptable to study other thin film deposition processes, provided that an empirical potential is used which appropriately describes the interactions between the different atomic species. For example, experimental evidence of an amorphous interface layer consisting of Al, Si, and O between the bottom aluminium contact and the silicon substrate has been reported 61 . It may be possible to observe the development of this interface layer by including the silicon substrate in the simulation and performing an iterative oxidation calculation.
The growth of ultra-thin oxide layers is relevant to the manufacturing of many different devices. Single-barrier junctions which use superconductors such as aluminium or niobium can be used as Josephson junctions 5,62 . Double-barrier junctions constructed with aluminium and aluminium oxides are used in magnetic tunnel junctions (MTJs) 63 . Other materials are often used in MTJs such as in CoFeB-MgO-CoFeB junctions 63 . While some concepts for creating magneto-resistive random access memory (MRAM) have even more exotic geometries, all of these devices make use of at least one thin oxide layer in their design 63 .

Molecular dynamics
Most of the calculations in this work were performed with the General Lattice Utility Program (GULP) 64 . This program includes an empirical potential developed by Streitz and Mintmire (S-M) which describes the interactions between aluminium and oxygen atoms 65 . As the aluminium oxide is an ionic material, charge transfer between atoms is an important component of this potential. We solve the equations of motion and the distribution of charge every 1 fs. While the S-M potential was parameterised to describe bulk Al and Al 2 O 3 crystals, we note that S-M and other similar charge transfer ionic potentials have previously been used to study the oxidation of aluminium nanoparticles [15][16][17][18]66,67 . Additionally, a recent assessment of empirical potentials for the description of alumina nanoclusters found that S-M compares favourably to DFT calculations at intermediate sizes (1-4 nm) 68 .
The system is simulated as an NVT ensemble (i.e. with constant particle number N, volume V, and temperature T) and held at the chosen temperature by using the Nosè-Hoover thermostat 69,70 . The coupling between the system and the heat bath is an important parameter of the thermostat. A detailed description of how we choose this value for various simulations is provided in Supplementary Methods B.
For comparison, we have also performed MD with LAMMPS 71,72 . The parameters of these simulations (temperature, timestep, duration, etc.) are identical to those we use in GULP. Interactions between atoms-including the charge equilibration processes-are described by the ReaxFF force field 73,74 using parameters for aluminium and oxygen published by Hong and van Duin 19 . These parameters were found to reproduce the dependence of the limiting oxide thickness on temperature reported by Jeurgens et al. 52 .
Aluminium substrates are prepared by creating supercells with the experimentally reported lattice constant of 4.041386 Å 75 . An optimisation of the geometry is performed in the MD software (either GULP or LAMMPS) where the atomic positions and the dimensions of the supercell are allowed to change to find the lowest energy configuration. Vacuum is then added to increase the z dimension of the supercell to 20 nm before a second optimisation which allows for aluminium layers to expand at the metal-vacuum interfaces. This is the direction in which the oxide will grow as oxygen atoms are deposited. The lattice constants of substrates optimised with S-M and ReaxFF potentials differ very slightly, however both are within ±0.2% of the experimental value.
There are two steps in our methodology for each atom added to the surface, each involving a different MD simulation. First the atom is positioned 2.4 nm from the existing surface with a randomised position in x and y. The initial velocity is obtained from the Maxwell-Boltzmann distribution and constrained to be directed towards the surface. The system is then allowed to evolve for 15 ps. If the atom has been reflected from the surface or is not bonded for any other reason, the iteration is discarded and a new atom added. If the atom is bonded, then a relaxation calculation is performed where the system equilibrates for 2 ps. This technique serves to separate the individual atomic depositions so that they can be considered to be independent events. In order to simulate the addition of enough atoms to form the surface oxide and the second electrode-approximately 300 atoms for a substrate with a side length of  x = y = 16 Å-a total simulation time of 5 ns or more is required. While this forms a substantial computational challenge, it is many orders of magnitude shorter than the minutes of oxidation in experiments. Calculating the deposition of 300 atoms corresponds to~100-200 h in real time using 2 × 24-core Intel Xeon Cascade Lake Platinum 8274 3.2 GHz CPUs.
The methodology for aluminium deposition is the same as for the oxidation, excepting that the velocities are selected from a normal distribution with a mean of~600 m/s and a standard deviation of 20 m/s. These values are representative of the evaporation method of thin-film deposition which is used experimentally 76 . The second aluminium electrode is grown until it is of a similar thickness to the initial aluminium contact region.

Calculation of structural properties
The calculation of structural properties as a function of position (along the z axis) is achieved by taking a Gaussian window with a full-width halfmaximum (FWHM) of 2.4 Å (σ ≃ 1 Å) and moving it along z in increments of 0.05 Å. The FWHM of the Gaussian is taken from the position in the radial distribution function g(r) between the first and second peaks. This same distance is used to determine coordination numbers and reflects the average distance between nearest-neighbour atoms. Based on the position of each atom in space relative to the window, a weighting between 0 and Fig. 10 Density profiles for various simulation parameters. a-d The density of Al/AlO x /Al junction models as a function of z (averaged over the x and y dimensions). Low densities are consistently observed at the Al-AlO x interfaces. The blue and orange shaded regions are used to determine the minimum and mean densities, respectively, for Fig. 11. This calculation was performed with the S-M potential. Fig. 11 Histogram of the minimum and mean densities in the junction models. The distribution of minimum and mean densities in junction models deposited on 16 × 16 Å 2 Al(100) substrates with different potentials. The minimum density is the minimum value in the blue shaded regions in Fig. 10. The central density is the mean value of the orange shaded region. As shown in Fig. 10, the minimum densities occur at the Al/AlO x interfaces. 1 is thereby allocated. The calculated density at a given position is then given by the average of the weighted atomic masses for the different atomic species. The stoichiometry at a position is the weighted ratio of oxygen to aluminium atoms in the Gaussian window. The coordination number is likewise calculated for all atoms before the mean of the weighted values is assigned to the location in z.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.

CODE AVAILABILITY
GULP is freely available to download from Curtin University (https://gulp.curtin.edu. au/gulp/). LAMMPS is freely available under the GNU General Public License from https://lammps.sandia.gov/. The code used as a wrapper for these packages to perform the iterative deposition calculations was written by M.J.C. and is available upon reasonable request. The bond angle analysis code was written by G.O. and forms a part of NCPac which is currently under development in alpha testing at DATA61, CSIRO (https://research.csiro.au/mst/ncpac/). The junction structure as grown on Al(100). c The junction structure as grown on Al(111). d Bond angles in the bottom 24 × 24 Å 2 substrates. This calculation was performed with the ReaxFF potential.