Bounds to electron spin qubit variability for scalable CMOS architectures

The interface between silicon and its oxide permeates most of the technology that enabled the digital era, and as such it is one of the most studied materials in human history. The practical process engineering advantages of silicon dioxide contrast with the complex chemistry of this material, and the decades of research that has underpinned the development of the CMOS industry. As we approach a new era of quantum technology, this know-how is largely considered a major advantage for technologies such as silicon-based spin qubits Maurand2016 . In particular, the similarity between quantum dots defined by gate electrodes on top of a silicon/silicon dioxide interface and the MOSFET transistors in materials, design, and fabrication enables the integration of manufacturing techniques exclusive to semiconductor foundries onto the scaling of quantum processors Veldhorst2017 .

Here, we specifically treat the case of such qubits formed in quantum dots at the Si/SiO₂ interface, which are compatible with the high-yield integration of on-chip electronic components, and refer to these as CMOS qubits. We note, however, that other forms of silicon-based quantum dots can be manufactured, for instance by leveraging a Si/SiGe quantum well philips_universal_2022 . These materials present their own complex challenges and advantages and impose significantly different architectural choices compared to the Si/SiO₂ interface, and hence are left out of this investigation.

Despite a relatively late start Veldhorst2014 , the performance of silicon qubits has led to fidelity levels comparable with more well-established quantum technologies like superconducting or ion-trap qubits noiri_fast_2022 ; mills_two-qubit_2022 ; Yang2019 . The recent demonstration of repeatable high fidelity two-qubit operations across three nominally identical CMOS devices tanttu2023 signals the beginning of an age focused on extensive repeatability of the high performance to achieve scaled architectures. However, despite the improvements in gate uniformity demonstrated by the integration of foundry-level manufacturing techniques zwerver_qubits_2022 , new concerns are stirred by the fragility of spins to effects that are ignored in classical transistor technology, such as variations in spin-orbit coupling (impacting one-qubit frequencies), exchange interaction (two-qubit couplings) and valley splitting (nearest excitation energy).

It is only by understanding this variability quantitatively that it is possible to develop a sensible scalable quantum processor architecture. Enabling the breakthrough applications of quantum computation requires millions of qubits to perform error correction gidney_how_2021 ; beverland_assessing_2022 , and it is infeasible to address all of these qubits with pulses catering to their particular parameters with wires individually running from the room temperature controllers. Instead, this variability must be embraced and corrected through the combination of on-chip electronics operating at cryogenic temperatures Veldhorst2017 ; vandersypen_interfacing_2017 , and robust quantum control pulses that will be shared among several qubits hansen_pulse_2021 ; hansen2023entangling . Designing control cells that can offset qubit properties across a large range and with sufficient accuracy is only possible if the electric tunability compensates for the range of variations. As we will discuss in this paper, this interplay between variability and tunability differs between qubit parameters (one-qubit frequencies and two-qubit couplings, etc), so a different control strategy must be adopted in each scenario.

Qubit variability is caused by the same factors that affect conventional CMOS technology, such as strain, fabrication defects, accidental introduction of charged impurities in the oxide, interface roughness, and so on asenov_simulation_2009 . Industrial foundries focus most of their efforts in addressing these issues zwerver_qubits_2022 ; elsayed_low_2022 ; Intel2019 , with a central role played by the choice of materials for substrate, dielectrics, gate metals, etc Saraiva2022 .

In the centre of the discussion is the dielectric interface where the electrons are confinedLawrie2020 ; Burkard2023_semiconductor_qubits . High fidelity silicon qubits have been measured in Si/SiGe and Si/SiO₂ heterostructures philips_universal_2022 ; mills_two-qubit_2022 ; noiri_fast_2022 ; Yang2019 ; tanttu2023 . In the case of Si/SiGe heterostructures, a thin layer of uniaxially strained silicon binds the electron due to the conduction band shift caused by the strain when compared to the relaxed Si_xGe_1-x alloy. The Si/SiGe interface provides, in general, reduced levels of interfacial disorder compared to that of Si/SiO₂ Lawrie2020 . Potential shortcomings of Si/SiGe technology are the reduced gate control when compared to MOS devices and the limited tolerance of the material stack to high temperature annealing processes, commonly adopted in the CMOS industry. More information on this material and its comparison to oxides can be found in Ref. Saraiva2022 . In the context of spin qubit variability, comparing an oxide interface to Si/SiGe is hard because the nature of disorder in the two materials is different - alloy disorder and miscut angles in SiGe, compared to amorphous oxidation in SiMOS. Moreover, the dominant effect of spin-orbit coupling being studied here is masked by the presence of micromagnets which are commonly adopted in SiGe qubit architectures philips_universal_2022 .

In this paper, we focus on qubits at the Si/SiO₂ interface. SiO₂ does not have a regular lattice structure when thermally grown on the silicon surface, so the interface is atomically rough. The higher levels of interfacial disorder when compared to Si/SiGe are attributed to this roughness and to the presence of fixed charge defects that can be either at the Si/SiO₂ interface, in the bulk of SiO₂ or at the metal-oxide interface Lawrie2020 ; elsayed_low_2022 ; Intel2019 ; Müller_2019_traps_solids . Potential advantages are in the higher electrical tunability and compatibility with conventional CMOS technology, which benefits the integration with on-chip electronics Veldhorst2017 .

The roughness of the Si/SiO₂ interface is one of the most critical sources of disorder for CMOS spin qubits jock_siliconMOS_sandia_2018 ; Veldhorst2015_spin_orbit ; Tanttu2019 ; Ruskov2018 ; Ferdous2018 ; Gamble2016 . A more recent paper indicated that a second source - charged impurities in the oxide - could dominate over interface roughness martinez_variability_2022 , at least in the case where a micromagnet is integrated to allow electron spin qubits to be driven electrically. Electric driving requires a large spin-orbit coupling, which exposes the spins to the impact of charge impurities. In this paper we focus on spin qubits driven magnetically Veldhorst2014 ; vahapoglu_coherent_2022 , which do not have this requirement. These qubits can be controlled coherently without the inclusion of micro-magnets, thus preserving the low spin-orbit coupling of electrons in silicon and protecting the spin from electric fluctuations. We will show in this paper, that under these conditions the remaining spin-orbit variability is interface-induced with charge impurities only affecting the qubit by shifting the quantum dot formation against the roughness profile of the interface.

