Atomic fluctuations lifting the energy degeneracy in Si/SiGe quantum dots

Electron spins in Si/SiGe quantum wells suffer from nearly degenerate conduction band valleys, which compete with the spin degree of freedom in the formation of qubits. Despite attempts to enhance the valley energy splitting deterministically, by engineering a sharp interface, valley splitting fluctuations remain a serious problem for qubit uniformity, needed to scale up to large quantum processors. Here, we elucidate and statistically predict the valley splitting by the holistic integration of 3D atomic-level properties, theory and transport. We find that the concentration fluctuations of Si and Ge atoms within the 3D landscape of Si/SiGe interfaces can explain the observed large spread of valley splitting from measurements on many quantum dot devices. Against the prevailing belief, we propose to boost these random alloy composition fluctuations by incorporating Ge atoms in the Si quantum well to statistically enhance valley splitting.

the spin coherence 24 . Furthermore, small valley splitting may affect Pauli spin blockade readout 25 , which is considered in large-scale quantum computing proposals 5,6 . Therefore, scaling up to larger systems of single-electron spin qubits requires that the valley splitting of all qubits in the system should be much larger than the typical operation temperatures (20−100 mK).
It has been known for some time that valley splitting depends sensitively on the interface between the quantum well and the SiGe barrier 26 . Past theoretical studies have considered disorder arising from the quantum well miscut angle 27 and steps in the interface [28][29][30][31][32] demonstrating that disorder of this kind can greatly decrease valley splitting in quantum dots. However, a definitive connection to experiments has proven challenging for a number of reasons. At the device level, a systematic characterisation of valley splitting in Si/ SiGe quantum dots has been limited because of poor device yield associated with heterostructure quality and/or device processing. At the materials level, atomic-scale disorder in buried interfaces 33 may be revealed by atom-probe tomography (APT) in three-dimensions (3D) over the nanoscale dimensions comparable to electrically defined quantum dots. However, the current models employed to reconstruct in 3D the APT data can be fraught with large uncertainties due to the assumptions made to generate the threedimensional representation of the tomographic data 34 . This results in limited accuracy when mapping heterointerfaces 35 and quantum wells [36][37][38] . These limitations prevent linking the valley splitting in quantum dots to the relevant atomic-scale material properties and hinder the development of accurate and predictive theoretical models.
Herein we solve this outstanding challenge and establish comprehensive insights into the atomic-level origin of valley splitting in realistic silicon quantum dots. Firstly, we measure valley splitting systematically across many quantum dots, enabled by high-quality heterostructures with a low disorder potential landscape and by improved fabrication processes. Secondly, we establish a new method to analyse APT data leading to accurate 3D evaluation of the atomiclevel properties of the Si/SiGe buried interfaces. Thirdly, incorporating the 3D atomic-level details obtained from APT, we simulate valley splitting distributions that consider the role of random fluctuations in the concentration of Si and Ge atoms at each layer of the Si/SiGe interfaces. By comparing theory with experiments, we find that the measured random distribution of Si and Ge atoms at the Si/SiGe interface is enough to account for the measured valley splitting spread in real quantum dots. Based on these atomistic insights, we conclude by proposing a practical strategy to statistically enhance valley splitting above a specified threshold as a route to making spin-qubit quantum processors more reliable-and consequently-more scalable.

Results
Material stacks and devices Figure 1 overviews the material stack, quantum dot devices, and measurements of valley splitting. To increase statistics, we consider two isotopically purified 28 Si/Si 0.7 Ge 0.3 heterostructures (quantum wells A and B) designed with the same quantum well width and topinterface sharpness (Methods), which are important parameters determining valley splitting 23,26  quantum well A (Fig. 1a) has a sharp 28 Si → Si-Ge heterointerface at the top and a diffused Si-Ge → 28 Si heterointerface at the bottom, whereas in quantum well B (Fig. 1b) the growth process was optimized to achieve sharp interfaces at both ends of the quantum well. These heterostructures support a two-dimensional electron gas with high mobility and low percolation density (Supplementary Figs. 1 and 2), indicating a low disorder potential landscape, and high-performance qubits 10,39 with single-and two-qubit gates fidelity above 99% 10 . We define double-quantum dots electrostatically using gate layers insulated by dielectrics (Methods). A positive gate voltage applied to plunger gates P1 and P2 (Fig. 1c) accumulates electrons in the buried quantum well, while a negative bias applied to other gates tunes the confinement and the tunnel coupling between the quantum dots Q1 and Q2. All quantum dots in this work have plunger gate diameters in the range of 40-50 nm ( Fig. 1d and Supplementary Table 1), setting the relevant lateral length scale for atomic-scale disorder probed by the electron wave function.

