Simulating the fabrication of aluminium oxide tunnel junctions

Aluminium oxide (AlO$_\mathrm{x}$) 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-AlO$_\mathrm{x}$-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.


I. 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][8] 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 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 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 emerge as oxygen and aluminium atoms are consecutively added to the surface. In Fig. 1(a)-(c) three stages of oxide growth are shown for a typical simulation. The aluminium surface is partly then completely covered with oxygen as more atoms are deposited. After the oxide layer is formed aluminium is added (metallisation) to completes the tri-layer junction structure [ Fig. 1(d)-(f)].
In Sec. II we outline the methodology we have used to simulate the oxidation process. Structural properties such as density and stoichiometry are studied over the course of the simulated oxidation. We also study 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. In Sec. III we 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.

A. Methodology
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 other 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 high vacuum limit. A derivation is provided in Appendix A which supports this approximation.
While experimentally molecular oxygen (O 2 ), 5,26 ozone (O 3 ), 27 and both charged (O + ) 28,29 and neutral atomic oxygen (O) 30,31 can be used, it is known that O 2 has a large dissociation energy which the empirical potentials have not been designed to describe. 32 Campbell et al. simulated the oxidation of nanoclusters and found similar behaviour for both neutral atomic oxygen and O 2 molecules excepting a change in the temperature at the surface. 15 We therefore deposit individual neutral oxygen atoms in the present work for simplicity (except where otherwise stated).
Most of the calculations in this work were performed with the General Lattice Utility Program (GULP). 33 This program includes an empirical potential developed by Streitz and Mintmire (S-M) which describes the interactions between aluminium and oxygen atoms. 34 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. The system is simulated as an NVT ensemble held at the chosen temperature by using the Nosè-Hoover thermostat. 35,36 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 Appendix B.
For comparison we have also performed MD with LAMMPS. 37,38 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 39,40 using parameters for aluminium and oxygen published by Hong and van Duin. 19 Aluminium substrates are prepared by creating supercells with the experimentally reported lattice constant of 4.041386Å. 41 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 constant 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.

B. Results
From the positions of the atoms we calculate the material density as a function of z (the direction of the oxide growth). Density is reported in units of ρ where ρ = 1 corresponds to the density of crystalline Al 2 O 3 (3.97 g cm −3 ). 42 This 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. Figure 2(a) shows an example of a density profile calculated in this way for the junction depicted in Fig. 2(b). 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. gen 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. 3(b). The spatially varying density at each time is calculated in the same way as for Fig. 2(a). 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 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 techniques. 43,44 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. 4(a) 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 occursuch as the crystalline-to-amorphous transition depicted in Fig. 3 (the second dashed line at N = 144 atoms)are marked with dashed lines.
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. 4 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. 4(b) we note that the AlOx/vacuum interface is oxygen rich compared with the Al/AlOx interface. Figure 4(c) 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.
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.
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 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. 45 We also examine how the charge distribution differs between the two empirical potentials we have used [compare Fig. 7(a) and (b)]. In both cases the net charge is neutral in the bulk of the oxide and tends to become negative at the metaloxide 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. 8(b)] which is filled in as the oxide continues to grow (see Supplemental 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. 46 The islands proceed to grow laterally and merge to cover the remaining exposed aluminium.

C. Discussion
The oxidation of aluminium is known to self-terminate when a thin amorphous oxide layer has been formed 22,47 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. 48 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. 49,50 In our calculations with ReaxFF we observed selflimiting 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, 48 Zeng et al. measured the ox-ide thickness in this way at hundreds of positions for three Al-AlO x -Al samples. 14 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). 51 Measurements of the limiting thickness made in this way find values in the range 5-10Å for Al(100) and Al(111) substrates. 52,53 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 nm. 22 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 Pa. 54 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 Pa. 46 Sankaranarayanan et al. simulate the oxidation of Al(100) with both O atoms and O 2 molecules by maintaining a particular number density of oxygen around a aluminium crystal structure. 55 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 measurements. 51 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. 22,54 The stoichiometry at the surface has been found to be lower than that in the centre of the oxide. 46 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). 47 Fritz et al. report stoichiometries for oxides grown using four different techniques: thermal oxidation with and without UV illumination, plasma oxidation, and physical vapor deposition achieved by heating Al 2 O 3 -pellets with an electron beam. 56 The stoichiometries were determined using STEM electron energy loss spectroscopy 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. 43 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 simulations. 55 Jeurgens et al. observe a change from amorphous to crystalline morphology in oxide layers at the temperature was raised from 573 to 773 K. 57 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. 58 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 amorphouscrystalline 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.

