Bounds to electron spin qubit variability for scalable CMOS architectures

Spins of electrons in silicon MOS quantum dots combine exquisite quantum properties and scalable fabrication. In the age of quantum technology, however, the metrics that crowned Si/SiO2 as the microelectronics standard need to be reassessed with respect to their impact upon qubit performance. We chart spin qubit variability due to the unavoidable atomic-scale roughness of the Si/SiO2 interface, compiling experiments across 12 devices, and develop theoretical tools to analyse these results. Atomistic tight binding and path integral Monte Carlo methods are adapted to describe fluctuations in devices with millions of atoms by directly analysing their wavefunctions and electron paths instead of their energy spectra. We correlate the effect of roughness with the variability in qubit position, deformation, valley splitting, valley phase, spin-orbit coupling and exchange coupling. These variabilities are found to be bounded, and they lie within the tolerances for scalable architectures for quantum computing as long as robust control methods are incorporated.


Introduction
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 [1].In particular, the similarity between quantum dots defined by gate electrodes on top of a silicon/silicon dioxide interface and the MOS-FET transistors in materials, design, and fabrication enables the integration of manufacturing techniques exclusive to semiconductor foundries onto the scaling of quantum processors [2].
Here, we specifically treat the case of such qubits formed in quantum dots at the Si/SiO 2 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 [3].These materials present their own complex challenges and advantages and impose significantly different architectural choices compared to the Si/SiO 2 interface, and hence are left out of this investigation.
Despite a relatively late start [4], the performance of silicon qubits has led to fidelity levels comparable with more well-established quantum technologies like superconducting or ion-trap qubits [5][6][7].The recent demonstration of repeatable high fidelity two-qubit operations across three nominally identical CMOS devices [8] 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 [9], 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 [10,11], 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 [2,12], and robust quantum control pulses that will be shared among several qubits [13,14].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 [21].Industrial foundries focus most of their efforts in addressing these issues [9,22,23], with a central role played by the choice of materials for substrate, dielectrics, gate metals, etc [20].
In the centre of the discussion is the dielectric interface where the electrons are confined [24,25].High fidelity silicon qubits have been measured in Si/SiGe and Si/SiO 2 heterostructures [3,[5][6][7][8].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 x Ge 1−x alloy.The Si/SiGe interface provides, in general, reduced levels of interfacial disorder compared to that of Si/SiO 2 [24].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 (blue) and T2 d (cyan), and the computer generated surface in a (see also Supplementary Fig. 6 and methods).The data is plotted as a function of λ = 2π q , where q is the wave-number.This allows us to compare λ with most relevant length scales, namely the silicon lattice parameter (a Si =0.543 nm), the dot diameter (10-15 nm), the double dot length (80-100 nm) and the lateral length of the simulation cell containing all 7x7 dots (500 nm).i, Average RMS of segments of length λ for the random surface generated numerically and for device T1 to T6. Error bars indicate the standard error estimated from repeated measurements across multiple TEM images.(see methods and Supplementary Table .2).Source data of figures e,h,i are provided in the Source Data file.
in Ref. [20].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 spinorbit coupling being studied here is masked by the presence of micromagnets which are commonly adopted in SiGe qubit architectures [3].
In this paper, we focus on qubits at the Si/SiO 2 interface.SiO 2 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 2 interface, in the bulk of SiO 2 or at the metal-oxide interface [22][23][24]26].Potential advantages are in the higher electrical tunability and compatibility with conventional CMOS technology, which benefits the integration with on-chip electronics [2].
The roughness of the Si/SiO 2 interface is one of the most critical sources of disorder for CMOS spin qubits [27][28][29][30][31][32].A more recent paper indicated that a second source -charged impurities in the oxide -could dominate over interface roughness [33], 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  [8,15] Table 1 | List of devices used in qubit measurements.The devices are identical except by differences indicated in here.The data in device A is taken at three electronic configurations: (1,1) -for 1 electron under gate 1 and 1 electron under gate 2, (3,1) and (1,3).The double quantum dots are formed under the gates P1, P2 or P2, P3 depending on the device.In one of the devices the qubits were driven magnetically with a dielectric resonator instead of an antenna [16,19].
A vector magnet enabled rotations of the magnetic field for the measurements in two of the devices.The next column refers to the gate material.We use a combination of palladium and atomic layer deposition (ALD) alumina for some devices and for others we form the gates with aluminium and isolate them with thermally formed alumina.Differences between these situations are discussed in [20].The last column shows the publications associated to each device .Device F is the only one with a single dot configuration instead of double dot.We use this device only to provide additional data on the valley splitting.
focus on spin qubits driven magnetically [4,16], 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 interfaceinduced 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 2 roughness on qubit properties at different lengthscales [34] (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).