Valley splitting measurements
We perform magnetospectroscopy measurements of valley splitting E v in dilution refrigerators with electron temperatures of about 100 mK (Methods). Figure 1e shows a typical charge stability diagram of a double quantum dot with DC gate voltages tuned to achieve the few electron regime, highlighted in Fig. 1f. We determine the 2-electron singlet-triplet energy splitting (E ST ) by measuring the gate-voltage dependence as a function of parallel magnetic field B along the (0,1) → (0,2) transition ( Fig. 1g) and along the (1,1) → (0,2) transition ( Supplementary Fig. 4). In Fig. 1g, the transition line (black line) slopes upward, because a spin ↑ electron is added to form a singlet ground state S 0 . Alternatively, a spin down electron can be added to form a T −state, with a downward slope. A kink occurs when the S 0 -state is energetically degenerate with the T − -state, becoming the new ground state of the two-electron-system. From the position of the kink (B ST = 1.57 T) along the theoretical fit (red line) and the relation E ST = gμ B B ST , where g = 2 is the electron gyromagnetic ratio and μ B is the Bohr magneton, we determine E ST = 182.3 μeV for this quantum dot. E ST sets a lower bound on the valley splitting, E v ≥ E ST 21,40 . Due to small size, our dots are strongly confined with lowest orbital energy much larger than E ST ( Supplementary Fig. 3), similar to other Si/SiGe quantum dots 14,18,22 . Therefore, we expect exchange corrections to have negligible effects 40 and here take E v ≈ E ST .
Here we report measurements of E v in 10 quantum dots in quantum well A and 12 quantum dots in quantum well B (Supplementary Figs. 5 and 6) and compare the measured values in Fig. 1h. We observe a rather large spread in valley splittings, however we obtain remarkably similar mean values and two-standard-deviation error bars E v ± 2σ of 108 ± 55 μeV and 106 ± 58 μeV for quantum wells A and B, respectively. The quantum dots all have a similar design and hence are expected to have similar electric fields across the devices with a small influence on valley splitting under our experimental conditions. We argue that quantum wells A and B have similar E v ± 2σ because the electronic ground state is confined against the top interface, which is very similar in the two quantum wells.

