Magnetic stability of oxygen defects on the SiO₂ surface

Defects in silica are known to degrade the performance of many devices, including conventional electronic¹ and quantum^2,3 circuits. Defects are introduced into the silica by both device processing and operation,^4–8 such that the concentration of defects depends signssificantly on the fabrication process. For example, in the case of superconducting circuits, noise attributed to dielectric loss can be reduced by employing a thermal SiO₂ process instead of plasma-enhanced chemical vapor deposition,⁹ which creates a significantly higher concentration of E′ centers in the film.¹⁰ The E′ center is one of the most extensively studied defects in silica and α-quartz^11,12 and is generally defined as a paramagnetic three coordinated Si atom (III-Si). Significant research efforts over the last five decades have focused on E′ centers in irradiated silica, since irradiation can generate both oxygen deficiency centers (ODC),⁶ which are precursors to E′ centers, and also E′ centers directly.^13,14

First-principles calculations have enabled the assignment of atomistic defect structures to the measured electron paramagnetic resonance (EPR) and optical signatures of a wide variety of E′ centers and ODCs in silica and α-quartz.^15–17 Calculations of the EPR of E′ centers remains an active field.¹⁸ Most calculations have focused on bulk silica, but similar defects exist on silica and α-quartz surfaces.^19–21 Surface defects are particularly relevant in circuit processing because they can be buried at interfaces by overlaying other materials on top of the silica, and because their concentrations and properties can be controlled through the fabrication process.

Coupling of circuit components to low lying energy levels introduced by defects in amorphous materials, like silica, can result in random fluctuations in quantities relevant to device performance such as charge and flux states.^22,23 Models of this noise typically consist of a bath of (possibly correlated) two level quantum systems (TLS). In fact a simple model consisting of a bath of TLS with a distribution of energy barriers between the two quantum states can reproduce much of the observed noise spectral density in superconducting quantum circuits.^9,24–27 The microscopic origins of the TLS, however, remain speculative.^28–31 TLS arising from paramagnetic defects in particular are believed to be a significant source of decoherence and dissipation in quantum circuits.³² 

In this study we focus on the magnetic stability of E′ centers on the surface of silica as motivated by previous research on superconducting quantum interference devices (SQUIDs) and superconducting qubits (SQs).^32–39 These works proposed that magnetic flux fluctuations are caused by unpaired electrons hopping between defect sites⁴⁰ or flipping spin.²⁸ Both of these studies suggested that paramagnetic III-Si are a possible source of hopping or flipping electrons, but do not explicitly compute electronic or magnetic properties of these defects in the bulk or on the surface of silica. Subsequent experimental efforts that investigated the dependence of the noise spectral density on device geometry^33,35 suggested the source of 1/f noise originates either on the substrate surface, at the interface between the superconducting material and the substrate, or on the superconductor surface’s native oxide.

The development of scalable solid-state quantum computers based on certain superconducting circuit designs,⁴¹ as well as ultrasensitive magnetometers, could be greatly impacted by identifying the sources of magnetic defects so that the materials and fabrication processes may be engineered to mitigate flux noise.⁴² Significant effort can be devoted, for example, to removing the oxidized surface layer of silicon substrates and/or replace it with less lossy SiN_x,^2,36 but knowledge of the specific atomic source of the noise would simplify the engineering of the silicon substrate. For instance, some TLS on sapphire, another commonly used substrate, have been calculated to be associated with changes in the magnitude of magnetic moments of OH adsorbates on the surface,²⁹ and other DFT calculations also suggested adsorbed oxygen molecules may be a source of paramegnetic defects with low spin anisotropy.³¹ Displacing OH or O₂ with nonmagnetic or more magnetically stable surface passivants²⁹ and surface annealing combined with improved vacuum hygiene⁴³ are promising routes for addressing problems arising from adsorbates. It is likely, however, that additional (intrinsic) magnetic defects will continue to cause decoherence in these systems. Thus, characterizing the magnetic stability of E′ centers on silica surfaces will help assess the practicality of using silica in SQUIDS and SQs.

