Formation of intrinsic and silicon defects in MoO3 under varied oxygen partial pressure and temperature conditions: an ab initio DFT investigation

Molybdenum trioxide (MoO3) is a promising material for energy conversion applications, including recent uses as a hole selective contact in silicon photovoltaic devices. The electrical and chemical properties of MoO3 are known to be strongly sensitive to the presence of intrinsic and extrinsic defects, which in turn are dependent on the fabrication route and processing conditions used to form the device layers. Of particular interest to this study were intrinsic defects comprising oxygen vacancies and extrinsic defects involving possible contaminant silicon atoms. Density functional theory simulations were used to predict defect concentrations as a function of processing temperature and oxygen partial pressure. A rigorous method is outlined to calculate defect formation energies for all intrinsic defects in MoO3, resolving conflicting information arising from previous studies. Brouwer diagrams were constructed and used to show that the charge neutral oxygen vacancy is dominant under most of the temperature and oxygen partial pressure conditions investigated. It was also shown that at commonly-used processing temperatures and oxygen partial pressures, silicon interstitials in MoO3 can introduce a spin-polarised defect state 0.5 eV above the MoO3 valence band maximum. Their concentration in MoO3 may reach 1.3 ppm with processing conditions of 700 K and 10 6 atm oxygen partial pressure, and this concentration is predicted to increase dramatically with higher temperatures and/or lower oxygen partial pressures. Our findings highlight the possibility of silicon contamination in hole-selective contact layers for silicon photovoltaic devices, with a potential increase in the parasitic absorption due to silicon defects in the contact layers reducing energy conversion efficiency.