Recent advances in CMOS quantum dot fabrication allowed for sufficient yield to create a number of small-scale quantum processing units and measure their variability. This work combines measurements of 12 qubits across 6 different CMOS quantum dot devices, transmission electron microscopy (TEM) images of cross-sectional cuts of 6 other quantum dot devices and theoretical analysis of quantum properties of electrons in simulated quantum arrays of 49 quantum dots (Fig. 1a). Such a geometry allows us to study the impact of the self-affine scaling of the Si/SiO₂ roughness on qubit properties at different length-scales Yoshinobu1995 (Fig. 1b). Starting from a detailed view to the device architecture, quantum dot formation, and materials interface down to the atomic scale, we predict the bounds to spin qubit variations. This prediction is compared with data from qubit devices , some of which have led to manuscripts and publications (see Table 1). All devices were fabricated with geometrically identical designs (see structure depicted in Fig. 1c) and differ only in material stack compositions and spin control methods (Table 1).

Si/SiO₂ roughness The Si/SiO₂ interface has an intrinsic fractal structure Yoshinobu1995 , that has not been considered in previous variability studies (Ferdous2018 ; Gamble2016 ; martinez_variability_2022 ) . Our decision to simulate a 7$\times$7 dot array is motivated by understanding how this fractal scaling of the interface roughness impacts qubit properties at different length scales (see Fig. 1b). One-qubit properties, for instance, depend on the roughness obtained within the quantum dot diameter ($\sim$ 15 nm), while two-qubit properties are related to the interdot distance (30-60 nm). The roughness amplitude can differ significantly between these two scales due to this fractal structure, so a quantitative characterization of the interface is necessary.

Our model of the Si/SiO₂ (Fig. 1.b) is based on the roughness observed in TEM images (Fig. 1d-e) Goodnick1985 ; Yoshinobu1995 , allowing us to include more realistic features. By convolving these images with the expected face-centred cubic lattice of monocrystalline silicon we can mathematically discern the interface and analyse the roughness at different scales Goodnick1985 , quantified through its power spectral density (PSD) as a function of the in-plane correlation length scale $\lambda$ Yoshinobu1995 . Our work incorporates a theoretical study of the realistic scale-dependent fractal structure of the roughness Yoshinobu1995 and the development of theoretical tools capable of capturing this multiscale physics in a realistic device model (Fig. 1f-g).

Combining multiple TEM images, we obtained in Fig. 1h a consistent roughness pattern characteristic of a fractal scaling down to the silicon lattice parameter of the form $\text{PSD}^{\text{1D}}(\lambda)\propto\left(\frac{2\pi}{\lambda}\right)^{-1-2H}$. We estimate a Hurst exponent of $H$=0.28 Jacobs2017a (details in Supplementary Fig. 6 and methods section). The root-mean-square roughness also scales up with the lateral region $\lambda$ as $\text{RMS}\left(\lambda\right)\propto\left(\frac{2\pi}{\lambda}\right)^{-H}$. As shown in Fig. 1i this roughness pattern is consistent across all devices measured, and extends up to half a micrometre . We note that these levels of roughness are typical for industry-standard interfaces Intel2019 . Our computer-generated interface in Figs. 1b,f and g was also designed to mimic these features.

A direct conclusion from Fig. 1h-i is that the size of the dots (approximately 15 nm for all devices studied) and the separation between dots (approximately 50 nm) will have a large impact on the qubit exposure to surface roughness, and that the interface distortions within a quantum dot are smaller than those between neighbouring dots. The root-mean-square (RMS) is approximately 0.15 nm within a dot (approximately 1 monolayer of the Si lattice), and almost twice for the double dot (RMS is 0.3nm, approximately 2 monolayers of the Si lattice).

Variability of quantum dot structure and excitation energy

The consistency of this roughness pattern across devices allows to theoretically forecast its impact on qubit performance at scale. We simulate the spin qubit variability of quantum dots formed in different subsections of the computer-generated interface in Fig. 2a-b. Unless indicated differently, all quantum dots are simulated using the same electrostatic potential, which is based on COMSOL simulations of realistic digital models of our devices (see Fig. 1c-g and methods section). The only difference between simulations is the location of the quantum dots in the surface profile. Such a model allows us to incorporate the self-affine characteristic of the Si/SiO₂ interface in our simulations.

In this paper, the quantum dots are simulated in a 7$\times$7 grid array. Other architectures are also being explored SpiderWeb2022 , in part, because the practical design of such a dense array requires a sophisticated fabrication process with multiple metal layers to route the signals to the gates Veldhorst2017 . While dense wiring in multiple metal layers is routinely integrated in front-end-of-line industrial processes, qubit demonstrations using this integration have only recently been explored HRL_Sledge_2022 . Moreover, this dense array leaves no space for interspersed readout devices such as single-electron transistors, and would be dependent on a gate-based readout approach Gonzales_Readout_2021 or would require quantum information to be shuffled to the edges of the array for readout Peta_swap_sigillito_coherent_2019 . Our results for qubit variability are not drastically affected by the choice of a grid array, except for a small degree of nearest neighbor correlations (see Supplementary Fig. 11), which once simulated across the full 450 $\times$ 450 nm Si/SiO₂ computer-generated interface, provides us with sufficient sampling to obtain statistical analyses accurately.

To understand how this roughness affects the quantum behaviour of electrons, it is necessary to focus on their wavefunction at the atomic scale (Fig. 2b). We use an atomistic tight binding model of Si and SiO₂ which incorporates relativistic effects and the impact of a magnetic field, yielding eigenstates with realistic spin and valley structure. It can be used to calculate the ground state wavefunction and a few excited states.

Despite SiO₂ not having a regular lattice structure, we simulate it by assuming an atomistically ordered virtual crystal approximation (see Methods). The material is endowed with the same lattice structure as Si and the tight-binding parameters are set to emulate the electronic structure of the interface kim_full_2011 ; Bersch2008 . This approximation allows us to simulate interface disorder atomistically as seen in Fig. 2a. In addition, we developed techniques to extract this structure and calculate properties of the disordered quantum dot that would not be obtained with a purely spectral analysis, such as the valley phase and the spin g-tensor (see Methods).