We study the $Eγ′$ (puckered configuration) and Si–Si dimer (⁠$Eδ′$⁠) surface defects (Fig. 1), which are the most common charged ODC-related paramagnetic defects in silica,^44,45 and also the peroxy radical, which is important for oxygen-rich conditions. We find that these defects are magnetically stable with magnetic moments of 1$μB$⁠. Thus, these defects will not change their magnitude of magnetization easily on the surface, in contrast to, e.g., OH on sapphire. This of course does not preclude the defects from being trapped at interfaces or at the surface and contributing to noise through previously proposed spin-flip (magnetization rotation) or spin-hopping mechanisms associated with coupling to nearby tunneling TLS or phonons.^28,40

We employ density functional theory (DFT)⁴⁶ with the local spin density approximation (LSDA) to study the thermodynamic and magnetic properties of the oxygen surface defects. The Vienna Ab-Initio Simulation Package (VASP) with the projector augmented-wave method is used.^47–50 The calculations are converged with respect to plane wave cutoff and k-point sampling using the α-quartz unit cell, which contains 2 SiO₂ formula units (f.u.); a 6 × 6 × 6 Monkhorst-Pack⁵¹ k-point grid centered at Γ gives energy convergence to $2.2×10−5$ eV for the unit cell. A plane wave cutoff of 700 eV gives convergence to 1 meV. The forces on all ions are relaxed to less than 0.01 eV/Å.

Positively charged states are created within the Ewald formalism by removing electrons from the system and compensating with a homogeneous negative background charge in the standard way by setting the DC (G = 0) component of the potential in reciprocal space to zero.

The magnetic stability of each defect is estimated by calculating $ΔEmag$⁠, the energy difference between its ground state and the first excited state with a different magnetization. In order to find the ground state, no constraints are placed on the magnetization. Then, starting from the ground state geometry and magnetic configuration, the difference between the populations of the spin up and spin down channels is constrained to have a different (integer) magnetization from the ground state. For example, a magnetization of 1 $μB$ indicates the difference between the number of spin up and spin down electrons is one, while the non-magnetic state has the same number of spin up and down electrons.

Unfortunately, DFT does not provide enough accuracy to predict energy differences at or below 10⁻⁶ eV, which is the thermal energy at milliKelvin temperatures (the typical operating temperature range of the superconducting devices of interest). The relative accuracy of DFT in this context is believed to be about 1 to 10 meV, which corresponds to a thermal energy of ∼11 K. Higher levels of theory, such as hybrid density functionals like HSE,⁵² can improve the estimation of the band gap in materials, but still present difficulties distinguishing precise energy differences less than 1 meV. So, in this context, we consider ΔE_mag≲10 meV to be small $∼0$⁠. The LSDA, compared to the generalized gradient approximation (GGA) or hybrid functionals, gives an overestimate (upper bound) to the energy differences computed in this study. We tested the GGA and HSE06 functionals on a few defects and find no significant qualitative differences in the results for ΔEΔ_mag, and indeed find that the LSDA yields an upper bound. Thus, we report the most conservative LSDA values. In addition, we point out that the athermal noise observed to limit coherence in measured devices suggests a higher effective temperature in the relevant device regions than the cryostat base temperature, but this effective temperature is still below the expected accuracy of DFT (effective temperatures up to ∼200 mK correspond to energies about 0.02 meV).

Realistic surface models, which are periodic in two dimensions with a vacuum region in the third, require a slab thick enough to contain a bulk-like center. We converge our calculations with respect to the slab thickness and the vacuum size in order to confirm that the defects on the surfaces are screened by both the bulk-like region in the center of the slab and the vacuum. Finite size convergence studies of defect formation energies and magnetic stabilities showed that a simulation cell with a surface area of 100–330 Ų/side and a vacuum spacing of 15 Å is adequate to model the dilute limit. The dimensions of the supercells used in our calculations are summarized in the  Appendix. The cell shape is allowed to change when relaxing slabs, but the volume is held fixed.