Introduction
The worldwide search for efficient and cost-effective energy conversion and storage systems has resulted in the investigation of new uses for lesser known materials. Molybdenum trioxide and its sub-stoichiometric form (MoO 3Àx ; 0 < x < 1) is a material with proposed applications in batteries, 1-3 electrochemical capacitors, 4-8 sensors, 9,10 catalysis, [11][12][13] photo-catalysis, 14,15 electrochromism and photochromism. [16][17][18] It has also been found to be an effective hole transport layer in semiconductor devices [19][20][21] and a range of photovoltaic cells. [20][21][22] Sub-stoichiometric molybdenum trioxide was rst reported as a contact layer for organic, [23][24][25] kesterite 26 and perovskite [27][28][29][30] solar cells, however recently it has also found application as a hole-selective contact for higher efficiency silicon heterojunction cells. It has a large workfunction of 5.7 eV and wide band gap of 3.3 eV, which allows an electrical contact to achieve high hole selectivity without inducing signicant parasitic absorption. 22 Battaglia et al. 22 found that by replacing the p-type amorphous silicon hole contact layer with amorphous MoO 3Àx , the short circuit current density could be increased by 1.9 mA cm À2 , although a reduction in ll factor degraded the overall efficiency. The low ll factor was subsequently addressed by Geissbühler et al. 31 by limiting the processing temperature for devices to 130 C. This enabled them to achieve an energy conversion efficiency of 22.5% for silicon heterojunction cells. However, this achievement still falls short of the record efficiency of 24.7% 32 for a two-side contacted heterojunction cell, and the low processing temperature is incompatible with other manufacturing processes such as silver screen-printing. In order to overcome these problems, Cauduro et al. 33 recently proposed that controllable crystallisation of reactive-sputtered lms of MoO 3Àx can be used to form lms with a high workfunction and a highly ordered nanocrystalline sublayer structure that may enhance carrier mobility and allow for higher processing temperatures. Sook et al. 34 have also reported that crystalline MoO 3Àx has improved mobility and lower resistivity than amorphous MoO 3Àx , and has material properties that are highly dependent on preparation conditions. Consequently, in order to realise the full potential of MoO 3Àx as a low-cost hole transport layer for silicon solar cells, it is necessary to understand how preparation conditions affect the material's structural and electrical properties, as well as it's compatibility with silicon substrates.
Sub-stoichiometric molybdenum trioxide can be prepared by a number of techniques such as thermal oxidation, 35 spray pyrolysis, 36 or reactive sputtering. 34,37 The most common preparation method involves a two-step process, where MoO 3 powder is rst thermally-evaporated onto a substrate, and then the resulting lm is annealed at high temperature. 38,39 The conditions of preparation affects the defect chemistry, in particular by the formation of oxygen vacancies 35 and, if enough vacancies form, then the material can amorphize. 40 The stoichiometry of MoO 3Àx has been shown to inuence many important properties of the material, such as work function, 35 battery capacity, 41 band gap, 38,42,43 conductivity, 44 and optical reectance. 45 Point defects are also believed to create a defect band within the band gap of MoO 3Àx , which is thought to be crucial for selective carrier extraction. 22 A mechanistic understanding of the sensitivity of defect chemistry to deposition method and post-deposition processing is needed in order to tailor the properties of MoO 3Àx for device applications.
Crystalline a-MoO 3 has a structure consisting of layers of distorted octahedral, joined along the b-axis by van der Waals forces, as shown in Fig. 1. Three distinct oxygen sites exist in MoO 3 , with different coordination. The O1 oxygen atom, sometimes referred to as the terminal oxygen, is connected to only one Mo atom and extends along the b-axis into the van der Waals planes. The O2 oxygen atom, also known as the asymmetric oxygen, has two bonds to molybdenum atoms, of different length. The O3 atom is triply coordinated and has been called the symmetric oxygen, as the in-plane bonds with molybdenum atoms are the same length. Oxygen vacancies may form at each of these sites.
Density functional theory (DFT) calculations have been used to determine the oxygen vacancy formation energies at each site in previous studies, 46-49 however these studies have used different methodologies and yielded results that are highly inconsistent with each other. The formation energy of molybdenum vacancies was considered in one study, 48 but the formation energy of molybdenum and oxygen interstitials is still not known. In this study, a well-dened and consistent method is used to analyse the formation energy of all intrinsic defects, and then these are related to real-life defect concentrations at various experimental conditions of temperature and oxygen partial pressure through the construction of Brouwer diagrams.
Since many of the selective contact applications involve MoO 3Àx in contact with silicon surfaces, 20,22 silicon defects within bulk MoO 3 are also studied in this paper. A study by Anwar et al. 50 on Si-MoO 3 mixtures suggested that silicon atoms may penetrate MoO 3 and that these atoms could reduce the electrical conductivity. Previously, extrinsic defects in MoO 3 crystals such as Na, 47 Li 3,51 and H 52,53 have been studied with DFT, but so far DFT studies of the interaction between MoO 3 and silicon have only examined molybdenum oxides on silica surfaces. [54][55][56] In this paper, the possibility of silicon diffusion in crystalline MoO 3 is investigated and its effect on electronic properties is explored.

Computational details
Density functional theory calculations were performed with the projector augmented wave approach using the Vienna ab-initio simulation package (VASP). 57 Simulations in DFT require the choice of an exchange-correlation functional among other important parameters in order to balance accuracy with computational complexity. For this study, the PBE-GGA functional 58 was used with the addition of an on-site Coulomb correction (+U) of 6.3 eV in order to improve the localisation of d-electrons in molybdenum using the simplied rotationally invariant formalism of Dudarev et al. 59 in line with previous work on MoO 3 . 52,60,61 In order to account for the signicant effect of van der Waals forces in this material, the semi-empirical DFT-D3 energy correction by Grimme et al. 62 was used. While a number of different parameters have been used in the modelling of MoO 3 in the past, as is outlined in the Model verication subsection, the parameters chosen for this paper provide accurate results without excessively long simulation times.
The Brillouin zone was sampled with a 3 Â 1 Â 2 k-point mesh for the supercell of choice, 2 Â 2 Â 3. The choice of supercell shape and size is justied in the subsection Charge correction. Simulations were performed with a plane wave basis set energy cut-off of 500 eV, and a Gaussian smearing of the partial occupancies was applied with a width of 0.1 eV to speed up the simulations. The density of states were obtained by performing an additional single point calculation aer geometry optimisation using a denser k-point mesh of 6 Â 2 Â 4 kpoints and interpolating with the tetrahedron method with Blöchl corrections. The static ion-clamped dielectric tensor and elastic constants were calculated with density functional perturbation theory. For elastic constant calculations, both lattice and atoms were displaced in steps of AE0.005Å and AE0.01Å.