We find that the geometry explored here always leads to the successful formation of quantum dots, regardless of the local roughness profile. This is consistent with the yield of measurable quantum dots – all devices with functional gate electrodes (as determined by their influence on the charge sensing single electron transistor) could form controllable pairs of dots. Roughness mostly alters the quantum dot effective shape and centre position (Fig. 2c), with location of the electron departing from the potential minimum by less than 5 nm, with a standard deviation of 1.4 nm (Fig. 2d). The dot position in the geometry studied here is highly tunable by biasing lateral gates (approximately 5 nm/V, see Supplementary Fig. 7 and Supplementary Table 3), so that this disturbance can be corrected.

The excited orbital states are more impacted by the interface roughness. We are particularly interested in the first excited state, which corresponds to a valley excitation for a [001] Si/SiO₂ interface. The conduction band valleys along $\pm z$ crystal directions are energetically favourable due to the effective mass anisotropy, and the degeneracy between these two valleys is lifted by the sharp interface. The performance of spin qubits is strongly impacted by the interface-induced valley coupling, which creates a superposition between the two valley quantum states. This superposition creates an oscillatory behaviour at the atomic scale, which can be seen in simulations in Fig. 2e. These oscillations are known to cause variability in valley structure even for interfaces with low levels of disorder. We refer to the relative phase between valleys in this superposition as valley phase and the energy separation between the two states as valley splitting.

To have pure spin systems, valley splittings exceeding the Zeeman energy are desirable. In Fig. 2f, we show how the valley splitting can be controlled tuning the vertical electric field $E_{z}$, comparing measurements in two devices and the results of the simulations. The surface roughness causes variability in valley splitting of over one order of magnitude for a fixed electric field. The full range of valley splittings spreads from tens of $\mu$eV to a few meV, compatible with observed experimental values in CMOS devices Yang2013 ; Gamble2016 .

When comparing these valley splittings to the spin splitting, we may ignore the variability in Zeeman energy (which is only of a few parts per thousand). Therefore, if we set a relatively high electric confinement ($\approx$28 meV ${\rm nm}^{-1}$ is sufficient in our simulation) and we tune the magnetic field low enough (¡700 mT in our study), all the 49 qubits in the simulation will obey the condition of valley splitting larger than the Zeeman splitting (See inset Fig. 2f). However, the fitted distributions in Fig. 2f show that the valley splitting can sometimes be very small, even at high electric confinements. This is an exceptional event, and for low magnetic fields none of the 49 dots simulated here had valley splittings too small. However, it could potentially affect the development of large-scale quantum processors with millions of qubits as it is expected that a number of quantum dots will have the valley splitting clashing with the spin splitting. A possible solution would involve changing the number of electrons in the dot Leon2020 . In worst case scenarios, the dot must be discarded from the processor at the firmware level. A thorough discussion of the impact of this decision on error correction is presented in Ref. Nagayama_2017 .

Even with a consistently high valley splitting, qubit performance can still be impacted by variations in valley phases between neighbouring dots. The electron density will present Bloch oscillations in the $z$ direction, which are barely visible in Fig. 1f and were enhanced in Fig. 2e by taking the difference between the electron densities of both valley states (Methods). Notice that the $z$ oscillations have the same phase across the whole dot instead of conforming to the roughness of the oxide. This implies that the valley phase is well defined even in the presence of surface disorder. This phase has an impact on operations that involve two dots, such as electron tunnelling and exchange coupling, because it determines whether these valley oscillations interfere constructively or destructively tariq_impact_2022 ; Tagliaferri2018 . Fig. 2g shows the valley phases across the 49 simulated dots, revealing that dots with larger valley splittings (typically above 300 $\mu$eV) tend to have similar valley phases, near zero in our definition. We speculate that this behaviour is a consequence of most dots with high valley splittings being formed in a preferred atomic layer, with similar z (see inset Fig. 2f), such that they have aligned valley phases.

Qubit frequency variations Spin-orbit effects lead to variability in qubit frequencies of the order of $\sim$100 peV, which appears in the form of a variable $g$ factor for the spin subject to an external magnetic field. This variation represents less than 1% of the qubit frequency. The mean value of the ratio of Zeeman frequency and external magnetic field $f_{\rm Zeeman}/B_{0}=g\mu_{\rm B}$ is 27.9 GHz T^-1, with the differences occurring only at the order of tens of MHz T^-1. This is a particularity of silicon electrons, whose spin-orbit coupling is among the smallest in all quantum dot technologies. This provides CMOS spin qubits with special protection against disorder and electric fluctuations. However, these very small g-factor variations are still important for qubit operations aimed at higher than 99.9% fidelity, and they are directly linked to the roughness of the interface at the atomic scale.