Atom probe tomography
We now characterise the atomic-scale concentration fluctuations at the quantum well interfaces to explain the wide range of measured valley splittings with informed theoretical and statistical models. To probe the concentrations over the dimensions relevant for quantum dots across the wafer, we perform APT on five samples each from quantum wells A and B, with a field of view of approximately 50 nm at the location of the quantum well (Methods). First, we show how to reliably reconstruct the buried quantum well interfaces, then we use this methodology to characterise their broadening and roughness. Figure 2a shows a typical point-cloud reconstruction of an APT specimen from quantum well B. Each point represents the estimated   Supplementary Fig. 14). b, c Voronoi tessellation of the APT reconstructions for quantum wells A and B, respectively, and extracted isosurfaces corresponding to 8% Ge concentration. z is the average position of the 8% Ge concentration across these particular samples. We limit the lateral size of the analysis to ≈ 30nm × 30 nm, reflecting the typical lateral size of a quantum dot (Fig. 1d). d Average germanium concentration depth profiles across quantum wells A (magenta) and B (green). Shaded areas mark the 95% confidence interval over each of the sets of five APT samples. e Statistical analysis of the top interface width 4τ determined by fitting the data for quantum wells A (magenta) and B (green) to sigmoid functions. Thick and thin horizontal black lines denote the mean and twostandard-deviation error bars for the different APT samples. Dotted black lines show 4τ results from the HAADF-STEM measurements ( Supplementary Fig 13). position of an ionized atom detected during the experiment 34 . Qualitatively, we observe an isotopically enriched 28 Si quantum well, essentially free of 29 Si, cladded in a SiGe alloy. To probe the interface properties with the highest possible resolution allowed by APT and differently from previous APT studies on Si/SiGe 38 , we represent the atom positions in the acquired data sets in form of a Voronoi tessellation 41,42 and generate profiles on an x − y grid of the tessellated data, as described in Supplementary Note 2c. A sigmoid function ½1 + expðz À z 0 Þ=τ À138 is used to fit the profiles of each tile in the x − y grid. Here, z 0 is the inflection point of the interface and 4τ is the interface width. As the Voronoi tessellation of the data set does not sacrifice any spatial information, the tiling in the x − y plane represents the smallest lateral length scale over which we characterise the measured disorder at the interface. Note that we do not average at all over the z axis and hence maintain the inherent depth resolution of APT. We find that for tiles as small as 3 nm × 3 nm the numerical fitting of sigmoid functions to the profiles converges reliably. Although each tile contains many atoms, their size is still much smaller than the quantum dot diameter, and may therefore be considered to be microscopic. We use the sigmoid fits for each tile stack to visualise and further characterise the interfaces (Supplementary Figs. 8-10). Importantly, Ge concentration isosurfaces as shown in Fig. 2b, c are constructed by determining the vertical position for which each of the sigmoids reaches a specific concentration. Note, that we oversample the interface to improve the lateral resolution by making the 3 nm × 3 nm tiles partially overlap (Supplementary Note 2c).
In Fig. 2d, we show the average Ge concentration profile and measurement to measurement variations from the tessellated volumes (Supplementary Note 2b, c) of all samples for both quantum wells A and B. APT confirms HAADF-STEM results in Fig. 1a, b: quantum wells A and B have an identical sharp top interface and quantum well A has a broader bottom interface. Furthermore, the shaded colored areas in Fig. 2d reveal narrow 95% confidence levels, pointing to highly uniform concentration profiles when averaged across the wafer. Strong disorder fluctuations emerge at the much smaller tile length scale. In Fig. 2e we show for all samples of a given quantum well the interface width mean value with two standard deviations 4τ ± 2σ, obtained by averaging over all the tiles in a given sample. The results indicate uniformity of 4τ, and further averaging across all samples of a given heterostructure (μ 4τ , black crosses) yields similar values of μ 4τ = 0:85 ± 0:32 nm and 0.79 ± 0.31 nm for quantum wells A and B, consistent with our 4τ analysis from HAADF-STEM measurements (black dotted lines). However, the two-standard-deviation errors (2σ) of each data point can be up to 30% of the mean value 4τ.
To pinpoint the root cause of atomic-scale fluctuations at the interface, in Fig. 2f, g we utilize the 3D nature of the APT data sets, calculate, and compare the root mean square (RMS) roughness of the interfaces (solid green lines) as measured by APT on quantum well B to two 3D models (Fig. 2f, g) mimicking the dimensions of an APT data set. Both models are generated with random distributions of Si and Ge in each atomic plane (Supplementary Note 2d). The first model (solid black lines) corresponds to an atomically abrupt interface where the Ge concentration drops from~33.5% to 0% in a single atomic layer. It hence represents the minimum roughness achievable at each isoconcentration surface given the in-plane randomness of SiGe and the method to construct the interface. The second model (dashed black lines) is generated with the experimentally determined Ge concentration profile along the depth axis ( Supplementary Fig. 11). As shown in Fig. 2f, g, the roughness extracted from the second model fits well with the measured data, suggesting that the RMS roughness measured by APT is fully explained by the interface width and shape along the depth axis. Furthermore, as the deviation of each isosurface tile position from the isosurface's average position also matches that of the measured interfaces from the second model (Supplementary Movie 1) the APT data are consistent with a random in-plane distribution of Ge perpendicular to the interface in all data sets of quantum well B. For 2 out of 5 samples on quantum well A that we analyzed, we observe features that are compatible with correlated disorder from atomic steps (Supplementary Fig. 13). In the following, the alloy disorder observed in the APT concentration interfaces is incorporated into a theoretical model. As shown below, the calculations of valley splitting distributions associated with the 3D landscape of Si/SiGe interfaces can be further simplified into a 1D model that incorporates the in-plane random distribution of Si and Ge atoms.