Results
Si/SiO 2 roughness.The Si/SiO 2 interface has an intrinsic fractal structure [34], that has not been considered in previous variability studies ( [31][32][33]) .Our decision to simulate a 7×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 (∼ 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 2 (Fig. 1.b) is based on the roughness observed in TEM images (Fig. 1de) [34,35], 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 [35], quantified through its power spectral density (PSD) as a function of the in-plane correlation length scale λ [34].Our work incorporates a theoretical study of the realistic scale-dependent fractal structure of the roughness [34] 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 PSD 1D (λ) ∝ 2π λ −1−2H .We estimate a Hurst exponent of H=0.28 [36] (details in Supplementary Fig. 6 and methods section).The root-mean-square roughness also scales up with the lateral region λ as RMS (λ) ∝ 2π λ −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 [23].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 2 interface in our simulations.
In this paper, the quantum dots are simulated in a 7×7 grid array.Other architectures are also being explored [37], 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 [2].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 [38].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 [39] or would require quantum information to be shuffled to the edges of the array for readout [40].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 × 450 nm Si/SiO 2 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 2 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 2 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 [41,42].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 [mT] Ref. [43] Spin  z  is the wavevector of the conduction band minima in the silicon crystal.f, Valley splitting distribution of the 49 quantum dots versus electric field.Box plots indicate median (middle line), 25th, 75th percentile (box) and 5th and 95th percentile (whiskers) as well as outliers (single points).We compare our results with experimental data measured in device F and with measurements in Ref. [43].Error bars indicate the standard deviation for the measured value.The electric fields are obtained from COMSOL simulations.The inset figure shows the correlation between the logarithm of the valley splitting versus the centre of the dot in the z axis.The valley splitting is reconverted to magnetic field units in the right axis to compare it with the Zeeman splitting at different fields.g, Distribution of valley phases versus the valley splitting.For convenience, we define ϕv = 0 • as the point with the highest density of valley phases.The colour code represents the value of Ez as in f.Source data of f-g are provided in the Source Data file.
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 2 interface.The conduction band valleys along ±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 µeV to a few meV, compatible with observed experimental values in CMOS devices [32,43].
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 (≈28 meV 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 largescale 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 [44].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. [45].
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 [46,47].Fig. 2g shows the valley phases across the 49 simulated dots, revealing that dots with larger valley splittings (typically above 300 µ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 ∼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 Zeeman /B 0 = gµ 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 2 interfaces -Rashba (α) and Dresselhaus (β) [30].These two can be experimentally differentiated by measuring the g-factor dependence on the in-plane magnetic field orientation φ [29].The dependence is sinusoidal with the form g (ϕ) ≈ g 0 + α + β sin (2ϕ), 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 ( [110] [110] [110] A (1,1) Atomistic Simulations Simulations Fig. 3 | Variability of the spin-orbit coupling.a-b, Comparison between g-factor variability in atomistic simulations and measurements in devices A to E under a varying magnetic field angle.On each device and configuration we measured two qubits.In a we compare the difference between the Larmor frequencies of the two-qubit qubits measured in each device (g 1 − g 2 )µ B vs the differences between the frequency of neighbouring dots simulated atomistically (Fig. 1b).The marker for experimental data is associated with the device (A to E).In b, we compare the top gate Stark shift dg/dV measured in the two qubits of each device, with atomistic simulations of the dots (methods).Two qubits in a device have a different Stark shift dg/dV due to variations in surface roughness.Filled markers represent the data for the first qubit and empty makers for the second qubit (e.g. the empty purple triangles represents the Stark shifts measure on the second qubit of device E). c-d, Distribution of simulated Dresselhaus β (c) and Rashba α (d) spin-orbit terms versus vertical electric field Ez.Box plots indicate median (middle line), 25th, 75th percentile (box) and 5th and 95th percentile (whiskers) as well as outliers (single points).e-f, Schematic table showing that the sinusoidal dependence of the g-factors versus in-plane magnetic angle (e) follows from the anisotropy in the silicon lattice near the interface shown in (f ).In an ideal flat surface, the border of silicon must end in one of the two possible sublattices A (black circles) or B (gray circles) and the interface looks different when observed from the [110] and the 1 10 lattice orientations.In a realistic rough surface the border is a mixture of both A and B sub-lattice terminations, which explains the observed g-factor variability g, Distribution of qubit frequencies for two magnetic field orientations: [110] and [100].Qubits affected by near-valley degeneracies in the simulation are not included in this figure.The bars show an estimate for the maximum gate tunability of the g-factors with the top gate and a lateral gate (methods).The typical range of tunability for this gate is about 0.2 V for these devices.A higher potential bias could induce a charge transition, thus ruining the two-qubit system.h, Zoom to the [100] data in g.Source data are provided in the Source Data file 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 spinorbit effect and its variability -with the exception of the configuration (1,1) in device A (Fig. 3.a).The Rashba coefficient α is typically one order of magnitude smaller than β (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 α and β 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 [15,48,49], 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 β 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 [31].Fig. 3e shows that in this case the Dresselhaus parameter β 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 β 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 [4].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 [13,14].Qubits will be driven in resonance with the global field (forming a dressed qubit [50]), 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 ≈100 mT [51].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. [13]).The Stark shift tunability also decreases, such that control strategies might need to cope with the inability to tune qubits completely out of resonance [13].
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 [22,26].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 [18,51,52].
Previous works have modeled charged traps as negative electron charges e − directly at the Si/SiO 2 interface [33].However, factors such as the positive correlation between the SiO 2 thickness with charge mobility in Hall bar devices [23,53], and measured large Dingle ratio [22] 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 2 with density −4 × 10 10 nm −1 [23] (see Fig. 4a).Each trap is assumed to lead to a deformation of the potential of We set Ez = 28 meV nm −1 in these simulations (Fig. 2f ).In both cases, simulations are performed with the same surface profile (Fig. 1b).e, Dispersion of dot centres under the presence of interface roughness and charged traps.where the trap at r t and ϵ Si/SiO 2 = 0.5(ϵ Si + ϵ SiO2 ) = 7.5ϵ 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 2 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 [54].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 [46].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 [56,57].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 [58].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 2 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 [59].The full method is described in Ref. [55] 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 P1 +V 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 [60].
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. [8], observing larger control rates and variability due to the absence of the ALD oxide and irregularities in the gate structure [20].
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 [55].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.

Discussion
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 2 interfaces.The main conclusion is that most of the qubit variability in current devices is explained by the roughness of the Si/SiO 2 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 gfactor variations that can affect important qubit metrics, such as the phase coherence T * 2 [18,29].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 [61].
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 spinorbit 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 [62].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 2 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 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 [13].These results outline the minimum demands for a CMOS spin qubit architecture.

Methods
Fabrication process.The SiO 2 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 [63].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 [64,65].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 [66].
Image analysis to filter the Si/SiO 2 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 2 regions.These convolutions are defined by four matrices M cc , M sc , M cs , M ss with dimensions n x × n y and entries M ss ij = sin(2πi/n x ) sin(2πj/n y ) where i ∈ {1 . . .n x } and j ∈ {1 . . .n y }.The result in Fig. 13b is obtained after taking the average of the images generated by four convolutions M cc * I, M sc * I, M cs * I, M ss * I , where 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 ∈ (0, 1) to segment both regions (Fig. 13c).We tune the parameters N and T until the fit clearly corresponds to the Si/SiO 2 interface in the original TEM (Fig. 13d).
Si/SiO 2 characterization of interface from TEM images.We characterize the Si/SiO 2 roughness from TEM images of devices T1 to T6 (see Supplementary Fig. 6) using a similar procedure to [35].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 2 roughness is characteristic of a fractal self-affine interface [34,36].We obtained an average Hurst exponent of H = 0.28 ± 0.2 and a roughness amplitude parameter of C 0 ≈ 1.4nm 3 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 2 interface from TEM images.To compute the scaling of the RMS(λ) over a certain interface length characterised by the wavelength of its Fourier decomposition λ, we took separated segments of the fitted line surface from TEM image of width λ and computed the average RMS for each device.If a TEM has a width of L =40 nm, we would compute the RMS(λ) for λ = L/2, L/4, . . .until reaching the atomic scale a 0 [34].We can select more subsegments when λ is small, which explains the smaller uncertainty of the average RMS for small λ.We found that the RMS scales with an exponent of −H as which is expected for the PSD profile in equation ( 6) [36].
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 [67], 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 electronbeam 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.
Electrostatic simulations and quantum dot model.We import the digital twin device into COMSOL Multiphysics and perform electrostatic potential simulations with the integrated Poisson solver.We fit this to a harmonic model with a vertical electric field where (x c , y c ) is the centre of the parabolic potential, E z is the electric field in the z-axis (typically 8 to 40 meV nm −1 ) and c x , c y are the lateral curvatures (approximately 0.3 meV nm −2 ).We simulate potential sweeps from different gates to characterize their impact on the quantum dots (see Supplementary Fig. 7).The harmonic model allows to transform these actions to more intuitive parameters such as shifts in the electric confinement, dot movement, ellipticity, etc.A summary of this impact is included in Supplementary Table 3.
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-orbitinteractions [68].In all simulations the magnetic field magnitude is set to 1T, 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 2 (Si), as observed in Figure 2a.The SiO 2 region is modeled with a sp 3 tight-binding (TB) model under a virtual crystal approximation (VCA).The SiO 2 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 0 .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 2 atom.This is a standard way to model alloyed (SiGe) or disordered (SiO 2 ) materials under the VCA in the atomistic TB technique.At the interface region whether an atom is marked as a Si atom or SiO 2 atom then creates the atomistic disorder profile.The details of the model with parameter values can be found in [41,69].
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 2 .Despite this limitation, we are able to simulate effectively a Si-SiO 2 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 (∥Ψ v− (x, y, z))∥ 2 and ∥Ψ v+ (x, y, z)∥ 2 ) [70].Their difference may be interpreted in terms of an envelope function Ψ Env multiplied by cosine oscillations with the wavevector of the conduction band minima k A colour plot of ψ Osc is shown in Fig. 2e.To compute the valley phase we average ψ Osc over the x − y plane and perform a Fourier transform in the z-axis.There is a distinctive peak with frequency 2κ 0 , as it is expected for valley oscillations.The valley phase ϕ v is obtained from the complex phase of the transform at this frequency.
G matrix computation from atomistic tightbinding.Atomistic simulations with NEMO3D include intrinsically spin-orbit interactions in the sp 3 d 5 s * 20-band model for silicon [68].For all simulations, we set the magnetic field magnitude amplitude to 1T and vary only the field orientation 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 [15,48].
For any magnetic field B, the spin part of the Hamiltonian of the system is where G is the G-matrix and B ef f = 1 g0 GB is the effective magnetic field after including spin-orbit effects.The atomistic tight binding software outputs the Zeeman splitting E Zeeman = g 0 µ B ∥B ef f ∥ and also the full wavefunction of the ground state Ψ ↓ written in a base of atomic positions, orbitals and spins |R, α, s⟩.From this we can estimate the mean spin vector ⟨σ ↓ ⟩ = ⟨Ψ ↓ |σ|Ψ ↓ ⟩.This spin aligns anti-parallel to the effective field g 0 B ef f by definition.In total, we have obtained the magnitude and orientation of the vector g 0 B ef f .If we perform this computation for three linearly independent magnetic fields B 1 , B 2 , B 3 , we will obtain the effective fields B eff,1 , B eff,2 , B eff,3 .Then, the linear system B eff,i = GB i can be inverted to compute G.
A typical G-matrix obtained from atomistic simulations is where the basis is aligned with the lattice orientations {[100], [010], [001]}].In here, g 0 ≈ 1.9937 is the bulk g-factor, as obtained from the asymptotic behavior of the simulations at low electric field.The entries α ′ ∼ −10 −3 and β ′ ∼ ±10 −2 determine the in-plane spin-orbit-coupling and the parameters g 13 ∼ ±10 −3 , g 23 ± ∼ 10 −3 and g 33 ∼ 2.00192 − O(10 −4 ) describe the out of-plane components.The two remaining entries were calculated to be smaller than 10 −5 at all cases and hence approximated to 0. To obtain the g-factor at any magnetic field orientation we compute g(r) = ∥Gr∥.If r(ϕ) = [cos ϕ sin ϕ 0] T is an in-plane normal vector Supplementary information.2 Table of devices where transmission electron microscopy images (TEM) were taken.The TEMs capture the Si/SiO2 interface, showing the interface roughness.All oxides were grown under the same conditions except for device T4 consisting of a 7.5 nm SiO2 layer grown on isotopically purified 50pmm 29 Si.The PSD was obtained from TEMs with varying numbers and ranges, with more TEMs yielding a higher degree of precision in the PSD and RMS estimate.

Fig. 1 |
Fig. 1 | Modeling CMOS spin qubits.a, Example scaled architecture of a 49 qubit device.b, The quantum dots are formed below the computer generated rough surface.c, Model of the three-dot devices measured in this paper.The metallic gates are coloured by their order of deposition in different layers.d-g, Comparison between device cross sections in TEM images and computer model.d, TEM image of device T1, showing a cross section of the device located approximately at the position of the violet rectangular region in c. e, TEM of device T2 with a focus on the silicon oxide interface.We highlight in d a square region with the same size.f, Cross-section of model at the green rectangular region in c. g, Atomistic simulation showing the electronic wavefunction of a quantum dot below rough Si-SiO 2 .h, Average power spectral density (PSD) of the Si/SiO 2 interface comparing the interface from transmission electron microscopy (TEM) images of device T1(blue) and T2 d (cyan), and the computer generated surface in a (see also Supplementary Fig.6 and methods).The data is plotted as a function of λ = 2π q , where q is the wave-number.This allows us to compare λ with most relevant length scales, namely the silicon lattice parameter (a Si =0.543 nm), the dot diameter (10-15 nm), the double dot length (80-100 nm) and the lateral length of the simulation cell containing all 7x7 dots (500 nm).i, Average RMS of segments of length λ for the random surface generated numerically and for device T1 to T6. Error bars indicate the standard error estimated from repeated measurements across multiple TEM images.(see methods and Supplementary Table.2).Source data of figures e,h,i are provided in the Source Data file.

Fig. 2 |
Fig.2| Quantum dot variability.a, To simulate a Si/SiO 2 interface atomistically, we use a virtual lattice approximation (see methods).This allows us to emulate the electronic properties of SiO 2 -which does not have a regular lattice structure -in a simulated material with the same lattice structure as silicon.We then define the interface between the silicon lattice and its oxide with a simulated rough surface dividing the atomic sites between the two materials.b, Three-dimensional visualization of a CMOS quantum dot at the Si/SiO 2 interface simulated atomistically.The blue arrow represents the vector spin ⟨σ⟩ averaged across the spin-orbitals of all the atoms in the quantum dot.c, In-plane visualization of the variability in the 7 quantum dots inside the purple rectangle in Fig.1b.The 5 nm diameter cyan circle is a static reference to compare the wavefunctions in different simulations.Black asterisks represent the center of each quantum dot ⟨r⟩ = ⟨ψ|r|ψ⟩.d, Variability distribution of dot centres.e, Visualization of the valley oscillations parallel to the [001] lattice orientation.(see Methods section).The oscillation wavelength is 4π k 0 , where k 0 = 0.82 2π a 0 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 consideredc

Fig. 4
Fig.4Impact of charged impurities in quantum dots below a rough Si/SiO 2 a, A random distribution of negative charge impurities with typical industry-level densities of ∼ 4 × 10 10 nm −1[23].b, Variability of the quantum dot wavefunctions of the same dots in Fig.2c.There is a negative trap close to dot number 2. c-d, Atomistic simulations of valley splittings (b) and g-factors (c) of each quantum dot in the 49 dot array, comparing simulations with and without charged traps.We set Ez = 28 meV nm −1 in these simulations (Fig.2f).In both cases, simulations are performed with the same surface profile (Fig.1b).e, Dispersion of dot centres under the presence of interface roughness and charged traps.f Comparison between the amplitude of the dispersion of dot centres in both scenarios.Box in all figures indicate median (middle line), 25th, 75th percentile (box) and 5th and 95th percentile (whiskers) as well as outliers (single points).Source data of a,c-f are provided in the Source Data file.

f
Fig.4Impact of charged impurities in quantum dots below a rough Si/SiO 2 a, A random distribution of negative charge impurities with typical industry-level densities of ∼ 4 × 10 10 nm −1[23].b, Variability of the quantum dot wavefunctions of the same dots in Fig.2c.There is a negative trap close to dot number 2. c-d, Atomistic simulations of valley splittings (b) and g-factors (c) of each quantum dot in the 49 dot array, comparing simulations with and without charged traps.We set Ez = 28 meV nm −1 in these simulations (Fig.2f).In both cases, simulations are performed with the same surface profile (Fig.1b).e, Dispersion of dot centres under the presence of interface roughness and charged traps.f Comparison between the amplitude of the dispersion of dot centres in both scenarios.Box in all figures indicate median (middle line), 25th, 75th percentile (box) and 5th and 95th percentile (whiskers) as well as outliers (single points).Source data of a,c-f are provided in the Source Data file.

Fig. 5 |
Fig. 5 | Variability of the exchange coupling.a, Diagram depicting exchange interactions between neighbouring CMOS spin qubits in the presence of interface roughness and negative charged traps.The green and yellow paths between the quantum dots represent the two electrons exchanging.These are simulated with a path integral Monte Carlo algorithm [55].The interstitial exchange gate (J1) controls the interdot barrier.At high J-gate voltage the spin interaction is high enough and the two qubit frequencies split.This is measured in device A at the (3 electron, 1 electron) configuration, b .c, Representation of the 3D path of the two electrons moving inside a double quantum dot generated with path integral Montecarlo (PIMC).The electric potential configuration is plotted in gray and the electron density appears in the xy-plane.d, Exchange versus J-gate detuning, comparing PIMC simulations in the presence of interface roughness (Fig. 1.b) and charge traps (Fig. 4.a) with data from devices A and C. The colors of the simulations (small circular markers with error bars) represent the input potential values for the J-gate (1 V, 1.5 V, 2 V).Error bars are associated to the standard error of Path interal Monte Carlo simulations[55].The inset figure shows a cut over the x-axis of the potential configuration for each value of J.This potential is obtained from finite element simulations in COMSOL MULTIPHYSICS (see Supplementary Fig.7 and Methods).e Double dot wavefunctions simulated with path integral Monte Carlo (PIMC) for V J = 1 V and V J = 2 V[55].f, Correlation between exchange coupling and inter-dot distance in PIMC simulations.g, Histogram of the exchange controllability rates (dJ/dV J ) of device simulations, compared with experimental values (see markers).Source data of b,d,f,g are provided in the Source Data file.

Fig. 7 |
Fig.7| Potential simulations and gate-impact on quantum dots.a, Horizontal view of the 3D device model that we input in COMSOL for potential simulations.The green square shows the region where the dots are formed.b, Potential landscape of the device simulated in COMSOL.The set of gate potentials is derived from experiments.A double quantum dot is isolated in the green square region bellow the gates P 1 and P 2.c, Zoom into the potential profile at the green region in b.Single dots are formed inside the white rectangles.We fit the potentials inside these regions to the harmonic model V (x, y, z) = cx(x − xc) 2 + cy(y − yc) 2 + zEz.d,f, Evolution of the potential profile over the cyan line in c under the tuning of gates P 1 (d) and J1 (f ).e,g, Characterization of the impact of gates P 1 (e) and P 2 (g) on each quantum dot.We evaluate this impact over 5 variables: Displacement of the dot position from the mean (δxc, δyc), transversal electric field Ez and curvatures (cx, cy).The numbers inside the plots show the slope of the curve.Their units depend on the variable evaluated, for instance , the unit of dxc/dV is [nm/V].

Fig. 9 |Fig. 10 |
Fig.9| Dependence of qubit parameters on the surface RMS: a-b, We generated two additional surfaces with higher (a) and lower (b) RMS that the one used in the paper to understand the potential benefits of improving the surface quality.c-d, 1D RMS and 1D PSD of the original surface and the new surfaces plotted in a and b.The profiles of devices T1 and T2 are also included in both figures for comparison with the dispersion of the measured data.e-f, Spin-orbit coupling dependence on the RMS.Dresselhaus values are slighly more dispersed for smoother interfaces as the results approach to the flat surface limits.g-h, RMS vs valley splitting(g) and two-dot exchange coupling (h).Smoother surfaces lead to higher mean valley splittings g and to smaller variability in exchange coupling h.

Fig. 11 2 Fig. 12
Fig. 11 Quantum dot parameters simulated with atomistic tight-binding of each dot in the 7×7 qubit array: a Valley spliting (in logarithmic scale).b g-factor with magnetic field pointing to [110].The simulated rough surface is plotted behind the dots for reference.Each dot was simulated with atomistic tight-binding in a simulation cell of 40 nm ×40 nm containing a section of the global rough surface.Ez = 28 meV nm −1 in these simulations.

Fig. 13 a
Fig. 13 a High resolution TEM image of device T2.b-d, Steps to filter the edge of the Si-SiO 2 interface using the periodic pattern recognition method.b, Image created after periodic convolution.c, Discretization setting a threshold parameter on the intensity of the figure b. d, Original TEM with fit of the Si/SiO2 interface.