Interface-induced spin-orbit coupling has two flavours in Si/SiO₂ interfaces – Rashba ($\alpha$) and Dresselhaus ($\beta$) Ruskov2018 . These two can be experimentally differentiated by measuring the g-factor dependence on the in-plane magnetic field orientation $\varphi$ Tanttu2019 . The dependence is sinusoidal with the form $g\left(\phi\right)\approx g_{0}+\alpha+\beta\sin{\left(2\phi\right)}$, where we take $g_{0}$ to be the theoretical bulk g-factor $g_{0}=1.9935$ calculated from atomistic simulations including relativistic spin-orbit effects. The difference between the frequencies of any two qubits (Fig 3.a), as well as the electric field dependence $dg/dV$(Fig 3.b have the same behaviour. All 12 qubits measured in devices A to E show behaviours consistent with this description. Qubits in the same quantum dot but with different electron numbers are considered as different in the total count, as they have substantially different g-factors (see data from device A in Fig 3.a-b). This is most likely due to different exposition to the atomic profile of the interface for electrons in different valley states. In most cases, Dresselhaus dominates both the total spin-orbit effect and its variability – with the exception of the configuration (1,1) in device A (Fig. 3.a). The Rashba coefficient $\alpha$ is typically one order of magnitude smaller than $\beta$ (Fig. 3.c-d).

We explore theoretically this variability by extracting the g-tensor of electrons in disordered quantum dots from the eigenfunctions calculated by tight binding. The results of simulations, shown as solid lines in Figs. 3a and b, are then used to extract the dependence of $\alpha$ and $\beta$ on the vertical electric field (Figs. 3c and d). The Dresselhaus effect emerges from breaking the lattice inversion symmetry near an interface, which explains why it is strongly dependent on the electric field that confines the electron against the oxide (Fig. 3c). The interface-induced Rashba effect is also dependent on the electric field, but at a smaller scale (Fig. 3d ). Notice that a couple of simulations escape the overall trend. These are valley-spin degeneracies, occurring in these simulations because some valley splittings clash with the Zeeman splitting (Fig. 2.f ) at the input magnetic field of 1 T. While these degeneracies can be used for fast electrical driving Bourdet2018 ; Corna2018 ; gilbert_-demand_2023 , they significantly deviate from the target parameters for pure spin qubits. In practice, it is possible to either tune the valley splitting out of this regime or reduce the magnetic field to reduce the probability of accidental degeneracies, as discussed in the previous section.

We can also observe in Fig. 3c that the Dresselhaus parameter $\beta$ is bounded between two extreme values for each electric and magnetic field. This can be understood by analysing the perfectly flat [001] interface model, which introduces a distinction between the two sublattices in the diamond structure Ferdous2018 . Fig. 3e shows that in this case the Dresselhaus parameter $\beta$ is maximally positive for one sublattice termination and inverts for the other – the reason for this inversion can be understood viewing the [001] terminations in Fig. 3f. In comparison, a rough interface will contain terminations in both sublattices, and the value of $\beta$ will then lie between these two limits (see also Supplementary Fig. 8).

Strategies for qubit control need to be designed according to these results in order to tolerate the natural statistical dispersion in qubit frequencies introduced by the oxide interface. Individual addressability is a particular challenge. The most common strategy explored so far for addressing a specific spin qubit relies on exploiting this g-factor variability, and driving a variable microwave field in resonance with its Larmor resonance frequency in order to induce spin rotations Veldhorst2014 . However, the fact that the spin-orbit effect has a maximum natural spread results in frequency crowding at large qubit numbers Fig. 3g, making it hard to address a given qubit without impacting other qubits with similar frequencies.

Instead, a more scalable pathway relies on a global microwave field acting on all qubits simultaneously hansen2023entangling ; hansen_pulse_2021 . Qubits will be driven in resonance with the global field (forming a dressed qubit seedhouse_quantum_2021 ), thus requiring that all qubits have identical frequencies aligned with the resonator frequency. Individual addressability is then achieved electrically, by tuning the qubits in and out of resonance via Stark shift. However, in the case of a magnetic field along [110], the range of electric control of the g-factors is insufficient to tune them into the same frequency as can be seen in Fig. 3g . The variability in Larmor frequencies can be reduced by pointing the magnetic field along [100] (Fig. 3h) and by reducing the field magnitude to less than $\approx$100 mT zhao_single-spin_2019 . These modifications lead to a standard deviation of less than 200kHz, which is ideal for driving degenerate qubits with high fidelity (assuming a 2MHz Rabi frequency in Ref. hansen_pulse_2021 ). The Stark shift tunability also decreases, such that control strategies might need to cope with the inability to tune qubits completely out of resonance hansen_pulse_2021 .

The role of charge impurities in the presence of interface disorder Semiconductor devices are in general exposed to the presence of charge impurities in the oxide layer, originated from dangling bounds, Pb-centres, oxide vacancies, and other defects Müller_2019_traps_solids ; elsayed_low_2022 . Some of these are fixed charged traps in the oxide, while others fluctuate over time, which is the origin of 1/f noise in semiconductor devices. The most concerning charges from a variability perspective are charged traps. While two-level fluctuators are still important to understand qubit noise and decoherence, their absolute impact on the variability of qubit parameters occurs at a much smaller scale stuyck_real-time_2023 ; zhao_single-spin_2019 ; ModelingShehata2023 .

Previous works have modeled charged traps as negative electron charges $e^{-}$ directly at the Si/SiO₂ interface martinez_variability_2022 . However, factors such as the positive correlation between the SiO₂ thickness with charge mobility in Hall bar devices yagi_oxide_thickness_2008 ; Intel2019 , and measured large Dingle ratio elsayed_low_2022 are possible signs that these charges might be dominantly located closer to the metal-oxide interface in some cases.

We simulated a uniform charge distribution distributed across bulk SiO₂ with density $-4\times 10^{10}$ nm^-1 Intel2019 (see Fig. 4a). Each trap is assumed to lead to a deformation of the potential of

where the trap at $\mathbf{r}_{t}$ and $\epsilon_{\text{Si/SiO}_{2}}=0.5(\epsilon_{Si\rm}+\epsilon_{\text{SiO}_{2}})=7.5\epsilon_{0}$. Here we investigate the impact of these charges on spin qubits in the presence of interface disorder (Fig. 1b).

In most cases in Fig. 4b there are no charges trapped inside the dot region, so the quantum dot wavefunction is similar to the simulations without charged impurities (see Fig. 2c). The shifts in the potential profile induce small displacements of the quantum dots below the rough Si/SiO₂ interface. This leads to variations in valley spittings and g-factors with respect to initial simulations (Fig. 4c-d). An average increase in the valley splitting is observed, which is presumably due an overall enhancement of the lateral confinement of the dots in the presence of the negative traps. This increase is typically larger for quantum dots that already had a large valley splitting (See Supplementary Fig. 12) due to their higher susceptibility to electric fluctuations. The dispersion of the dots centers is slightly larger (see Fig. 4e-f), but still small in comparison with the quantum dot size.

Only when the charge is inside the dot region, the quantum dot wavefunction is significantly affected (see, for instance, dot 2 in Fig. 4b). While the g-factor of this dot remains within range (Fig. 4c-d), the valley splitting enhancement is larger than what it is expected for typical values (See Supplementary Fig. 12). The probability of this effect happening is 5.6% for this particular charge density, as estimated from a Poisson distribution Tong_2020_variabilityCharges . This situation could also affect large scale architectures. However, it is unclear whether a qubit exposed to this type of defect would necessarily be unusable – the comparison between open and closed symbols in Fig. 4c-d indicates that the presence of traps does not increase the statistical dispersion of one-qubit parameters significantly. We will show next that the impact of traps on the two-qubit exchange coupling is also manageable with voltage offsets.

Variability of exchange interactions

Besides the single qubit gate control, roughness and charge impurities also limit the homogeneity of the exchange coupling between neighbouring quantum dots. The differences in valley phase between dots, addressed in Fig. 2g, are the first source of variability in exchange coupling tariq_impact_2022 . In the worst case, the valley phases would be random and the resulting exchange coupling would be impacted by the destructive interference of the valley oscillations. The probability of a completely destructive interference is, however, negligibly small. The typical valley interference causes at worst an offset in the exchange coupling of one or two orders of magnitude, which is easily corrected by an offset in the exchange control gate voltage given that the tunability ranges from 6 to 10 decades per Volt. We find in simulations, however, that the disorder in quantum dot position has a stronger impact.

The fact that exchange coupling is a contact interaction means that any effect impacting how the wavefunction tails off from one dot into its neighbouring dot, such as interdot distance and potential barrier height, has an exponential impact Li2010 ; Saraiva2007 . In the devices investigated here, exchange is controlled by the action of the interstitial J gate (see Fig. 5.a). This gate induces a lateral displacement of the two dots toward each other reducing the interdot distance at a rate of approximately $10$ nm V^-1 (see Supplementary Table 3). This is contrary to the picture frequently evoked in this scenario which assumes that the J gate controls the electron penetration length into the classically forbidden region between dots without affecting much its position.

In Fig. 5b we can see a method of extracting the exchange coupling by measuring the qubit resonant frequency for a randomly initialised pair of spins as a function of the J gate voltage.

To simulate this system, we use a path integral Monte Carlo approach Ceperley1995 . This method is relatively fast and intrinsically includes the effect of interactions in the electron dynamics. For each exchange simulation, we sample realizations of likely paths of two electrons in a three-dimensional double quantum dot potential obtained from a finite elements simulation (see Fig. 5c). Interface roughness can be readily included by defining a 3.1 eV step potential barrier to simulate the conduction band offset between silicon and the SiO₂ layer. The paths of these electrons are allowed to exchange between the dots, and the impact on the path action is used to estimate the exchange coupling Pedersen2010 . The full method is described in Ref. CifuentesPIMC2023 and simulation details are included in the methods section.

In Fig. 5d we compare the J-gate tuning of the exchange coupling with random surface and trap realizations against the experimental data obtained from devices A and C. We plot the exchange $J$ against $V_{J}-0.5(V_{\rm P1}+V_{\rm P2})$, which is an approximate measure of the voltage bias between the J-gate and the plunger gates (our dots are fully isolated from the reference voltage set by the reservoir). All devices were tuned to the symmetric operation point reed_reduced_2016 .

Our dataset combines qubits at the outer shell of quantum dots with one electron and 3 electrons. The exchange interaction is larger for higher electron numbers, due to the increasing overlap between the wavefunction of the quantum dots. This is observed in device A, where the exchange baseline of the (3,1) configuration is higher than in (1,1). These situations are so far not considered in the simulations, which are performed in the (1,1) electron configuration. The gate stack also has a significant role in the exchange control, affecting the effectiveness of the electric fields generated by the J gate in the channel. Here we compare devices A and C, which are both made with Pd/Ti gates with ALD oxides (see Table 1 and device architecture in Fig. 1f). Exchange control was also measured in devices E and F with Al gates in Ref. tanttu2023 , observing larger control rates and variability due to the absence of the ALD oxide and irregularities in the gate structure Saraiva2022 .

We found good agreement between exchange simulations and experimental data from devices A and C (Fig. 5d). An important part of this agreement was the inclusion of negative charged impurities in the model, together with interface roughness (see Fig. 4a). The charges that most significantly impact the exchange coupling between neighbouring dots are the ones located in the interdot channel CifuentesPIMC2023 . A single trap can reduce the exchange baseline a few decades. In total, this adds up to an average decrease of the exchange baseline of 1.5 decades in comparison with simulations without charge impurities, that have a higher baseline than experiments, as seen in Supplementary Fig. 10.

Importantly, the variability in exchange coupling can be compensated with more tunability. As seen in (Fig. 5e), the J-gate tunes the double dot potential inducing a displacement of both quantum dots to the centre by almost $5$ nm each. We can observe in Fig. 5f that this exchange dependence on displacement is consistent over multiple realizations, even when it is perturbed by the variability of the dot centres caused by surface disorder and charged traps. Because of the strong correlation between the inter-dot distance and the exchange coupling, we can associate the three orders of magnitude of exchange variability to these shifts in the interdot distance. Besides, the J-gate controls the interdot distance at an approximate rate of 10 nm per Volt, resulting in tunability ranges from 5 to 8 decades per volt in simulations and experiments (Fig. 5 g). This is large enough to compensate for the disorder and consistently hit a target “on” exchange rate, as indicated by the yellow region in Fig. 5d.

The design of scalable quantum processor cells must be guided by a precise understanding of qubit variability. Besides demonstrating a complete strategy for diagnosing qubit variations for a given choice of qubit design, fabrication process and materials, this study leads to some general conclusions about the general physics of spins under Si/SiO₂ interfaces. The main conclusion is that most of the qubit variability in current devices is explained by the roughness of the Si/SiO₂ interface roughness. The presence of charge impurities is also significant in regard to two-qubit properties. Other effects (such as strain inhomogeneity and geometrical deformation of the gates) can be mitigated down to levels that are, at most, comparable with these intrinsic mechanisms.

Another conclusion is that electric tuning of qubit frequencies (using spin-orbit effect) and exchange coupling (using barrier gates) both rely to a large extent on moving the quantum dot laterally, dragging it against the rough interface. This may lead to considerations in future designs of quantum dots and the methods for characterising the interface .

Charge noise coming from two-level fluctuators (TLFs) can also couple to the qubits through a similar mechanism. The induced displacements of the dots in the rough interface lead to small g-factor variations that can affect important qubit metrics, such as the phase coherence $T_{2}^{*}$ stuyck_real-time_2023 ; Tanttu2019 . The charge noise couples directly to the qubit Stark shift, so the methods of this paper can be applied to investigate the microscopic nature of this effect cifuentes2023Crosstalk .

One remaining question is how much improvement can be realistically expected in the interface quality, and how it impacts the qubit performance. We address the last question in Supplementary Fig. 9, showing that the main benefits would be an enhancement of the average valley splitting and a smaller exchange variability. The spin-orbit coupling is not significantly affected due to its intrinsic atomistic dependence (Fig. 3f). Methods to improve the interface quality would involve replicating previously observed correlations between roughness amplitude and different fabrication parameters, such as growth time, oxide thickness, etc fang_evolution_1997 . The characterization methods developed in this study can help assess if the improvement in roughness amplitude occurs at the length scales that are more relevant for CMOS quantum processors.

This work focused on the case of Si/SiO₂ interfaces, which has been studied enough that we are able to draw firm conclusions on its impact on qubits. Other oxides or dielectrics might also be understood adapting the methods developed here, providing pathways to shortcut the qualification of material stacks for quantum processor fabrication with the assistance of theoretical calculations.

Finally, this study realistically sets the ultimate variability of spin qubit parameters in CMOS devices. We may extract, for instance, the voltage offset that would be necessary to bring a qubit parameter to approximately the same value for all qubits in the architecture, which we call the voltage offset deviation $\text{VOD}(x)$ for each parameter $x$ (see Methods section). For example, the typical voltage offset to bring valley splittings to the same range is 0.58 V, while doing the same for g-factors requires 0.23 V if the magnetic field is pointed towards [100] and 9.1 V for a field along [110]. The smallest value is for the exchange coupling variability which can be corrected with only 0.09 V voltage offsets. That clearly reveals that some parameters can be electrically tuned to a target system-wide value while others will require the implementation of strategies to circumvent the variability with a combination of locally generated control signals and robust global control strategies hansen_pulse_2021 . These results outline the minimum demands for a CMOS spin qubit architecture.

Fabrication process The SiO₂ gate oxide (7.5-8.0 nm) was thermally grown on the silicon surface in a custom built high quality oxide furnace as part of a standard MOS device fabrication process. The gate fabrication process was iterated multiple times to improve yield, which was an enabling feature for this study, leading to the successful formation of several devices with nominally identical layouts.

Spin spectroscopy To measure the Stark shift and the difference in Zeeman splittings between the spins, we have to be able to measure the Larmor frequency of the two qubits at a given operation point. In these experiments we have double quantum dots which we use as two electron spin qubits. To begin, we initialise both electrons in the same dot forming a singlet state. We then separate the electrons and, depending on the rate of the separation, they will either end up in a $T^{-}$ state or having an odd parity, for instance an up-down state. To find the Larmor frequency we apply a fixed pulse or adiabatic microwave pulse with an antenna or resonator LauchtAdiabatic2014 . This pulse will flip a spin only if the applied frequency corresponds to the Larmor frequency. If the spin is flipped, the parity of the two spins will change too. We then bring the two electrons to a position where they are allowed to tunnel to the same dot only if their total parity is odd. If the parity is even, the electrons stay in their respective dots. We call this the Pauli-spin blockade-based parity readout Yang2020 ; Seedhouse2021PSB . The resulting charge state is read using a nearby charge sensor. Here we used a single electron transistor. The frequencies where we measured a flip in the parity of the two-spin states will correspond to the Zeeman splitting of one of the two qubits.

Measurements of valley splitting A tunnel rate-based spectroscopy method is used as an approximate measure of the valley splitting of quantum dots. The technique is as follows; 1) A repeated square-wave is applied to a dot gate, centred around a dot-to-reservoir charge transition. 2) The frequency of the square-wave is set equal to or faster than the dot-to-reservoir tunnel rate. In this mode, an electron will move in and out of the quantum dot in sync with the square wave, but for only a proportion of the square-wave repetitions due to the tunnel rate. 3) The amplitude of the square wave is swept from zero to 20 mV. The increase in amplitude causes excited state transitions to be accessed from the reservoir, which typically have increased tunnel rates to the quantum dot. 4) The changes in tunnel rate are fitted, and the splitting in amplitude between the ground state transition and excited state transition is multiplied by the dot gate lever-arm to retrieve the excited state energy. Here we assume the first excited state to be the valley excited state. This measurement technique is also explained in detail in yang2012 .