Formation energies
Defect formation energies were calculated according to eqn (1) where E defect and E host are the DFT energies of the defective and perfect supercells, the summation term refers to the chemical potential m of removed (positive) or added (negative) atoms, E VBM is the energy of an electron at the valence band minimum, m e is the Fermi level relative to E VBM , and E chgcor is the energy correction applied for charged defects to counter the selfinteraction across the periodic boundary conditions of nite supercells.

Chemical potentials
The chemical potentials used in eqn (1) have a large effect on defect formation energies, and must be determined carefully. The chemical potentials of molybdenum and oxygen atoms within the MoO 3 crystal are derived from the chemical potentials of the same atoms in reservoirs of BCC-Mo(s) crystal and O 2 (g) respectively. To overcome the well-known limitation of GGA-DFT inaccurately modelling oxygen gas molecules, 63 the method of Finnis et al. 64 is used to derive the chemical potential of oxygen in O 2 gas, m O 2 (g) , from the equation for standard formation energy of MoO 3 : where DG MoO3 f is the experimental standard formation energy of MoO 3 with a value of À668.079 kJ mol À1 , 65 and the chemical potentials m MoO 3 (s) and m Mo(s) are the DFT total energies of MoO 3 and BCC-Mo respectively, normalised per formula unit. The assumption is made that the chemical potentials of the MoO 3 and BCC-Mo(s) crystals are not strongly dependent on partial pressure and temperature, but that the potential of the oxygen gas will change according to the rigid-dumbbell ideal gas law: To test these assumptions, the expected free energy of formation at 700 K was calculated with eqn (2)-(4) to be À588 kJ mol À1 , close to the experimental value of À567 kJ mol À1 . 65 The total chemical potential of a formula unit of MoO 3 is the sum of the chemical potentials of its component parts inside the MoO 3 crystal: At equilibrium, the chemical potentials of atoms inside the crystal cannot be larger than that of the reservoirs outside. Within these bounding cases, the assumption is made that the chemical potential of oxygen in MoO 3 is equal to that of oxygen gas, which better describes the case of O-rich condition. Thus, the chemical potentials of oxygen and molybdenum in MoO 3 are: Where silicon impurities are considered, the chemical potential of silicon was calculated with the same procedure from DFT simulations of a-quartz SiO 2 , and assuming that it was in equilibrium with oxygen gas.

