Introduction

Superconducting quantum computers often use aluminium oxide tunnel junctions as Josephson junctions to introduce the required nonlinearity1,2,3,4,5,6,7. The tunnel barrier in such junctions is formed by a thin dielectric film of amorphous aluminium oxide (AlOx) 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 AlOx 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 devices7.

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–AlOx–Al tunnel junctions are most commonly produced using the double-angle evaporation process pioneered by Dolan9. 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 barrier10. Other fabrication methods which modify or even remove the standard Dolan bridge structure also employ a low-pressure oxidation step11,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).

Fig. 1: Overview of the junction formation simulation.
figure 1

An amorphous oxide is formed on an Al(100) surface as oxygen atoms are iteratively added to the simulation cell. Oxygen and aluminium atoms are shown as orange and grey spheres, respectively. a–c Views along and perpendicular to the z-axis after 30, 90, and 150 oxygen atoms have been introduced. d–f The growth of the second aluminium contact during the metallisation simulation. This calculation was performed with the S–M potential.

Full size image

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 reports14 (1.4–1.6 nm) to simulate the growth of Al–AlOx–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.

Results

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 resources15,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 atm6,20,21,22. Another way to develop models of Al–AlOx–Al junctions is to use a simulated annealing approach. In this case, a crystalline Al2O3 structure is simulated at a temperature above the melting point to generate disorder before being cooled to lock the atoms into a particular configuration23,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 Methods A which supports this approximation.

While experimentally molecular oxygen (O2)5,26, ozone (O3)27, and both charged (O+)28,29 and neutral atomic oxygen (O)30,31 can be used, fully describing the adsorption of oxygen molecules on aluminium surfaces is an open problem. There is an ongoing discussion in the literature regarding the magnitude of the O2 dissociation energy32,33,34,35 and it is not known whether charge or spin dynamics play the dominant role36. Campbell et al. simulated the oxidation of nanoclusters and found similar behaviour for both neutral atomic oxygen and O2 molecules excepting a change in the temperature at the surface15. 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 Al2O3 (3.97 g/cm3)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.

Fig. 2: Density analysis of a junction model.
figure 2

a The material density in the junction model as a function of position. b The final atomic structure of a Al–AlOx–Al junction model where oxygen and aluminium atoms are depicted as orange and grey spheres, respectively. This calculation was performed with the S–M potential.

Full size image

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-Carlo-based techniques38,39.

Fig. 3: Structural change in the oxide due to a deposited oxygen atom.
figure 3

a The path of the deposited oxygen atom is overlaid on the final atomic positions for this iteration. b The density profile in the structure over the duration of the simulation. The position of the deposited oxygen atom is shown as an orange line. c The path of the oxygen atom from a perspective above the structure. d The structure at the beginning of this simulation at t = 0 ps. e The structure at the end of this simulation at t = 15 ps. This calculation was performed with the S–M potential.

Full size image

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/AlOx 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.

Fig. 4: Structural properties during the oxidation simulation.
figure 4

a Density of the aluminium–aluminium oxide structure shown as a function of position and deposition progress. The dashed lines mark the points where structural changes occur as the surface oxide changes form or a layer of the aluminium substrate is consumed by the oxidation process. b Variation of the stoichiometric ratio between oxygen and aluminium over the course of the oxidation simulation. c Coordination of oxygen atoms about aluminium over the course of the oxidation simulation. This calculation was performed with the S–M potential.

Full size image

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 AlOx/vacuum interface is oxygen rich compared with the Al/AlOx 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 α-Al2O3, 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 distribution—was 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.

Fig. 5: The effect of temperature on the density (left) and stoichiometry (right) during oxidation simulations.
figure 5

Development of the structure during simulated oxidation at temperatures of (a, b) 77, (c, d) 300, (e, f) 370, and (g, h) 470 K. These calculations were performed with the S–M potential.

Full size image

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.

Fig. 6: The effect of temperature on the bond angles in the completed oxide structures.
figure 6

Bond angles in the surface oxide following simulated oxidation at temperatures of (a) 77, (b) 300, (c) 370, and (d) 470 K.

Full size image

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 molecules40. 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).

Fig. 7: Partial charges in the structure at varying stages of oxidation.
figure 7

The partial charges of oxygen and aluminium atoms are shown as orange and grey dots, respectively. The black lines show the net charge in the structure as a function of z. a, b A comparison between the S–M and ReaxFF potentials for a 1-nm-thick surface oxide. c–e The development of a charge gradient after continued oxidation using the S–M potential.

Full size image

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 calculations18. 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) surface41. The islands proceed to grow laterally and merge to cover the remaining exposed aluminium.

Fig. 8: Formation and filling of a hole in the oxide.
figure 8

Formation and filling of a hole in the oxide during the oxidation process, showing snapshots after the addition of (a) 60, (b) 120 and (c) 180 oxygen atoms. A hole is observed to form and then close during the oxidation of a 24 × 24 Å2 substrate. This calculation was performed with the S–M potential.

Full size image

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 Al2O3 between two metallic contacts of either pure aluminium or niobium and do not include any disorder in the oxide layer42,43. Junction models developed using a simulated annealing method provide a more accurate representation of the real oxide layer which is known to be amorphous25. Amorphous models have also been constructed by alternating the addition of O2 and Al layers to an Al substrate with intervening ab initio MD calculations at 300 K to progressively form the oxide layer44. We have recently reported transport properties of junction models formed with simulated annealing45.

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.