Image analysis to filter the Si/SiO₂ interface from TEM images

In the TEM in Fig. 13a it is possible to discern the periodic pattern caused by the electron diffraction on the Si crystal. We apply a set of periodic convolutions to create a different background between Si and SiO₂ regions. These convolutions are defined by four matrices $M^{cc},M^{sc},M^{cs},M^{ss}$ with dimensions $n_{x}\times n_{y}$ and entries

where $i\in\{1\ldots n_{x}\}$ and $j\in\{1\ldots n_{y}\}$. The result in Fig. 13b is obtained after taking the average of the images generated by four convolutions $M^{cc}*\mathcal{I},M^{sc}*\mathcal{I},M^{cs}*\mathcal{I},M^{ss}*\mathcal{I}$ , where $\mathcal{I}$ is the original image. Here the parameters $n_{x}$ and $n_{y}$ indicate the periodicity of the convolution in pixels in the $x$ and $y$ axis respectively. We define a parameter $N$ to be equal to $n_{x}$, and $n_{y}=0.8N$. This parameter defines the lower bound scale at which the extracted power spectral density (PSD) of the Si/SiO2 interface is still not heavily affected by the convolution.

After the convolution, the lower part of the image corresponding to the silicon region should be darker than the upper side. Then we define a threshold parameter $T\in(0,1)$ to segment both regions (Fig. 13c).We tune the parameters $N$ and $T$ until the fit clearly corresponds to the Si/SiO₂ interface in the original TEM (Fig. 13d).