Valley splitting simulations
We begin by considering an ideal laterally infinite heterostructure with no concentration fluctuations, and we denote the average Si concentration at layer l by x l . Due to the finite size of a quantum dot and the randomness in atomic deposition, there will be dot-to-dot concentration fluctuations. We therefore model the actual Si concentration at layer l by averaging the random alloy distribution weighted by the lateral charge density in the quantum dot, giving x d l = x l + δ x l , as described in Supplementary Note 3c. Here, the random variation δ x l is computed assuming a binomial distribution of Si and Ge atoms. We find that these fluctuations can have a significant impact on the valley splitting.
We explore these effects numerically using 1D tight-binding simulations. We begin with the averaged fitted concentration profiles obtained from the APT analysis in Fig. 2d, which enable us to directly measure the average Ge concentration in a given layer x l (Fig. 3a). The variance of the concentration fluctuations is determined by the size of the quantum dot, which we assume has an orbital excitation energy of ℏω = 4.18 meV and corresponding radius ffiffiffiffiffiffiffiffiffiffiffiffiffiffi _=m * ω p , as well as the average Si concentration x l . Here, m * is the effective mass of Si. Together, x l and the variance determine the probability distribution of weighted Si and Ge concentrations. Concentration profiles are sampled repeatedly from this distribution, with a typical example shown in Fig. 3b. The valley splitting is then determined from a 1D tight-binding model 43 . The envelope of the effective mass wavefunction ψ env (z) is shown in Fig. 3c (grey curve) for an electron confined in the quantum well of Fig. 3b. The procedure is repeated for 10,000 profile samples, obtaining the histogram of valley splittings shown in Fig. 3e. These results agree very well with calculations obtained using a more sophisticated three-dimensional 20-band sp 3 d 5 s* NEMO tight-binding model 44 (Supplementary Note 3b) and confirm that concentration fluctuations can produce a wide range of valley splittings. For comparison, at the top of Fig. 3e, we also plot the same experimental valley splittings shown in Fig. 1h, demonstrating good agreement in both the average value and the statistical spread. These observations support our claim that the valley splitting is strongly affected by composition fluctuations due to random distributions of Si and Ge atoms near the quantum well interfaces, even though the experiments cannot exclude the presence of correlated disorder from atomic steps in quantum dots.
Analytical methods using effective mass theory may also be used to characterise the distribution of valley splittings. First, we model the intervalley coupling matrix element 26 as Δ = R e À2ik 0 z l UðzÞ|ψ env ðzÞ| 2 dz, where k 0 = 0.82 × 2π/a 0 is the position of the valley minimum in the Si Brillouin zone, a 0 = 0.543 nm is length of the Si cubic unit cell, ψ env (z) is a 1D envelope function, and U(z) is the quantum well confinement potential. The intervalley coupling Δ describes how sharp features in the confinement potential couple the two valley states, which would otherwise be degenerate. In general, Δ is a complex number that can be viewed as the sum of two distinct components: a deterministic piece Δ 0 , arising from the average interface concentration profile, and a random piece δΔ, arising from concentration fluctuations. The latter can be expressed as a sum of contributions from individual atomic layers: δΔ = ∑ l δΔ l , where δΔ l is proportional to δ x l |ψ env ðz l Þ| 2 (see Methods). To visualize the effects of concentration fluctuations in Fig. 3c, we compute δΔ l using the randomized density profile of Fig. 3b (blue curve). We see that most significant fluctuations occur near the top interface, where |ψ env (z l )| and the Ge content of the quantum well are both large. In Fig. 3d we plot Δ values obtained for 10,000 quantum-well realizations using this effective mass approach. The deterministic contribution to the valley splitting Δ 0 (black dot) is seen to be located near the center of the distribution in the complex plane, as expected. However, the vast majority of Δ values are much larger than Δ 0 , demonstrating that concentration fluctuations typically provide the dominant contribution to intervalley coupling.
The total valley splitting is closely related to the intervalley coupling via E v = 2|Δ|, and therefore exhibits the same statistical behavior. In Fig. 3e, the orange curve shows the Rice distribution whose parameters are derived from effective-mass calculations of the valley splitting (see Methods), using the same concentration profiles as the histogram data. The excellent agreement between these different approaches confirms the accuracy of our theoretical techniques (Supplementary Note 3d).