Fig. 9: Structural properties during the junction formation simulation.
figure 9

a Evolution of the density as the growth of a complete Al–AlOx–Al junction is simulated. The change from oxygen to aluminium deposition is indicated by the vertical dashed lines. b Evolution of the stoichiometry. c Evolution of the coordination. This calculation was performed with the S–M potential.

Full size image

The development of low-density regions at the AlOx/Al interfaces is visible in Fig. 9a. In Fig. 9b, we observe an oxygen-rich 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 AlOx/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 O2 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 bounds of these regions are found by first taking the positions of the outermost oxygen atoms OL and OR to determine the thickness of the oxide layer d = OR − OL. We then use four points along the z-axis (OL − d/4, OL + d/4, OR − d/4, OL + d/4) to determine the boundaries of three regions with width d/2: one in the centre of the oxide and one bridging each interface. Considering the regions in Fig. 10 by eye, this appears to be a good heuristic approach to defining the central and interfacial regions.

Fig. 10: Density profiles for various simulation parameters.
figure 10

a–d The density of Al/AlOx/Al junction models as a function of z (averaged over the x and y dimensions). Low densities are consistently observed at the Al–AlOx 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.

Full size image

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 Al2O3) 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.

Fig. 11: Histogram of the minimum and mean densities in the junction models.
figure 11

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/AlOx interfaces.

Full size image

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.

Fig. 12: Partial charges in the structure at varying stages of metallisation.
figure 12

The partial charges of oxygen and aluminium atoms are shown as orange and grey dots, respectively. The black lines show the net charge in the structure is as a function of z. a The initial oxidised aluminium surface. b–d Continued metallisation with the S–M potential. e Partial charges in the final junction structure formed with the ReaxFF potential.

Full size image

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.

Fig. 13: Bond angle comparison between the constructed and deposited aluminium regions.
figure 13

a Bond angles in the top aluminium contact grown by deposition. b 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.

Full size image

Discussion

The oxidation of aluminium is known to self-terminate when a thin amorphous oxide layer has been formed22,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 design47. 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 studied48,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:

  • Al(100): 7.87 ± 0.80 Å.

  • Al(111): 8.54 ± 0.56 Å.

A calculation was determined to have reached the limiting stage when 25 oxygens atoms in a row had been reflected from the surface without bonding. By comparison, limiting behaviour was never observed in the simulated oxidation of aluminium surfaces with the S–M potential. New oxygen atoms continued to bond to the surface as other oxygen atoms are displaced deeper into the structure. This proceeded until all of the aluminium atoms in the initial contact were incorporated into the growing aluminium oxide.

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,47. Zeng et al. measured the oxide thickness in this way at hundreds of positions for three Al–AlOx–Al samples14. Mean thicknesses of 1.66–1.88 nm are reported though the oxide thickness was measured to be as thin as 1.1 nm in some places and up to 2.2 nm thick in others.

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) substrates51,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.2 nm22. Another similar study reports self-limiting thicknesses of 0.49–1.36 nm on Al(100) and Al(111) surfaces for partial pressures from 1 × 10−5 Pa to 1.0 Pa53. Nguyen et al. find that Al(111) surfaces have slightly thicker oxides than Al(100) surfaces while the oxide thickness increases from 0.95 nm to 2.6 nm as the pressure changes from 4 × 10−5 to 4 × 10−3 Pa41.

Sankaranarayanan et al. simulate the oxidation of Al(100) with both O atoms and O2 molecules by maintaining a particular number density of oxygen around a aluminium crystal structure54. A self-limiting oxide thickness of 1.6 nm is reported as well as low densities at the AlOx/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 measurements50. On Al(100) and Al(111) substrates, oxides have been reported to be super-stoichiometric 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 oxide41. 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 Al2O3-pellets with an electron beam55. 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 γ-Al2O3 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 imaging38. Based on these measurements, an atomistic model of a possible tunnel barrier structure is then reconstructed which predicts oxygen deficiency at the Al/AlOx interfaces. Sankaranarayanan et al. also find higher oxygen concentrations at the AlOx/gas interface than the AlOx/Al interface in their simulations54.

Jeurgens et al. observe a change from amorphous to crystalline morphology in oxide layers at the temperature was raised from 573 to 773 K56. A change from amorphous to semi-crystalline “γ(-like)-Al2O3” structures was observed at between 400 and 550 K by Reichel et al. depending on the crystallographic orientation of the lower Al substrate51. 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 K8.

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 AlOx 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 methods57,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 oxides15,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 deficiency in the AlOx layer which has been predicted from reverse Monte-Carlo reconstruction of the atomic structure and, more recently, measured directly8,38.

Fritz et al. achieved epitaxial growth of an Al(111) layer on a clean Si(111) substrate49. In this case, thickness fluctuations in the AlOx 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–AlOx–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–AlOx–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 semi-crystalline 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–AlOx 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 reported61. 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 junctions5,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 junctions63. 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 design63.

Methods

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 atoms65. 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 Al2O3 crystals, we note that S–M and other similar charge transfer ionic potentials have previously been used to study the oxidation of aluminium nanoparticles15,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 thermostat69,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 LAMMPS71,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 field73,74 using parameters for aluminium and oxygen published by Hong and van Duin19. 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 experimentally76. 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 half-maximum (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 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.