Charge correction
When charge is introduced into a supercell with periodic boundary conditions, an energy shi results from both the charge interaction between the repeating cells 63 and from the background uniform charge jellium introduced to compensate for this. 66 The magnitude of the required correction is a function of the screened Madelung potential, n scr M , calculated in an anisotropic medium with the following equation: 67 where g is a convergence parameter, R and G are the real and reciprocal lattice vectors and 3 is the dielectric tensor. Thus, in an anisotropic medium, the Madelung potential can be zero or negative for certain lattice vectors. 66 To demonstrate the effect of Madelung potential, the +1 charged O1 oxygen vacancy (V $ O using Kröger-Vink notation 68 ) was simulated in different supercell sizes and shapes. The resulting formation energy is shown as a function of the Madelung potential in Fig. 2. The uncorrected formation energy when using a 3 Â 1 Â 3 supercell is 0.25 eV lower than that from a 2 Â 2 Â 3 supercell. The amount of correction required is proportional to the square of the amount of charge introduced, 69 so for simulations involving high charges such as V $$ O or V 000000 Mo it is essential to have an accurate method of correcting for charged defects.
The 3 Â 1 Â 3 supercell has nearly equal simulation lengths in the a, b, and c directions, and has been used by most previous studies to model point defects. 47,48 However, when considering the anisotropic character of the dielectric constant, it is apparent that the 2 Â 2 Â 3 supercell has near-zero Madelung potential, yielding a very small Coulomb self-interaction energy. For the very highly charged defect of V 000000 Mo , a rst order approximation to the charge correction yields a charge correction less than 0.5 eV, and it has been shown that this is an overcorrection that should be scaled down further to correct for local screening, 70 indicating that the use of this supercell shape removes almost all of the uncertainties associated with the use of approximate charge correction schemes.
It is important to note that nite size effects may also arise from elastic self-interaction, 71,72 which are dependent on the anisotropic elastic constant tensor. Since the elastic tensor is distinctly different from the dielectric tensor, a compromise must be struck between elastic and electrostatic self-interaction when choosing the optimal supercell shape. Nevertheless, the moderate permittivity and comparatively large stiffness of MoO 3 , shown in Tables 1 and 2, means that the magnitude of the charge self-interaction energy is considerably larger than that of elastic self-interaction. Therefore, the current work was carried out using the 2 Â 2 Â 3 supercell.
The V O2 defects, which were stable in the 3 Â 1 Â 3 supercell, were not stable in this supercell, collapsing onto the O1 site. Akande et al. 48 report nearly identical formation energies for vacancies on the O2 and O1 sites, suggesting that this collapse was not a unique feature of our methodology. Nevertheless, the V O2 defect formation energy was calculated in a larger supercells , and it was consistently found to either collapse into V O1 or be much less favourable than V O1 for all charge states and in all supercell sizes investigated, thus the omission of these defects does not have a signicant effect on the predicted defect concentrations.

Brouwer diagram
Brouwer diagrams provide a means for identifying the concentration of defects in a material to the environment's oxygen partial pressure and temperature. The defect formation energies calculated from the DFT simulation were used to construct Brouwer diagrams with the aid of a defect analysis package. 76 For each value of p O 2 and T, there is only one Fermi level for which charge neutrality is conserved. Using an initial estimate for the Fermi level, the concentration of electrons and holes is calculated using the experimental bandgap according to Maxwell-Boltzmann statistics: where N c and N v are the conduction band and valence band density of states. Then, for each defect, the value of DE f i is calculated using eqn (1), (3) and (6). The concentrations of each defect i is given by: where m i is the multiplicity of the defect. These resulting trial concentrations are then used to calculate the total charge, q tot , from eqn (12) and if q tot is positive, the Fermi level is increased and if negative it is reduced. The process is iterated until the self-consistent Fermi level yielding charge neutrality was found for that p O 2 and T. The algorithm is then repeated for a range of temperatures and pressures to build Brouwer diagrams.