Si/SiO₂ characterization of interface from TEM images We characterize the Si/SiO₂ roughness from TEM images of devices T1 to T6 (see Supplementary Fig. 6) using a similar procedure to Goodnick1985 . To filter the interface, we apply image convolutions that enhance the differences between both materials. We then perform a power spectral density (PSD) decomposition of the filtered interface

and confirm that the scaling of the Si/SiO₂ roughness is characteristic of a fractal self-affine interface  Jacobs2017a ; Yoshinobu1995 .

We obtained an average Hurst exponent of $H=0.28\pm 0.2$ and a roughness amplitude parameter of $\mathcal{C}_{0}\approx 1.4$nm³ for our devices. The PSD profiles can change for TEMs taken at different regions of the same device. This is normal as some line versions of a 2D surface can be smoother than others, and this behaviour is observed even in a surface generated numerically. We found that comparing average 1D parameters from multiple TEM images accounts for a better estimate of the global 2D profile (Supplementary Fig. 6).

RMS characterization of Si/SiO₂ interface from TEM images To compute the scaling of the RMS$(\lambda)$ over a certain interface length characterised by the wavelength of its Fourier decomposition $\lambda$, we took separated segments of the fitted line surface from TEM image of width $\lambda$ and computed the average RMS for each device. If a TEM has a width of $L=$40 nm, we would compute the RMS$(\lambda)$ for $\lambda={L/2,L/4,\ldots}$ until reaching the atomic scale $a_{0}$ Yoshinobu1995 . We can select more subsegments when $\lambda$ is small, which explains the smaller uncertainty of the average RMS for small $\lambda$. We found that the RMS scales with an exponent of $-H$ as