Discussion
Based on the results obtained above, we now propose two related methods for achieving large valley splittings (on average), with high yields. Both methods are derived from the key insight of Fig. 3c: due to random-alloy fluctuations, the valley splitting is almost always enhanced when the electronic wavefunction overlaps with more Ge atoms. In the first method, we therefore propose to increase the width of the interface (4τ) as shown in Fig. 3f, since this enhances the wavefunction overlap with Ge atoms at the top of the quantum well. This approach is nonintuitive because it conflicts with the conventional deterministic approach of engineering sharp interfaces. The second method, also shown in Fig. 3f, involves intentionally introducing a low concentration of Ge inside the quantum well. The latter method is likely more robust because it can incorporate both deterministic enhancement of the valley splitting from a sharp interface, and fluctuation-enhanced valley splitting.
We test these predictions using simulations, as reported in Fig. 3g, where different colors represent different interface widths and the horizontal axis describes the addition of Ge to the quantum well. For no intentional Ge in the quantum well, as consistent with the heterostructure growth profile of our measured quantum dots, the calculations show a significant increase in the valley splitting with increasing interface width. Here, the narrowest interface appears most consistent with our experimental results (magenta marker), attesting to the sharp interfaces achieved in our devices. As the Ge concentration increases in the quantum well, this advantage is largely overwhelmed by concentration fluctuations throughout the well. A very substantial increase in valley splitting is observed for all concentration enhancements, even at the low, 5% level. Here, the light error bars represent 5-95 percentiles while dark bars represent 25-75 percentiles. At the 5% concentration level, our simulations indicate that >95% of devices should achieve valley splittings > 100 μeV. This value is more than an order of magnitude larger than the typical operation temperature of spin-  A). b Typical, randomized Ge concentration profile, derived from a. c Envelope function ψ env (z), obtained for the randomized profile in b (grey curve), and the corresponding concentration fluctuations weighted by the envelope function squared: δ x l |ψ env ðz l Þ| 2 (blue). Here, the wavefunction is concentrated near the top interface where the concentration fluctuations are also large; the weighted fluctuations are therefore the largest in this regime. d Distribution of the intervalley matrix element Δ in the complex plane, as computed using an effective-mass approach, for 10,000 randomized concentration profiles. The black marker indicates the deterministic value of the matrix element Δ 0 , obtained for the experimental profile in a. e Histogram of the valley splittings from tight-binding simulations with 10,000 randomized profiles. The same profiles may be used to compute valley splittings using effective-mass methods; the orange curve shows a Rice distribution whose parameters are obtained from such effective-mass calculations (see Methods). f Schematic Si/SiGe quantum well with Ge concentrations ρ W (in the well) and ρ b = ρ W + Δρ (in the barriers), with a fixed concentration difference of Δρ = 25%. g Distribution of valley splittings obtained from simulations with variable Ge concentrations, corresponding to ρ W ranging from 0 to 10%, and interface widths 4τ = 5 ML (red circles), 10 ML (blue triangles), or 20 ML (orange squares), where ML refers to atomic monolayers. Here, the marker describes the mean valley splitting, while the darker bars represent the 25-75 percentile range and the lighter bars represent the 5-95 percentile range. Each bar reflects 2000 randomized tight-binding simulations of a quantum well of width W = 120 ML. The magenta diamond at zero Ge concentration shows the average measured valley splitting of quantum well A. In all simulations reported here, we assume an electric field of E = 0.0075 V/nm and a parabolic single-electron quantum-dot confinement potential with orbital excitation energy ℏω = 4.18 meV and corresponding dot radius ffiffiffiffiffiffiffiffiffiffiffiffiffiffi _=m * ω p .
qubits and is predicted to yield a 99% readout fidelity 25 . This would represent a significant improvement in qubit yield for Si quantum dots. A recent report of SiGe quantum wells with oscillating Ge concentrations provides the first experimental evidence that intentionally placing Ge in the quantum well leads to significant variability and some of the highest recorded values of valley splitting 45 .
In conclusion, we argue for the atomic-level origin of valley splitting distributions in realistic Si/SiGe quantum dots, providing key insights on the inherent variability of Si/SiGe qubits and thereby solving a longstanding problem facing their scaling. We relate 3D atom-by-atom measurements of the heterointerfaces to the statistical electrical characterisation of devices, and ultimately to underlying theoretical models. We observe qualitative and quantitative agreement between simulated valley splitting distributions and measurements from several quantum dots, supporting our theoretical framework. Crucially, we learn that atomic concentration fluctuations of the 28 Si → Si-Ge heterointerface are enough to account for the valley splitting spread and that these fluctuations are largest when the envelope of the wavefunction overlaps with more Ge atoms. Moreover, while we have only incorporated random alloy disorder into our theoretical framework so far, we foresee that APT datasets including correlated disorder, such as steps, will be used to further refine our theoretical understanding of valley splitting statistics. Since atomic concentration fluctuations are always present in Si/SiGe devices due to the intrinsic random nature of the SiGe alloy, we propose to boost these fluctuations to achieve on average large valley splittings in realistic silicon quantum dots, as required for scaling the size of quantum processors.
Our proposed approaches are counter-intuitive yet very pragmatic. The interface broadening approach seems viable for hybrid qubits, which require valley splitting to be large enough to be usable but not so large as to be inaccessible. For single-electron spin qubits, which don't use the valley degree of freedom, the direct introduction of Ge in the quantum well appears better suited for targeting the largest possible valley splitting. By adding Ge to the Si quantum well in small concentrations we expect to achieve on average valley splitting in excess of 100 μeV. Early calculations from scattering theories 46 suggest that the added scattering from random alloy disorder will not be the limiting factor for mobility in current 28 Si/SiGe heterostructures. However, an approximate two-fold reduction in electron mobility was recently reported when an oscillating Ge concentration of about 5% on average is incorporated in the Si quantum well 45 . We speculate that fine-tuning of the Ge concentration in the quantum well will be required for enhancing the average valley splitting while not compromising the low-disorder potential environment, which is important for scaling to large qubit systems. We believe that our results will inspire a new generation of Si/SiGe material stacks that rely on atomicscale randomness of the SiGe as a new dimension for the heterostructure design.