Model verication
The calculated lattice parameters were 3.84Å, 14.38Å and 3.77 A, for the a-axis, b-axis and c-axis respectively, which are all Fig. 2 Defect formation energy of V $ O as a function of Madelung potential, before (black solid) and after (red dashed) applying the approximate Lany and Zunger charge correction. 70 The labels indicate the dimensions of the supercell used (in terms of unit cell replicas).   Fig. 3a and b compare the axis length and bond lengths calculated in the current work with previous computational reports as a percentage of the experimental values. The largest deviation for the axis lengths is observed in the b-axis, where the inter-planar distance is dominated by van der Waals interactions. To overcome this, previous studies have xed the b-axis length to the experimental value 47,61,78 resulting in larger bandgaps and lattice parameters which are in closer agreement with experimental values. However, a constrained structure necessarily has internal strains; in turn, these are known to affect the formation energy of the defects considerably, as noted by Huang et al. 79 Therefore, all constraints were removed in this study when simulating the defect-free lattice. The other major differences in axis lengths and bond-lengths that are evident in Fig. 3a and b are caused by differing choices of exchange-correlation functional, on-site Coulomb correction (+U) value, and van der Waals forces treatment. This work used the D3 treatment by Grimme, 62 while earlier studies used the D2, 47,80 the optB88-vdW 46,81,82 or in the case of early studies, 53,60,61 no treatment at all. The most common level of theory was GGA + U with a U parameter of 6.3 eV, however there were studies that used smaller U values, 83 no U parameter 53 or the HSE06 hybrid functional. 48 The band gap of the simulated structure of this work was found to be 1.5 eV, which is lower than the experimental values that range from 3 to 3.3 eV. 22,36,38,84 The underestimation of the bandgap is a well-known issue with DFT simulations employing local and semilocal exchange correlation functionals, and the band gap predicted here is close to other values using the same functional, such as the 1.6 eV reported by Huang et al. 79 and the 1.7 eV reported by Peelaers et al. 85 for the PBE functional. Studies that xed the b-axis reported higher bandgaps of between 1.95 eV 61 and 2.4 eV, 47 in line with results showing that lattice distortion has a signicant effect on bandgap. 79 It has been shown that adding a +U correction only widens the computed gap by around 0.1 eV or so, 51,52 while the use of the HSE06 hybrid functional produces bandgaps consistent with experimental value, 47,48,85 at the cost of greatly increasing computational complexity.
With Bader analysis, 86 the Bader charges on the oxygen atoms were found to be À0.68e, À0.99e and À1.10e, for the O1, O2 and O3 atoms respectively, with a Bader charge on the Mo of 2.77e. These results match within 8% the recent works of Agarwal et al. 83 and Akande et al., 48 but have an up to 22% magnitude difference from older studies with no van der Waals treatment such as Coquet & Willock, 60 Sha et al., 53 and Lei & Chen, 52 indicating that van der Waals treatment has a noticeable effect on charge distributions in this material. Fig. 4 shows the wide discrepancies in oxygen vacancy formation energies in the literature, such as the 2.1 eV difference between the formation energy reported on the V Â O1 by Akande et al. 48 and that reported in Kim et al. 46 These variations are a considerably larger than what is acceptable for meaningful analysis of defect concentration, as there is an exponential dependence of concentration on formation energy. For example, using simple Boltzmann statistics at 300 K, a 1 eV drop in defect formation energy would multiply the predicted defect concentration by a factor of over 10 17 .
These large defect formation differences are caused by variations of the parameters in eqn (1). The differences in the formation energy of charge neutral vacancies can only be explained by the differing methods of determining the oxygen chemical potential and the variations in computational method already discussed. The differences in the Fermi level at which the cross-over between two charge states occur, are not affected by the oxygen chemical potential, but may be explained by other factors. For example, the use of the 3 Â 1 Â 3 supercell size without charge correction would make the charged states appear more favourable than they are, as was shown in Fig. 2, making them dominant for larger Fermi levels than are found in this work. An underestimation of the electron potential at the valence band minimum would have the same effect.