In SiO₂, dangling bonds can originate from oxygen vacancies or peroxy radicals, which are common in both stoichiometric and non-stoichiometric silica. In stoichiometric silica, the oxygen deficiency center (ODC) is the generic term for an oxygen vacancy. While the concentration of the E′ defects is known to be highest at the Si/SiO₂ interface in electronic devices,^5,53 the same paramagnetic defects are also known to exist at the surface of silica.^21,54 Since the structure and characteristics of dangling bonds on quartz surfaces are similar to those in bulk or surface amorphous silica, simulations of crystalline slabs are preferred to probe the magnetic stability of individual defects without convolution with the distribution of amorphous structural defects. Indeed, crystalline surfaces have been used successfully as experimental models for bulk amorphous silica.⁵⁵ Thus, we employ different cuts of the α-quartz surface with a variety of terminations as a proxy for the amorphous surface.

The propensity for amorphous materials to exhibit 1/f noise from a Lorenztian distribution of frequencies, as small changes in the local environment of a TLS change its frequency of oscillation, has been well established in the literature.²² Thus, a well-defined two level system on the surface of the crystalline α-quartz model will exhibit a spread in ΔE_mag on amorphous silica when averaged over a macroscopic area. The surface of α-quartz is an adequate model for silica in this context because the $Eγ′$ and $Eδ′$ have similar structure and characteristics on both the crystalline and amorphous surfaces. Since DFT and even higher level theories are not sensitive enough to give accurate ΔE_mag at millikelvin temperatures, using a more complicated amorphous model would not lead to more accurate calculations of energetics or any differences in the conclusions below. As a matter of nomenclature, since we are using the α-quartz surface as a proxy for the amorphous surface, we will identify the surface defects here by their amorphous designations rather than their crystalline designations. Note, the $Eγ′$ is geometrically similar to $E1′$ in bulk α-quartz.