A. Methodology
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. 59,60 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 We have recently reported transport properties of junction models formed with simulated annealing. 61 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. 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 approximately 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. 62 The second aluminium electrode is grown until it is of a similar thickness to the initial aluminium contact region.
B. Results 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 AlOx/Al interfaces is visible in Fig. 9(a). In Fig. 9(b) 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 Supplemental 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 interfa-  Fig. 10 by eye this appears to be a good heuristic approach to defining the central and interfacial regions. Fig. 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) 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. 12(a)-(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. Fig. 13(d) 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. 13(a) is 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 13(b) and (c) show 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.

C. Discussion
The apparent density in the centre of the oxide increases as it grows. This effect can be seen in Fig. 4(a), 5, and 9(a). 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 63,64 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,43,65 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 metaloxide 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." 66 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. 9(b).
Fritz et al. achieved epitaxial growth of an Al(111) layer on a clean Si(111) substrate. 50 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.

IV. CONCLUSIONS
Using our novel 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 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-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. 67 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 an aluminium or niobium can be used as Josephson junctions. 5,68 Double-barrier junctions constructed with aluminium and aluminium oxides are used in magnetic tunnel junctions (MTJs). 69 Other materials are often used in MTJs such as in CoFeB-MgO-CoFeB junctions. 69 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. 69

ACKNOWLEDGMENTS
The authors acknowledge support of the Australian Research Council through grants DP140100375, CE170100026 (MJC), and CE170100039 (JSS). The authors also acknowledge useful discussions with P. Delsing, S. Fritz, J. Gale, D. Gerthsen, N. Katz, E. Olsson, J. Pekola, and L. Zeng. This research was undertaken with the assistance of resources from the National 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 Traditional Custodians of the lands and waters across Australia and their Elders: past, present, and emerging. events. We are interested in the flux of atoms striking a small surface region dA over time We start by considering an arbitrary velocity distribution for the particles in the gas. This can be written as The number of atoms inside an infinitesimal part of the velocity space can then be written as where N is the total number of particles. By assuming a spherical symmetry and transforming to spherical coordinates we can equivalently write where G(v) is defined only by the particle's speed v.
In real space, we consider only those atoms approaching the surface from a particular direction defined by θ and φ. We narrow this definition to include only particles at a given velocity v. These particles are at a distance v dt from the surface. Together these quantities define an infinitesimal volume dV (depicted in Fig. 14): where dA is the surface element. The number of atoms inside the volume dV can then be calculated from the concentration of atoms n: By substituting this in Eq. A4 with N → N V we obtain dN = n v cos θ dt dA G(v) v 2 dv sin θ dθ dφ (A7) Returning to the expression for the atomic flux (Eq. A1) we have Finally we are able to integrate this equation over a hemisphere to account for all the incoming particles on one side of the plane: x y z dA v dt v cos θ dt θ FIG. 14. A diagrammatic representation of the volume dV as defined in Eq. A5. Conceptually this is a slant cylinder containing all of the particles at a velocity v which will reach the surface dA after a given time dt. The angle of the cylinder relative to the surface is defined as θ.
where we have introduced the function f (v) We define the average velocity of a particle in the distribution f (v) asv which allows us to express the flux of the particles simply as Due to the low pressure state of the system we consider the gas to obey the ideal gas law equation. Under these conditions the Maxwell-Boltzmann velocity distribution is well suited to describing the statistics of the particles and the function f (v) has the form x 2n+1 e −px 2 dx = n! 2p n+1 ∀ p > 0. (A20) This is the average velocity of a particle of a given mass m in a Maxwell-Boltzmann distribution with the temperature T . We can also rewrite the concentration n in terms of the gas pressure using the ideal gas law Substituting these expressions forv and n in Eq. A16 we obtain the flux as a function of the pressure, temperature and the mass of the particles in the gas: Multiplying the flux by the area of the surface in the simulation allows us to estimate the number of atoms which strike the surface per unit time. Considering as an example the oxygen pressure of 1.33 × 10 −4 Pa reported by Jeurgens et al. and a temperature of 300 K gives a rate of approximately 52 surface interactions per second on a surface with area A 32.0 × 32.0Å 253 . This is equivalent to one of the largest structures we simulate. The limitation on the size of the surface arises from our choice of simulation package, the complexity of the empirical potential, and the available computing power on current supercomputing facilities. From our estimate one atom interacts with the surface about every 20 ms. As the total simulation time for one oxygen atom to be deposited on the surface is of the order of tens of picoseconds, it is a good approximation to consider these atom strikes as independent events.
Appendix B: Magnitude of the coupling strength for the Nosè-Hoover thermostat As the molecular dynamics simulations are performed at a constant temperature a thermostatting process must be included to describe heat transfer in and out of the simulation cell. We use the Nosè-Hoover thermostat algorithm in the simulation package GULP to perform simulations in the canonical (or NVT) ensemble. 36,71 This algorithm maintains a constant temperature in the simulation by coupling the particles to a fictional external thermal reservoir.
To understand how this works mathematically, consider the relationship between the instantaneous kinetic temperature T of the simulation and the velocities of the simulated particles. By the equipartition theorem we equate the thermal energy with the sum of the kinetic energies of the particles: In this expression d is the number of dimensions, N is the number of particles and k B is Boltzmann's constant.
The atomic mass and velocity of particle i are given by m i and v i respectively.
The thermostat introduces two new quantities: a coupling constant ν and a heat flow variable ξ which is related to the effective "mass" of the fictional heat reservoir. 35 The equations of motion are modified for particle i such that where F i is the force on particle i and T 0 is the target temperature for the thermostat, i.e. the desired temperature for the simulation. The relationship between ξ and ν arises from a consideration of the canonical NVT dynamics the thermostat is designed to reproduce. 36 When the instantaneous temperature T exceeds the target temperature, the value of Eq. B3 is positive. This means that the value of ξ will increase in Eq. B2, arresting the acceleration of the particle. By modifying the acceleration of the particles, heat is effectively added and removed from the system. The total energy of the system is given by the sum of the kinetic energy K, the potential energy U , and a contribution from the heat flow variable ξ: where q are the positions and p are the momenta of the particles. The thermostat coupling parameter ν defines the strength of the coupling, i.e. how quickly the particle velocities will respond to a temperature either above or below the target temperature. Here we establish a method for determining the appropriate strength of this coupling so as to accurately describe a canonical ensemble.
While the Nosè-Hoover thermostat equations were originally designed to reproduce the statistics of the canonical ensemble, Holian et al. have shown that using an incorrect value for the coupling can produce nonphysical oscillations in the temperature. 72 They present a number of physically motivated approaches for setting the thermostat coupling, one of which is to consider statistical fluctuations in the temperature over the course of the simulation. When the particles are strongly coupled to the external heat bath temperature fluctuations about the mean are small. Conversely, weak coupling allows the simulated system to act independently of the heat bath, leading to large fluctuations. The dependence of the temperature oscillations on ν is demonstrated in Fig. 15. In the weakly coupled case the temperature varies quite dramatically about the mean; as much as ±100 K. The strongly coupled case constrains the temperature more closely to the mean value.
For a canonical (or NVT) system Holian et al. give an equation which predicts the variance of the temperature: where d is the number of dimensions,T is the mean temperature and N is the number of particles. By calculating the variance in the temperature we obtain a metric which represents the magnitude of the fluctuations.
In order to study how the temperature variance changes as a function of ν, we simulate aluminium supercells with a range of coupling strengths. Three cubic supercells are simulated with side lengths of 16Å, 24Å, and 32Å. Each supercell is simulated with periodic boundary conditions for a total of approximately 0.8 ns solving the equations of motion every 1 fs. The total duration of ∼0.8 ns is reached by running a series of independent 10 ps calculations where the initial velocities of the particles in each calculation are randomised. Running a large number of independent simulations allows a reliable estimate to be obtained from the statistical analysis. The temperature variance was calculated for each 10.0 ps simulation before being averaged to generate the data on Fig. 16. An exponential fit to the data (R 2 = 0.9897) was calculated with MATLAB 73 σ 2 T (ν) = a exp (bν) (B6) where the constants were found to be a =610 and b = −49.6. Using the relationship between the coupling parameter and the variance of the temperature we are able to set the coupling to obtain the canonical temperature variance.
To demonstrate how ν is chosen for a given simulation with this expression, we consider as an example a typical simulation of 1000 particles where the mean temperature isT = 300 K. Using Eq. B5 we calculate the expected canonical variance σ 2 N V T as a function of the temperature and number of particles. In other words, in a simulation of 1000 particles at 300 K the canonical variance is reproduced by setting the thermostat coupling to ν = 0.0468. From a statistical analysis of the simulation temperature over a large period of time for different values of ν an exponential relationship is found between the temperature variance and the thermostat coupling. This allows us to express the coupling strength ν as a function of the variance (Eq. B6). Using the expression for the temperature variance in a canonical ensemble (Eq. B5) we are able to calculate the expected variance for any given system. 72 For our molecular dynamics simulations, we determine the expected variance as per Holian et al. and use Eq. B8 to set the value of ν appropriately.  16. The calculated variance in the temperature as the Nosè-Hoover thermostat coupling ν is varied. As expected the temperature oscillates more about the mean for small values of ν, where the thermostat coupling is weak. The dashed grey line shows the exponential function fitted to the data.