Intrinsic defects
There are several types of intrinsic point defects that could occur in MoO 3 , consisting of vacancies and interstitials of molybdenum and oxygen atoms. There are three different types of V O2 , corresponding to the three sites in a-MoO 3 with different coordination numbers. Three stable sites for were found for Mo i , shown in Fig. 5a, and three stable sites were found for O i , shown in Fig. 5b. The Wyckoff positions for these interstitials are presented in Table 3. Fig. 6 shows the defect formation energy at a temperature of 300 K and p O 2 of 0.2 atm for the most favourable charge state of each intrinsic defect as a function of Fermi level. Oxygen vacancies are most easily accommodated on the O1 site, and may take all three charge states but the charge neutral state is predominant for the majority of electron energies in the bandgap. The i6 site was found to be the most favourable oxygen interstitial, although it is only slightly more favourable than O i5 . Oxygen interstitials are predicted to mostly exhibit the 0 and À2 charge state, with the O i5 only exhibiting the charge neutral state. All oxygen interstitials exhibit highly unfavourable formation energies compared to the other defects. Molybdenum vacancies adopt all negative charge states except À4, and À6 while the dominant Mo i3 adopts all positive charge states except +4 and +5.
The only other study that considered V Mo was by Akande et al., 48 who predict formation energies in excess of 10 eV higher than the current work. Such large difference cannot be accounted for with the differences in DFT methodology, instead it is likely due to different treatment of the chemical potential.  Fig. 7a shows that at atmospheric pressure, signicant numbers of defects do not occur unless the temperature is elevated. The dominant defect by far is the charge neutral oxygen vacancy, with near-negligible amounts of V $ O and V 00000 Mo (two orders of magnitude less). These results match the experimental evidence that oxygen vacancies are the dominant defect in MoO 3 lms, 89,90 and that the amount of oxygen vacancies is dependent on the pressure 37 and temperature 45 used in sample preparation. The current Brouwer diagrams provide a further clarication of this temperature and pressure dependence.    Fig. 7b shows that if a high temperature anneal is performed under low pressure or low oxygen conditions, a large number of V Â O will form, along with some Mo i . The number of oxygen vacancies predicted at the very highest temperatures in Fig. 7b, would make the crystalline structure completely unstable. This matches with the experimental observation that initial thermal  evaporations at low pressures in the range of 10 À9 atm always results in sub-stoichiometric amorphous materials. 21,91,92 These results also are in agreement with the fact that annealing MoO 3 at temperatures of above 550 K in argon atmosphere results in high levels of oxygen vacancies but similar anneals in oxygen atmospheres do not. 36,40 It's worth noting that although the V Â O is dominant, the charged defects still play a crucial role in determining the Fermi energy that gives yields charge neutrality, which varies within the range of 1.1 eV to 1.9 eV for the pressures and temperatures considered. Fig. 7c and d show how changing the oxygen partial pressure affects the concentration of defects signicantly. The dominance of the V Â O in the range examined gives us a simple formula for determining the effect of varying oxygen partial pressure: by combining eqn (1), (3) and (11), it can be determined that the concentration of vacancies is inversely proportional to the square root of oxygen partial pressure, eqn (13). This concentration will be the dominant factor affected the stoichiometry of the material.
The accuracy of the estimated concentration of vacancies is hard to determine as there has been little quantication of defect concentrations in the experimental literature. Sian and Reddy 91 found a degree of sub-stoichiometry x of 0.01 for material annealed in air at 623 K, four orders of magnitude higher than that predicted in Fig. 7. However, the discrepancy may be largely attributed to the reported incomplete crystallisation in the experiment. The Brouwer diagrams also indicate that trace levels of Mo i defects may be observed at high temperature and low p O 2 , but no O i defects will occur under any of the preparation conditions shown.

Silicon defects
Recent applications of MoO 3 , particularly for photovoltaics, have involved MoO 3 junctions with crystalline or amorphous silicon, raising the possibility of silicon penetration into MoO 3 layers. To investigate this, DFT simulations of silicon defects within the MoO 3 crystalline structure were performed. Si i were  found to be stable on the same interstitial sites that were found for Mo i . Fig. 8 shows the formation energies for silicon defects. This shows that Si i exhibit amphoteric behaviour within the MoO 3 crystal. For all interstitials, the charge state switches gradually from +4 to +1 as the Fermi energy increases, before shiing to negative charges such as À2 and À4 near the conduction band minimum. The Si Mo was found to be less favourable than any of the Si i for most of the Fermi level range.
To investigate the accommodation of Si in MoO 3 under different processing conditions, Brouwer diagrams were constructed with silicon defects included, presented in Fig. 9. The Brouwer diagrams show that at high temperature and low p O 2 , silicon is highly soluble in MoO 3 , and that the preferential presence of Si i defects prevents the formation of V Mo , which were observed in trace amounts in the intrinsic Brouwer diagrams. The Si Mo defects do not appear under the conditions shown. At all temperatures and pressures, Si i are more abundant than Mo i , in fact at all pressures below atmospheric conditions, Si Â i are the second most prevalent defects aer V Â O , although the concentrations are not impactful until low p O 2 . Si interstitials are primarily charge neutral, with some +1 and À1 defects also occurring in roughly equal quantities at low p O 2 .
These results imply there might be some concern about silicon contamination. Fig. 10 shows the predicted silicon defect concentration as a function of both p O 2 and T, indicating that high amounts of contamination may occur under low p O 2 and high T conditions. The thermal evaporation process is oen performed at pressures below 10 À9 atm, 21,92 and a crystallisation anneal performed with a temperature of 700 K in an N 2 or Argon ambient with p O 2 of 10 À6 atm is predicted to produce a silicon defect concentration of 1.3 ppm (or 2.5 Â 10 17 cm À3 ). These conditions are similar to what is reported in the literature, 20,34 with p O 2 estimated on the assumption of 6N purity gasses, indicating that careful consideration of Si interfusion should be taken into account when preparing MoO 3Àx , either by limiting the temperature, avoiding low pressures, or by reducing the exposure time (provided the kinetics of Si diffusion are slower than that of MoO 3 evaporation). This result suggests that a high amount of silicon penetration would have a signicant impact on the electronic properties of the MoO 3 layer and hence the performance of any MoO 3Àx based device. Of particular concern, is the increased parasitic absorption in the MoO 3Àx contact layers which can directly reduce the energy conversion efficiency of silicon photovoltaic devices.