α-quartz is hexagonal with the P3₂31 space group (#154) and we simulate the left-handed chirality. Our calculations give a lattice parameter of 4.884 Å, c/a ratio of 1.1023, and internal coordinates of $(0.4658,0.0,0.0)$ and $(0.4120,0.2740,0.1141)$ for Si and O, respectively, which match the experimental and previous computational values well.^56–58 

In contrast to defects in a crystal, in amorphous silica the local environment around each dangling bond will vary at each defect site. A variety of geometries can be simulated using different α-quartz surface terminations, although the defects that can be simulated on the (001) surfaces are limited compared to the defects on the other surfaces considered [(100), (101), and (10$1¯$⁠)], due to the unique surface reconstruction of the (001) surface. Figure 1 shows the geometries of the $Eγ′$ (puckered configuration, PC) and $Eδ′$ (dimer) defects. The peroxy radical is simulated by adding O₂ to the $Eγ′$ dangling Si bond, as discussed and shown later (Fig. 7).

The dimer configuration can be simulated with the (001) and (100) α-quartz surfaces, while the PC can be simulated on the (100) surface. The chains of two-Si rings on the (100) surface, as shown in Fig. 1 (a), are representative of the (101) and (10$1¯$⁠) surfaces, while the (001) surface has a unique honeycomb pattern of six-Si membered rings.

The magnetic stabilities of the surface E′ and oxygen defects are determined by calculating the excitation energy to change the magnetization of their ground state, $ΔEmag$⁠. For example, on the surface of sapphire, an OH adsorbate has very small $ΔEmag$ (⁠$<0.04$ eV with the LSDA) between the ground magnetic and non-magnetic excited state, due to the ease of spin density redistribution between oxygen atoms near the OH adsorbate.²⁹ We hypothesized that similarly small $ΔEmag$ could be associated with intrinsic paramagnetic oxygen defects on the silica surface, but find the energy differences to be much larger, as summarized in Table I.

In Sec. III A, we will present the details of the $ΔEmag$ calculations between ground and excited state magnetizations of the oxygen deficiency center defects in the neutral and 1+ charge states, both for the dimer and PC configurations, and also of the peroxy radical. The formation energies of these defect structures are compared to those reported in the literature in Sec. III B.

All the $VO1+$ on both the (100) and (001) surfaces have $ΔEmag$ smaller than 350 meV, including a geometry on the (100) surface without any dangling bonds (the 2III-O geometry, which is formed by two Si atoms that rebond with nearby oxygen to form three-coordinated oxygen). These small $ΔEmag$ values (but still large on the scale of milliKelvin cryogenic cooling) are reported in Table I and are the focus of the following section. The reconstructed bare surfaces (with no oxygen vacancies) as well as all neutral $VO$ are much more magnetically stable (with large $ΔEmag$⁠) than the $VO1+$ defects.⁵⁹ For example, of all the reconstructed vacancy-free surfaces, a metastable (100) surface with dangling oxygen bonds gives the smallest $ΔEmag=1.9$ eV between the non-magnetic ground state and one with a magnetic moment of 1 $μB$⁠.

The charged $VO1+$ defects in Table I are simulated by first removing an oxygen atom from representative oxygen environments on either the (001) or (100) surfaces and allowing the atoms to relax, then subsequently removing an electron to find the relaxed geometry in the charged state. All the $VO1+$ have a ground state magnetization of 1 $μB$⁠. The peroxy radical is simulated by adding two oxygen atoms onto the dangling Si bond on the (100) surface and removing an electron. It also has a ground state magnetization of ∼1 $μB$⁠.

The (001) surface is tested with a small 2x2 reconstructed surface and a larger 4x4 reconstructed surface. The charged $VO1+$ results in a surface $Eδ′$ defect. The $ΔEmag$ values reported in Table I are from the larger 4x4 surface, as the magnetization density was found to extend over the periodic boundary in the smaller (001) 2x2 cells. The $ΔEmag$ from two distinct surface oxygen vacancies [labeled 1 and 2 in Fig. 1(b)] were calculated; when these oxygen are removed they create Si–Si dimers with very similar formation energies and $ΔEmag$⁠. While sites 1 and 2 are indistinguishable in the honeycomb pattern on the surface, the site 1 surface oxygen is part of a three-membered Si ring in the subsurface [indicated by the dashed lines in Fig. 1(b)], but that of site 2 is not. In Table I only the site 1 results are listed for brevity.

Of the four sites simulated on the (100) surface, the $VO1+$ PC (⁠$Eγ′$⁠) formed from removing an O atom that bridges two 2-Si rings (labeled b in Fig. 1(a) and in Table I) has the smallest $ΔEmag$⁠. The PC formed from removing the lower oxygen in the 2-Si ring (⁠$VO1+$r) has a 58 meV larger $ΔEmag$ than the bridge site. The peroxy radical has a similar $ΔEmag$ as the PC defects.

The formation energy of neutral and charged dimer and PC oxygen deficiency centers in the bulk have been previously calculated in many first-principles studies of E′ centers. A variety of approximations are made in these simulations, such as using quartz as a model for silica or using clusters rather than extended solids in the simulation. In many cases, these approximations have been shown to be valid; for example, the similarity of the $E1′$ center in quartz and the $Eγ′$ in silica is so well established, that quartz surfaces are used to probe defects in silica even in experiments.⁵⁵ The calculations of dimers and PC on α-quartz surfaces presented in this study will add to the extensive computational literature of formation energies of E′ centers and oxygen vacancy related defects in SiO₂.

The formation energies, $ΔEform$⁠,⁶⁰ under oxygen-rich conditions are computed with the LSDA using

where $Edefect$ is the total energy of the slab containing the defect, $Eslab$ is the total energy of the stoichiometric reconstructed slab, $μO2$ is the chemical potential of molecular oxygen, q is the charge of the defect, and E_F is the Fermi energy. In our calculations, the valence band maximum, $EVBM$⁠, for the (100) slab is −1.94 eV and for the (001) slab is −2.78 eV; $μO2$ is calculated to be −10.486 eV (triplet configuration).

In Table I, we report $ΔEform$ in the limit of the chemical potential of oxygen equal to that of molecular O₂ in its standard state (⁠$μO=12μO2$⁠) and set the Fermi energy to the valence band maximum of the slab. This corresponds to conditions that would favor hole capture, but not oxygen vacancies. In oxygen rich environments, the peroxy radical is preferred. In fact, annealing the silica surface can convert E^′ centers to peroxy radicals.⁶¹ Other oxygen chemical potential conditions (including corrections for entropy) correspond to a uniform shift of all formation energies. Various external conditions give rise to shifts in the Fermi energy, which we will briefly discuss in Sec. IV.

In all our simulation cells, the neutral oxygen vacancy (NOV) creates a Si–Si dimer and is a color center called the ODC(I) in the spectroscopy community.^12,16 The ODC(I) has nearly identical optical characteristics in both silica and α-quartz,¹² making quartz a reasonable model for the ODC(I) on the silica surface. A NOV that captures a hole and retains the dimer geometry is associated with the $Eδ′$ defect in amorphous silica. The ODC(I) is understood to be a precursor for paramagnetic E′ centers in silica.^16,62

We calculate the formation energy for the diamagnetic NOV to be 5.18 eV on the (001) surface and 5.21 eV (bridge) and 5.32 eV (ring) on the (100) surface. Table I lists the calculated formation energies of the charged $VO1+$ defects, with relaxed geometries, on the (001) and (100) surfaces.

If the Fermi energy is at the center of the band gap, the NOV are thermodynamically preferred on both the (001) and (100) surfaces, while the $VO1+$ is a metastable defect. Our calculations of the formation energies of the charged and neutral ODCs are in agreement with previous calculations, showing that the neutral ODC is thermodynamically more stable.^16,17,20 However, some calculations of the formation energy of E′ defects set the Fermi Energy 2.25-3.3eV above the VBM^63,64 (much closer to the VBM than to the center of the gap). On the (100) surface, if the Fermi energy is taken to be ∼2 eV above the VBM or lower, the 1+ PC configuration is more stable than the NOV.

We find that the $Eγ′$ is more stable on the (100) surface than in the bulk of α-quartz and silica.²⁰ Figure 2 is a schematic of the equilibrium charge state diagram of the NOV and $Eγ′$ on the (100) surface. In our simulation of the surface, we find the formation energy difference between the NOV and $Eγ′$ is 2.0 eV for the ring site at the VBM, as indicated in Fig. 2. The energy differences between the NOV and E′ defects in bulk quartz and silica are also indicated schematically. The formation energy difference between the neutral dimer (NOV) and the $E1′$ in bulk α-quartz is about −0.35 eV.⁶³ The $Eγ′$ is lower in energy than the NOV in silica, but only by about 0.30 eV⁶⁴ when E_F=E_VBM. Thus, we can expect that $Eγ′$ defects are more likely to be found in higher concentrations on the surface than in the bulk, since they are 2 eV more stable.

Most previous calculations of $Eform$ are for the bulk defects and span a wide range. The formation energies of neutral defects on the surface fall within the range for bulk defects, which we summarize here. Values of the oxygen chemical potential vary in the literature, so for comparison, we shift those that differ to a common reference of molecular oxygen. The formation energy of NOV in α-quartz has been reported between 4.5 and 7.9 eV.^{16,20,64–66} The spread of energies is due to the way the formation energy is calculation²⁰ and the simulation model, in particular using clusters⁶⁶ versus supercells (and potentially poor reporting of the oxygen chemical potential and Fermi energy. Our calculated value for the NOV is 5.8 eV, which is slightly larger than the value calculated in a similar supercell simulation with 24 formula units.⁶⁴ The formation energy for the neutral PC has been reported as 7.4 eV in the bulk,²⁰ while we calculate it to be 5.9 eV on the (100) surface. The formation energy of NOV in amorphous silica has been reported as 6.1–6.8 eV for the PC¹⁶ and 4.2–5.71 eV for the dimer.^16,63,67

The formation energy of the charged defects depends on the Fermi energy. In order to compare our values to those in the literature, we shift the oxygen chemical potential to that of molecular oxygen and set the Fermi energy to $EVBM$⁠. The formation energies of positively charged defects in amorphous silica have been reported as 5.4–5.9 eV for $Eγ′$ and 5.8 eV for the $Eδ′$⁠.⁶³ Smaller supercell calculations of α-quartz reported the formation energy of the charged dimer to be 4.5–5.15 eV,^20,64 while 7.5 eV was reported for the charged PC.²⁰ The energy of hole and electron trap states from ODC have been also been calculated with first-principles in Refs. 68 and 69, but formation energies were not reported.

Formation energies of defects on hydroxylated α-quartz surfaces using embedded cluster calculations are given in Ref. 19, where $Eform$ for the PC configuration is calculated to be 6.36 eV. Using a periodic slab model, we find a lower formation energy and different electronic structure than using embedded clusters, as discussed below in Sec. IV.

While first-principles calculations seem to suggest prohibitively high formation energies of the ODCs and E′ defects in crystalline quartz and silica, these defects are well known to exist experimentally and are formed through high energy processing or irradiation. To characterize these defects in silica, high concentrations have been generated with gamma rays or high energy neutrons.^45,70 These defects also are found in thermal silica grown on α-quartz.¹ Higher concentrations of these defects are found in plasma-enhanced chemical-vapor-deposited silica.⁷¹ Additionally, high voltages across silica in electronic circuits, such as in metal-oxide-semiconductor field effect transistors (MOSFETs), can inject electrons and holes into the SiO₂, which converts diamagnetic ODCs into paramagnetic E′ centers.⁵³ In summary, positively charged defects have relative large formation energies on the surface and in the bulk of silica, but can be induced by charge trapping under various external conditions. These charged E′ centers persist long enough to be measured and affect circuit performance, indicating that there is an activation barrier to return to the lower energy neutral state. Our own calculations are consistent with the presence of such a barrier that makes the charged defects metastable, as discussed in detail in the next section.

While the energy differences between the magnetic ground state and non-magnetic excited state of $VO1+$ defects and the peroxy radical are somewhat small (less than 0.3 eV), all the defects tested are magnetically stable at milliKelvin temperatures. Even the smallest energy difference of 175 meV for the $VO1+$ defect on the (001) surface corresponds to a thermal energy of 2000 K. Other processes in the measurement system, such as microwave or infrared leakage and hot carrier effects in the devices, can contribute higher effective temperatures up to ∼1K, which is observed as athermal noise at lower temperatures.^72,73 However, fluctuations of the magnitude of magnetization of E′ centers is unlikely to be a source of the athermal noise even at these higher effective temperatures, because the $ΔEmag$ are still too large.

Due to their magnetic ground states, E′ centers and peroxy radicals may still be a source of magnetic flux fluctuations at milliKelvin temperatures via previously proposed mechanisms, which require the electron and neighboring atoms to tunnel between minima.^28,40 In order for these defects to be a source of noise, they must be at least metastable. Paramagnetic dangling bonds on III-Si are known to exists at interfaces between Si and SiO₂, stabilized in part by the low chemical potential of oxygen at the interface.²⁰ At the surface, the chemical potential of oxygen is much larger, so the concentration of ODCs and E′ centers would be lower. While we don’t make any quantitative predictions of the concentration of E′ on the surface, our results suggest that the (charged) $Eγ′$ defect on the surface is metastable and long-lived. Measurements with SQUID magnetometry of the magnetic behavior of silica surfaces have observed ferromagnetic behavior in mesoporous silica, which contains a very high surface-to-volume ratio; thus, experimental evidence suggests a high concentration of coupled paramagnetic defects on silica surfaces.⁵⁴ 

As noted above, the $Eγ′$ is more stable on the surface than in bulk α-quartz and the 1+ charge state is thermodynamically preferred, if the Fermi energy is slightly shifted from the center of the band gap towards the VBM. In fact, in the puckered geometry, the 1+ charge state is up to 3.26 eV lower in energy than the neutral charge state for these conditions. Furthermore, once the $Eγ′$ defects are formed by trapping a hole at an ODC, there is a significant barrier to convert back to neutral, leaving paramagnetic defects at the surface.

We probed this activation energy barrier to return to the neutral dimer configuration by attempting to relax the puckered configuration in the neutral charge state and also by analyzing the approximate configuration coordinate diagram as a function of Fermi energy to confirm the presence of a finite energy barrier. Starting from the puckered configuration, we added an electron back into the simulation cell (creating a neutral defect) and then allowed the atoms to relax. The puckered configuration remained after the forces on the atoms were relaxed to 0.05 eV/Åor less and the energy decreased by 1.7 eV, indicating a finite activation energy barrier to return to the dimer configuration. Thus, once a $Eγ′$ defect is formed on the surface, it will persist in that metastable charge state for significant time.

The paramagnetic electron in the E′ centers on the surface arises from a state deep in the valence band, indicated by the arrows in Figs. 3 and 4 for $Eγ′$ and $Eδ′$⁠, respectively. The red line in each of these figures plots the difference between the spin up and spin down density of state (DOS), magnified 5 times for clarity. The deep defect state indicates that the unpaired electron is more stable than the O-p band-like states that make up the top of the valence band. The majority of the occupied states above the defect band, up to the VBM, are delocalized surface states not associated with the localized ODC defect, which we determined by examining spatial projections of the partial charge density associated with each state in the DOS (not shown).

The low energy non-magnetic states for both the dimer and PC configurations arise from a low energy empty defect state in the gap (near the valence band maximum). In Figs. 3(a) and 4(a), the DOS show an empty state near the top of the valence band (within 0.3 eV of the VBM). The energies of these empty defect states above the VBM are similar in magnitude to the respective $ΔEmag$ for the defects. In the non-magnetic excited state, some charge density is transferred from the top of the valence band to the defect state, which explains the approximate correlation of these energies. The partial charge density associated with this empty state shows that the state is associated with the dangling bond, either in the $Eγ′$ defect (Fig. 5) or between the two Si atoms in the $Eδ′$ defect (Fig. 6).

A significant amount of energy would be required to remove the unpaired electron from the deep state in the valence band for both the $Eγ′$ and $Eδ′$ defects. Thus, the non-magnetic excited states do not simply differ from the magnetic ground states by shifting spin density from the filled defect state to the empty defect state, which would be a high energy process, but rather spin density is shifted from the plethora of nearby oxygen p orbitals forming the top of the valence band into the empty defect state. This is illustrated in Figs. 5 and 6, which show the differences in total charge density and magnetization density between the excited and ground magnetic states. Small values of the isosurfaces are plotted to exaggerate the surfaces for clarity.

Figure 5 shows the magnetization and charge density differences for the $Eγ′$ defect. The total charge in the dangling bond does not change: note the lack of charge density difference on the dangling bond in Fig. 5(b). In contrast, the magnetization difference between the ground and excited states is located on the dangling bond. Some spin-down density from the dangling bond is transferred into the nearby oxygen p orbitals, while at the same time some spin-up density is transferred from the nearby oxygen p orbitals into the dangling bond, giving exactly equal integrated up and down spin densities in the defect. Thus, the total charge in the dangling bond is still nearly 1e in the non-magnetic excited state, but there is no longer a magnetic moment associated with the defect.

The $Eδ′$ excited state shows some dissimilarities from the $Eγ′$ excited state. In contrast to the $Eγ′$ defect, where the total charge density in the defect stays constant, there is a slightly increased charge density in the Si–Si dimer in the non-magnetic excited state, as shown in Fig. 6(b), although we note the extent is exaggerated in the figure by plotting a small isosurface value. The magnetic moment in the Si–Si bond decreases in the excited state by removing spin-down density and increasing spin-up density, which goes into the empty defect state near the VBM. An increase in total charge density results because slightly more spin-up density is added to the Si–Si bond than spin-down density is removed. The up and down spin densities are transferred to and from, respectively, the surrounding oxygen p orbitals that make up the top of the valence band, similarly to the magnetization difference seen in the $Eγ′$ defect. This effect is highlighted in Fig. 6(a).

The E′ defects will be more likely formed in environments with low oxygen chemical potential. In environments with high oxygen chemical potential, peroxy radicals are more likely to exist on the surface, as we calculate the formation energy to be negative. Figure 7 shows the magnetization density difference and DOS for the peroxy radical on the (100) surface. The electronic structure of the peroxyl radical on top of the (100) quartz surface shows similarities with the $Eδ′$ and $Eγ′$ defects, in that the empty state is very close in energy to the top of the valence band. There are three other unpaired electrons with energy states deeper in the valence band. In total there is a magnetic moment associated with the peroxy radical defect, which can contribute to flux noise.

The mechanism of changing the magnetization of an oxygen-related defect in silica by transfer of spin density between the defect state and neighboring oxygen atoms is similar to that of the OH adsorbate on sapphire, which exhibits a very small $ΔEmag$⁠.²⁹ However, the energy scale for the change of magnetization of these intrinsic defects on the silica surface is much higher than the effective temperatures of operation of devices of interest, including SQUIDs and superconducting qubits. It is possible that E′ centers contribute to decoherence in these superconducting devices, since they have magnetic ground states and the 1+ charge states are at least metastable; experiments support this view since higher concentrations of E′ centers in plasma-enhanced chemical vapor deposited (PECVD) silica compared to thermal silica⁷¹ correlate with a higher loss tangent.^9,42 However, since the $ΔEmag$ is too large for milliKelvin temperatures, further work is required to establish the precise mechanism by which silica surfaces, interfaces, and bulk can contribute to flux noise through fluctuations of the magnitude of magnetization of defects versus rotation of paramagnetic moments. We examined intrinsic defects on the silica surface associated both with oxygen-poor (ODCs, E′ centers) and oxygen-rich (peroxy radical) conditions; while the defects are very different in atomic character, each stabilizes an unpaired electron with energy below the top of the valence band, giving a magnetic moment that can contribute to noise in both environments.

In summary, we have gained an atomistic understanding of the stability of the peroxy radical and E′ defects on α-quartz surfaces, as a proxy for silica surfaces. Our calculations of the formation energies of the peroxy radical and $Eδ′$ indicate that they are likely paramagnetic defects on the surface of silica. In particular, we find that the (charged) $Eγ′$ defect is metastable and is stabilized on the surface compared to the NOV in the bulk; when charged it may contribute to flux noise in devices.

This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract No. DE-AC52-07NA27344 with support from the Laboratory Directed Research and Development Program, tracking numbers 12-ERD-020 and 15-ERD-051.

The size and shape of the supercells used for the slab calculations are given in Table II, along with each of their computed surface energies.