Si/SiGe heterostructure growth
The 28 Si/SiGe heterostructures are grown on a 100-mm n-type Si(001) substrate using an Epsilon 2000 (ASMI) reduced pressure chemical vapor deposition reactor equipped with a 28 SiH 4 gas cylinder (1% dilution in H 2 ) for the growth of isotopically enriched 28 Si. The 28 SiH 4 gas was obtained by reducing 28 SiF 4 with a residual 29 Si concentration of 0.08% 47 . Starting from the Si substrate, the layer sequence for quantum well A comprises a 900 nm layer of Si 1−x Ge x graded linearly from x = 0 to 0.3, followed by a 300 nm Si 0.7 Ge 0.3 strain-relaxed buffer, an 8 nm tensily strained 28 Si quantum well, a 30 nm Si 0.7 Ge 0.3 barrier, and a sacrificial Si cap. The layer sequence for quantum well B comprises a 1.4 μm stepgraded Si (1−x) Ge x layer with a final Ge concentration of x = 0.3 achieved in four grading steps (x = 0.07, 0.14, 0.21, and 0.3), followed by a 0.45 μm Si 0.7 Ge 0.3 strain-relaxed buffer, an 8 nm tensily strained 28 Si quantum well, a 30 nm Si 0.7 Ge 0.3 barrier, and a sacrificial Si cap. In quantum well A, the Si 0.7 Ge 0.3 strain-relaxed buffer and the Si quantum well are grown at 750°C without growth interruption. In quantum well B the Si 0.7 Ge 0.3 strainrelaxed buffer below the quantum well is grown at a temperature of 625°C, followed by growth interruption and quantum well growth at 750°C. This modified temperature profile yields a sharper bottom interface for quantum well B as compared to quantum well A.