which is expected for the PSD profile in equation (6) Jacobs2017a .

The RMS roughness amplitude of the computer generated surface was obtained by averaging the results over line segments of the 2D surface profile, as in Supplementary Fig. 6c with the power spectral density.

Random surface generation The computer model of the surface in Figure 1b was generated with a Fourier-filtering algorithm implemented in Matlab MonaMahboobKanafi2021 , that takes as input the 1D PSD profile in equation (6) and outputs a random rough surface with similar spectral density. The output surfaces look visually similar to the ones in the TEMs. To calibrate the model we compare the average 1D PSD and 1D RMS scaling profiles, with the edge profile in the TEM images until we find a good match (Supplementary Fig. 6).

Modeling of digital twin of the devices We create a 3D structure of the device in Matlab with the software provided by the DFX library. We begin by using a physical quantum dot electron-beam lithography (EBL) design layout as the primary framework and construct it as a 3D model that closely resembles its appearance in SEM and TEM images. Our software takes into account the levels of oxidation and thermal expansion that occur during the fabrication process to finely imitate the device geometry. We tweak these variables until the 3D model looks similar to transversal and in-parallel views of SEM and TEM images.

Atomistic simulations with interface roughness We perform tight binding simulations in NEMO3D in the $sp^{3}d^{5}s^{*}$ 20-band model for Si, which intrinsically includes spin-orbit-interactions Klimeck2007 . In all simulations the magnetic field magnitude is set to $1$T, and we only vary the magnetic field orientation. To include surface disorder, we terminate the silicon lattice with the local section of the rough surface in Figure 1b. We then label all the lattice sites above (below) the surface as SiO₂ (Si), as observed in Figure 2a. The SiO₂ region is modeled with a sp³ tight-binding (TB) model under a virtual crystal approximation (VCA). The SiO₂ TB parameters are optimized to reproduce the electrical properties of the oxide, namely the bandgap of 8.9 eV, conduction band offset of 3.15 eV relative to silicon, and conduction effective mass of 0.44m₀. The VCA model in TB assumes a well defined crystal structure, in this case zincblende with a lattice constant having the same value of Si, but treats each atom as a fictitious SiO₂ atom. This is a standard way to model alloyed (SiGe) or disordered (SiO₂) materials under the VCA in the atomistic TB technique. At the interface region whether an atom is marked as a Si atom or SiO₂ atom then creates the atomistic disorder profile. The details of the model with parameter values can be found in kim_full_2011 ; Rahman_Roughness_engineered_valley_orbit .

By loading the surface roughness profile into NEMO3D, we are able to simulate quantum dot wavefunctions with atomic resolution under the correct local symmetries induced by the disordered surface. A limitation of the model is that we cannot simulate atomically disordered Si-O bonds in the oxide. Considering the various geometrical permutations of such bonds in an amorphous solid, it becomes a computationally challenging problem for large scale simulations. However, the VCA model does replicate the bulk electrical properties of SiO₂. Despite this limitation, we are able to simulate effectively a Si-SiO₂ interface that in general is hard to describe from a tight binding approach.

Valley Phase calculation with atomistic simulation The atomistic tight binding software outputs the electron densities of the two valley states ($\|\Psi_{v_{-}}(x,y,z))\|^{2}$ and $\|\Psi_{v_{+}}(x,y,z)\|^{2}$) Saraiva2009 . Their difference may be interpreted in terms of an envelope function $\Psi_{\rm Env}$ multiplied by cosine oscillations with the wavevector of the conduction band minima $k_{0}=0.822\pi/a_{0}$

A colour plot of $\psi_{\textbf{Osc}}$ is shown in Fig. 2e. To compute the valley phase we average $\psi_{\text{Osc}}$ over the $x-y$ plane and perform a Fourier transform in the z-axis. There is a distinctive peak with frequency $2\kappa_{0}$, as it is expected for valley oscillations. The valley phase $\phi_{v}$ is obtained from the complex phase of the transform at this frequency.

G matrix computation from atomistic tight-binding

Atomistic simulations with NEMO3D include intrinsically spin-orbit interactions in the $sp^{3}d^{5}s^{*}$ 20-band model for silicon Klimeck2007 . For all simulations, we set the magnetic field magnitude amplitude to $1$T and vary only the field orientation $\hat{B}$. The spin g-factor is independent of the magnitude of the magnetic field as long as the valley splitting is non-degenerate with the Zeeman splitting Bourdet2018 ; gilbert_-demand_2023 .

For any magnetic field $\mathbf{B}$, the spin part of the Hamiltonian of the system is

where $\mathbb{G}$ is the $\mathbb{G}$-matrix and $\mathbf{B}_{eff}=\frac{1}{g_{0}}\mathbb{G}\mathbf{B}$ is the effective magnetic field after including spin-orbit effects. The atomistic tight binding software outputs the Zeeman splitting $E_{\text{Zeeman}}=g_{0}\mu_{B}\|\mathbf{B}_{eff}\|$ and also the full wavefunction of the ground state $\Psi_{\downarrow}$ written in a base of atomic positions, orbitals and spins $|\textbf{R},\alpha,s\rangle$. From this we can estimate the mean spin vector $\langle\mathbf{\sigma_{\downarrow}}\rangle=\langle\Psi_{\downarrow}|\mathbf{\sigma}|\Psi_{\downarrow}\rangle$. This spin aligns anti-parallel to the effective field $g_{0}\mathbf{B}_{eff}$ by definition. In total, we have obtained the magnitude and orientation of the vector $g_{0}B_{eff}$. If we perform this computation for three linearly independent magnetic fields $\mathbf{B}_{1},\mathbf{B}_{2},\mathbf{B}_{3}$, we will obtain the effective fields $\mathbf{B}_{\text{eff},1},\mathbf{B}_{\text{eff},2},\mathbf{B}_{\text{eff},3}$. Then, the linear system $\mathbf{B}_{\text{eff},i}=\mathbb{G}\mathbf{B}_{i}$ can be inverted to compute $\mathbb{G}$.

Then