Conclusion
This study has explored the formation of both intrinsic (oxygen vacancies) and extrinsic (silicon interstitials) in crystalline a-MoO 3 under different preparation conditions using ab initio DFT. The results of the simulations are compared with previous studies and it was shown that, despite the axis and bond length properties of all previous works agreeing to within 9% of each other, the estimation of the defect formation energies of oxygen vacancies are highly inconsistent with each other. In order to rationalise these differences, a rigorous method for calculating formation energies was established, using an unconventional 2 Â 2 Â 3 supercell to minimise the electrostatic self-interaction between charged defects. All intrinsic defects were simulated, including the molybdenum and oxygen interstitials which had not been considered previously, and their formation energies were calculated as a function of the Fermi level. The calculated formation energies were then used to produce Brouwer diagrams that showed how the concentration of intrinsic defects varied with changing oxygen partial pressure and temperature. It was predicted (from the Brouwer diagrams) that the charge neutral oxygen vacancy is the most dominant intrinsic defect under the range of conditions considered, but other defects such as molybdenum interstitials could also be present at low p O 2 and high temperatures. Furthermore, it was  shown that the stoichiometry of the MoO 3Àx varied as a simple function of p O 2 and temperature. These ndings can help explain previous experimental results in terms of the effect of preparation conditions on material properties.
The methods used for intrinsic defects were then extended to silicon defects in MoO 3 . It was found that silicon is most readily incorporated as an interstitial species, where it introduces a spin-polarised defect state 0.5 eV above the MoO 3 valence band maximum and consequently reduces the bandgap of the metal oxide by 0.3 eV. The silicon defect can reach a concentration of 1.3 ppm at commonly-used processing conditions of 700 K and 10 À6 atm oxygen partial pressure, and signicantly higher concentrations at lower p O 2 and/or higher temperatures. This result suggests that silicon contamination of MoO 3Àx contact layers may represent a serious durability issue for silicon PV cells and modules and may impact device performance due to increased parasitic absorption.
The DFT method outlined in this paper could also be extended to other foreign defects such as hydrogen, which may be of particular interest for silicon solar cells as hydrogen can reduce interface state densities at a crystalline silicon surface. Further modelling could also determine the electronic effects of specic oxygen vacancy defect clusters in this material. The Brouwer diagrams developed through this study could be used to optimise fabrication processes to achieve a target stoichiometry and/or defect concentration. Since the defects have signicant effects on material properties, this tuning could have a large effect on the viability of crystalline MoO 3Àx contacts and hence the development and advancement of efficient and costeffective MoO 3Àx -based devices.

Conflicts of interest
There are no conicts of interest to declare.