Atom probe tomography
Samples for APT were prepared in a FEI Helios Nanolab 660 dual-beam scanning electron microscope using a gallium focused ion beam at 30, 16, and 5 kV and using a procedure described in detail in ref. 48. Before preparation, a 150-200 nm thick chromium capping layer was deposited on the sample via thermal evaporation to minimize the implantation of gallium ions into the sample. All APT analyses were started inside this chromium cap with the stack fully intact underneath. APT was carried out using a LEAP 5000XS tool from Cameca. The system is equipped with a laser to generate picosecond pulses at a wavelength of 355 nm. For the analysis, all samples were cooled to a temperature of 25 K. The experimental data are collected at a laser pulse rate of 200-500 kHz at a laser power of 8-10 pJ. APT data are reconstructed using IVAS 3.8.5a34 software and visualized using the AtomBlend addon to Blender 2.79b and Blender 2.92 software. For the Voronoi tessellation the reconstructed data sets were exported to Python 3.9.2 and then tessellated using the scipy.spatial.Voronoi class of SciPy 1.6.2. Note that in these analyses the interfaces are represented as an array of sigmoid functions generated perpendicular to the respective interface on 3 nm × 3 nm tiles that are 1 nm apart. This sacrifices lateral resolution to allow for statistical sampling of the elemental concentrations but preserves the atomic resolution along the depth axis that APT is known to provide upon constructing the interface as shown in Fig. 2a.

Device fabrication
The fabrication process for Hall-bar shaped heterostructure field effect transistors (H-FETs) involves: reactive ion etching of mesa-trench to isolate the two-dimensional electron gas (2DEG); P-ion implantation and activation by rapid thermal annealing at 700°C; atomic layer deposition of a 10-nm-thick Al 2 O 3 gate oxide; deposition of thick dielectric pads to protect gate oxide during subsequent wire bonding step; sputtering of Al gates; electron beam evaporation of Ti:Pt to create ohmic contacts to the 2DEG via doped areas. All patterning is done by optical lithography. Quantum dot devices are fabricated on wafer coupons from the same H-FET fabrication run and share the process steps listed above. Double-quantum dot devices feature a single layer gate metallization and further require electron beam lithography, evaporation of Al (27 nm)

Electrical characterisation of devices
Hall-bar measurements are performed in a Leiden cryogenic dilution refrigerator with a mixing chamber base temperature T MC = 50 mK 50 . We apply a source-drain bias of 100 μV and measure the source-drain current I SD , the longitudinal voltage V xx , and the transverse Hall voltage V xy as function of the top gate voltage V g and the external perpendicular magnetic field B. From here we calculate the longitudinal resistivity ρ xx and transverse Hall resistivity ρ xy . The Hall electron density n is obtained from the linear relationship ρ xy = B/en at low magnetic fields. The carrier mobility μ is extracted from the relationship σ xx = neμ, where e is the electron charge. The percolation density n p is extracted by fitting the longitudinal conductivity σ xx to the relation σ xx / ðn À n p Þ 1:31 . Here σ xx is obtained via tensor inversion of ρ xx at B = 0. Quantum dot measurements are performed in Oxford and Leiden cryogenic refrigerators with base temperatures ranging from 10 to 50 mK. Quantum dot devices are operated in the few-electron regime. Further details of the 2DEG and quantum dot measurements are provided in the Supplementary Note 1.