The oxide interface is defined by a smooth step of $V_{{\rm Si-SiO}_{2}}=3.1eV$ with the function $\sigma(z)=\frac{1}{1+e^{-4((z-z_{s}(\mathbf{r}))/a_{0}}}.$ Where $a_{0}=0.543$nm is the silicon lattice parameter. $z_{s}(\mathbf{r})$ defines $z$ coordinate of the rough surface at the position of $\mathbf{r}$ projected in the xy-plane $(r_{x},r_{y})$.

Path integral simulation of exchange coupling Hundreds of realizations of two-electron paths are sampled with path integral Monte Carlo (PIMC) with a Metropolis algorithm to minimize the partition function $\mathbb{Z}=e^{-S/\hbar}$ Ceperley1995 , where $S$ is the total action $S=\sum_{m=0}^{N_{t}}\tau H_{2e}(\mathbf{r}_{1}(m\tau),\mathbf{r}_{2}(m\tau))$. We simulate both, paths that remain in separate quantum dots for the entire simulation and paths that exchange a few times between the quantum dots. The action of electron paths that exchange between the two dots is higher than for non-exchanging paths by an amount $\Delta S$ which allows to compute the exchange energy as $J=\frac{2}{\beta}e^{-\Delta S/\hbar}$ Pedersen2010 . As the estimates of the operators are computed from average sampled path realizations, the method provides natural error bars determined by the standard deviation of $-\Delta S$. To adapt the method to 3D electrons in MOS quantum dots we had to do modifications to the main algorithm that are detailed in a separate paper (see Ref. CifuentesPIMC2023 ). For this simulations we used paths with $N_{t}=8000$ time slices and $\beta\hbar=4$ps. The path integral method converges at these values as shown in Ref. CifuentesPIMC2023 .

Voltage Offset Deviation (VOD)

For any parameter $\sigma$ this metric estimates the deviation of voltage offsets that are needed to bring each variable to the average value $\langle\sigma\rangle$. For the main parameters discussed in this paper we obtained:

1. 1.

   Valley Splitting: $\text{VOD}(\text{VS})=$ 0.58 V

2. 2.

   g-factor [110]: $\text{VOD}(g_{[110]})=$ 9.1 V
   Maximum among in-plane B-field angles

3. 3.

   g-factor [110]: $\text{VOD}(g_{[100]})=$ 0.23 V
   Minimum among in-plane B-field angles

4. 4.

   Valley Splitting: $\text{VOD}(J)=$ 0.09 V

The values of $\sigma$ are obtained from roughness variability distributions. Only tunings in the range of hundreds of mV can be performed in real experiments. If the VOD is way higher than that, it means that it is not possible to tune $\sigma$ to the same value for all qubits.

Estimating the tunability of qubit parameters To compute the tunability of each qubit $d\sigma/dV$ we use simulations for different voltage tunings for the specific gates have a significant impact on each one of the variables. For the exchange ($J$) we focus on the J-gate. In contrast, single qubit parameters like g-factors and valley splittings can be tuned with more than one gate. In here we assume that the action is performed with a top gate and an additional lateral gate so that ${d\sigma/dV=d\sigma/dV_{\rm Top}+d\sigma/V_{\rm Lat}}$ . At the same time, each of these tunings is divided as

Estimates for the impact of different gates on single quantum dot parameters can be found in Supplementary Table 3. The values for $d\sigma/dE_{z}$ and $d\sigma/dx$ where simulated from small changes in the simulation parameters as in Figs. 2f and 3 c-d.

The data for the figures generated in this study are provided in the Supplementary Information/Source Data file. Additional data and corresponding analysis code relevant to verify the results obtained here and generate the plots shown in this paper have been deposited in the Figshare database under accession code DOI: https://doi.org/10.6084/m9.figshare.23507439.

The path integral Monte Carlo algorithm to estimate exchange interactions with is described in Ref. CifuentesPIMC2023 . The algorithm to perform atomistic tight-binding simulations is described in Ref.Klimeck2007 . The original code used for processing and generating all the plots in this paper have been deposited in the Figshare database under accession code DOI: https://doi.org/10.6084/m9.figshare.23507439.

We acknowledge support from the Australian Research Council (FL190100167 and CE170100012), the US Army Research Office (W911NF-23-1-0092), the US Air Force Office of Scientific Research (FA2386-22-1-4070), and the NSW Node of the Australian National Fabrication Facility. The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the Army Research Office, the US Air Force or the US government. The US Government is authorized to reproduce and distribute reprints for Government purposes notwith- standing any copyright notation herein. J.Y.H., S.S., M.K.F. and J.D.C. acknowledge support from the Sydney Quantum Academy. This project was undertaken with the assistance of resources and services from the National Computational Infrastructure (NCI), which is supported by the Australian Government and includes computations using the computational cluster Katana supported by Research Technology Services at UNSW Sydney.

W.H.L. and F.E.H. fabricated the devices, with A.S.D.’s supervision, on isotopically enriched 28Si wafers supplied by K.M.I. (800 ppm), N.V.A., H.-J.P., and M.L.W.T. (50 ppm). T.T., W.G, J.H , E.V, R.C.C.L, S.S provided measurements of the qubit devices with supervision of C.H.Y., A.S., A.L and A.S.D.  The data was gathered by J.D.C., who used it to estimate the CMOS variability in comparison with one qubit and two qubit simulations.  J.D.C performed the one qubit simulations with atomistic tight binding in NEMO 3D with the supervision of R.R and A.S.  For two qubit systems, J.D.C performed the exchange simulations with an in-house Path Integral Monte Carlo code developed by P.Y.M , F.S and J.D.C. with the supervision of A.S.  D.O coded an initial framework for the fractal analysis of the Si/SiO₂ interface that was later improved by J.D.C with the supervision of A.S.  F.E.H took the TEM images of the devices.  J.Y.H and D.D coded the software that designs the digital twin of the devices with the supervision of C.C.E. and A.S.  These digital twins where imported for electrostatic simulations in COMSOL by J.D.C. and D.D. with the supervision of C.C.E. M.K.F provided theoretical support at all stages. J.D.C and A.S wrote the manuscript, with the input from all authors.

A.S.D. is the CEO and a director of Diraq. WTT, WG, EV, CCE, AL, CHY, WHL, FEH, AS, and ASD declare equity interest in Diraq Pty Ltd. Other authors declare no competing interests.

Supplementary information