Experimental realization of an extended Fermi-Hubbard model using a 2D lattice of dopant-based quantum dots

The Hubbard model is an essential tool for understanding many-body physics in condensed matter systems. Artificial lattices of dopants in silicon are a promising method for the analog quantum simulation of extended Fermi-Hubbard Hamiltonians in the strong interaction regime. However, complex atom-based device fabrication requirements have meant emulating a tunable two-dimensional Fermi-Hubbard Hamiltonian in silicon has not been achieved. Here, we fabricate 3 × 3 arrays of single/few-dopant quantum dots with finite disorder and demonstrate tuning of the electron ensemble using gates and probe the many-body states using quantum transport measurements. By controlling the lattice constants, we tune the hopping amplitude and long-range interactions and observe the finite-size analogue of a transition from metallic to Mott insulating behavior. We simulate thermally activated hopping and Hubbard band formation using increased temperatures. As atomically precise fabrication continues to improve, these results enable a new class of engineered artificial lattices to simulate interactive fermionic models.


Supplementary Eq. (1)
The electrostatic charge at site i can be expressed as = − + ∑ , , where ∈ [0, 1,2] is the integer number of electrons added onto the site via hopping, and , is the capacitively induced charge when applying a voltage on electrode . Here, e is the absolute charge value of an electron; and we neglect the background induced charges from the substrate. The classical charge at the array can be expressed as, The electrostatic potential V and the total electrostatic energy P of the array can be obtained from, Here is a function of integer numbers of electrons at each site and the voltage values at each electrode . When = 0 for all sites, we denote the total electrostatic energy as 0 . At a given set of voltage conditions, we denote and as the electrostatic chemical potential and on-site charging energy to add the first and second electrons, respectively, onto site , which can be numerically calculated from, Note, it can be numerically verified that and , are independent of V l . On the other hand, and , depend critically on the size of the metallic disk that is used to represent a quantum dot in the classical capacitance model. In each of the array devices in this study, we chose the disk size so that the average value of (see Supplementary Table 1(a)) matches (within +/-2meV range) the average expected on-site charging energy values of few-dopant clusters from previous atomic tight-binding calculations (See the last column in Supplementary Table 3).
We denote α l,i as the lever arm of an electrode l to a lattice site i, which is defined as the ratio between ∆ , , the shift in the electrostatic chemical potential at site i in response to ∆ × , the change in the electrostatic potential at electrode l. , = ∆ , ∆ × Supplementary Eq. (6) As described in the main text, the gate-gate slope at a charge addition boundary corresponds to the gate-gate level arm ratio of the added electron, that is, the conducting electron, and characterizes whether the single electron is added onto sites in the upper, middle, or lower row in the array. Here we define the effective lever-arm of gate to an ( + 1) ℎ addition-electron, which, with finite hopping, has the charge distribution ∑ ∆ , , , as, Here ∆ represents the change in occupation at site when a single electron is added onto the array, altering the charge distribution of the array from an eigenstate of a total charge number to that of an eigenstate of + 1. The charge number conservation condition requires that ∑ ∆ = 1. On a numerically simulated charge stability diagram, we can extract the added electron occupations at charge addition boundaries, see plots in Fig. 3c in the main text, for instance, and calculate the lever arm ratio 1, → +1 2, → +1 at each pixel point. In the right panels in the main text Fig. 4, we plot the histogram distributions of the lever arm ratio calculated from the simulated charge stability diagrams for the first, second, and third arrays, respectively.
On a conductance map spanned by voltages of the two in-plane gates, the gate-gate lever arm ratio equals the slope of the conduction line 1 2 = −∆ 2 ∆ 1 . Indeed, we observe qualitative agreement between the measured slope histograms and calculated lever arm ratio histograms in the main text Fig. 4  At each quantum dot in the STM images of the arrayed device hydrogen lithography patterns (see Supplementary Fig. 2.), we identify single dangling bonds, valid and invalid dimer sites for P incorporation sites using red dots, green ellipses, and red ellipses, respectively. Previous studies [2][3][4] have shown that the number of incorporated P atoms at STM-patterned quantum dot can be estimated by counting the number of hydrogen-desorbed dimers that are valid for P incorporation. We have carried out previous studies to verify the dopant atom incorporation probability in patterned devices by comparing the STM lithography to few atom quantum dot spectroscopy as described in more detail in Supplementary Note 3. Having at least three adjacent H-desorbed dimers within the same dimer row is a necessary condition to incorporate one P dopant into the surface silicon lattice. We define the best estimate of the number of incorporated P atoms at each site by requiring three contiguous dimers which are then counted as capable of incorporating a single P dopant. The upper bound of our estimate of the number of incorporated P atoms is set by 25% of the number of dangling bonds sites within allowed dimers, according to the 0.25 monolayer of P incorporation density in saturatedoped Si:P -layers. 5 To estimate the lower bound of incorporated P atoms we follow Fuchsle et al., 4 where it is found that, with a saturation dose of phosphine gas at room temperature, the P incorporation density in nm scale desorbed areas decreases to 0.09 monolayer of contiguous desorbed areas within the quantum dot pattern due to competition for dangling bond sites to lose H atoms from absorbed PHx (x = 1, 2) during incorporation. These estimates agree with our independent evaluation of incorporation probabilities using few dopant atom transistors as described in Supplementary Note 3. We round the estimated lower and upper bound dopant numbers to their nearest integers. Supplementary Table 1 lists the best estimated numbers of dopants at each quantum dot and the lower and upper bounds of the estimates. While it is possible to count the number of incorporated dopants directly by imaging the quantum dots after P incorporation using STM, we avoid this time-consuming step to minimize the exposure of the reactive patterned device to contamination before epitaxial encapsulation as well as the risk of unintentional tipsurface interactions that could introduce atom-scale contaminations and defects.  1  2  3  1  3  4  2  3  5  D12  1  2  3  2  3  5  2  3  5  D13  0  2  4  1  2  3  2  3  5  D21  0  1  1  1  2  3  2  3  6  D22  1  2  3  1  3  4  2  3  6  D23  1  2  3  1  2  3  0  1  2  D31  1  2  3  2  3  5  2  3  5  D32  0  1  1  1  2  3  2  3  5  D33  1  2  4  1  2  3  1  2  3 Supplementary  6 E b values are with respect to the conduction band edge. The variations for 2P and 3P clusters represent one sigma (70%) in cluster configuration distributions. 6 We account for the variations in the best-estimate numbers of dopant atoms per site (Supplementary

Supplementary Note 3. Characterizing an individual few-dopant quantum dot
central region of the few-dopant quantum dot after hydrogen lithography, but before PH3 dosing. (b) Close-up STM image of the H-desorbed quantum dot region. The number of exposed Si dangling bonds (DB) and dimers can be counted by overlaying the Si (100) 2x1 surface reconstruction lattice grids with the STM images after hydrogen lithography. The allowed and forbidden P incorporation sites are highlighted in green and red respectively. (c) the low-temperature (T=4 K) differential conductance charge stability diagram. As indicated by arrows at the ↔ +1 and -1↔ transitions, there appears symmetric resonant tunneling features at positive and negative biases, indicating approximately equal tunnel coupling between the dot and the drain and source leads. The occupation number of the dot is expressed using an integer . (d) Following the method as described by Weber and co-workers, 6 we extract the binding energy spectrum of the few-dopant quantum dot from (c) and overlay the extracted binding energy levels (blue and red triangles) on the experimental and theoretical binding energy d spectrum previously published by Weber and co-workers. 6 The overlay indicates that there are 3 dopant atoms incorporated in the dot shown in (b).

Supplementary Note 4. Impact of hopping, interactions, and disorder on simulated addition energy spectrum and charge stability diagrams
Supplementary Fig. 4. Impact of hopping, interactions, and disorder on simulated addition energy spectrum and charge stability diagrams. (a) Simulated charge addition energy spectrum of the third array as a function of hopping amplitude ( ) with all leads and gate electrodes at zero ground potential. The impact of long-range electron-electron Coulomb repulsion ( , ) and long-range electron-ion core Coulomb attraction ( , ), on the charge addition energy spectrum are illustrated by comparing plots in the first column calculated using the full Hamiltonian (Equation 1 in the main text) and plots in the subsequent columns where , and/or , are turned off. The panels in the upper row assume identical three-dopant quantum dots with no variation in binding energy, while the lower row of panels takes into account the estimated disorder in binding energy (see Supplementary Tables 2 and 3). We have calculated with different sets of disorder configurations and found that they do not alter the qualitative features as shown in the lower panels in (a). The 18 addition energy levels correspond to a total of 18 excess electrons modeled in the 3x3 array system. (b) Simulated charge stability diagrams of the third array illustrating the impact of long-range and on-site electron-electron Coulomb interactions. Diagrams in the panels of the first column are calculated using the full Hamiltonian; diagrams in the second and the third column panels are calculated respectively by turning off the long-range electron-electron interaction ( , ) only and by turning off both the long-range and on-site electron-electron interactions ( , and ) all together. Following the simulation conditions in (a), diagrams in the upper row of panels and lower row of panels are calculated assuming an ideal and a disordered array configuration, respectively.
In Supplementary Fig. 4(a) 9 In this study, we extend previous efforts in simulating dopant arrays by including two experimentally defined in-plane gates in the extended Hubbard model and explore their impact on the charge distributions and addition energy spectrum in the array.
In Supplementary Fig. 4(b), we illustrate the impact of long-range and on-site electron-electron Coulomb repulsion ( , and ) on the simulated charge stability diagrams. As expected, the absence of electron-electron interactions reduces the charge stability regions and the overall span of the charge stability spectrum over the gate-gate space. Eliminating the many-body interaction effects (no , and no ) all together results in charge stability diagrams of the single-particle energy spectrum in the noninteracting regime. The dominate charge stability regions in the noninteracting regime are determined by the overall chemical potential landscape within the array, a combined effect of on-site binding energy, long-range electron-ion Coulomb attraction, and electrostatic potential from electrodes.

Supplementary Note 6. Impact of disorder on simulated charge stability diagrams for the third array
Supplementary Fig. 6. Understanding the impact of disorder and gate potential gradient on the ground eigenstate charge stability diagram. The simulated charge stability diagram in (a) is based on an idealized array with no disorder; all lattice sites are of identical three-dopant cluster quantum dot (binding energy E b = −81 meV) and all nearest-neighbor hopping amplitudes t = 0.5 meV. Charge occupation numbers at representative regions are highlighted in the diagram. When the two in-plane gates are tied and sweeping together along the diagonal direction, a finite analog Mott gap is observed at half-filling (occupation = 9), separating the lower Hubbard band from the upper Hubbard band. When a differential voltage is applied across the two in-plane gates and swept from zero to a large differential value, the gap at half-filling closes, and two other prominent addition-energy gaps at one-third (occupation = 6) and two-thirds (occupation = 12) filling opens. This process occurs due to the lowering of the chemical potential of the first or the third rows of lattice sites causing the electrons to energetically favor doubly occupying those rows and overcoming the on-site charging energy (U), as opposed to half-filling of the entire array. In (b), we introduce disorder in the number of dopant atoms per site based on the best estimates from experimental lithographic patterns shown in the main text Fig. 1c and Supplementary  Fig. 2: one dopant at D23 (E b = −47 meV), two dopants at D33 (E b = −70 meV), and three dopants for the rest of the dots (E b = −81 meV). Because of the higher energy cost to add the first electron onto D23, the central Mott-like gap along the diagonal direction now appears with charge number 8. In (c), the dopant number at each site is the same as those in (b), however, random noise from a Gaussian distribution with zero mean and standard deviation of 10 meV has been added to E b to gain a better idea about the effect of variations in the number of dopants at each site. Note that (c) is not averaged over disorder realizations; it represents a typical case in which there is significant variation in E b across the array. Along the diagonal direction, the added disorder further broadens the upper and lower Hubbard bands and suppresses the gap. that can be added to the 3x3 array at gate voltage conditions along the dashed line in (a). We emphasize the apparent similarity between the charge addition boundaries in the gate-gate charge stability map and the charge addition energy levels in the addition energy spectrum, manifesting the same charge addition response to the chemical potential gradient across the array in the gate-gate direction; the charge addition energy is flat versus the potential gradient when an electron is added to the middle row sites; the charge addition energy increases (decreases) with increased negative voltage on the upper gate when an electron is added to the upper (lower) row sites. (c) Charge occupation configurations on the left and right sides of selected avoided crossings in the addition energy spectrum region that are highlighted in (b). The charge configurations become resonant at the avoided crossing. As discussed in (b), it is energetically favorable for an electron to be added to different lattice sites on the left and right sides of an avoided crossing, as indicated by the different addition energy slopes on either side. The energy separation at an avoided crossing in the addition energy spectrum characterizes (U m + 2t ′ ), where U m is the mutual charging energy between the two involved charge addition distributions and t ′ is the resonant tunnel coupling between the two charge distribution configurations. The subplot in (c) illustrates that the U m contribution at an avoided crossing can be obtained by extrapolating the linear sections on both sides of the avoided crossing, in analogy to standard practices in double-dot systems.
Since the system is in the strong interaction regime (U i , U i,j ≫ t), U m dominates the avoided crossing separation. (d) Schematic illustration of eigenstate charge distributions and single charge addition at the left and right sides of the avoided crossings as highlighted in (c). Note, the zero addition energy level in (c) corresponds to the Fermi level in the source and drain leads. Taking charge addition #7 as an example, the gray scale plots represent the charge distributions of many-body ground states when there are 7 electrons in the array. The red-blue plots represent the change in charge distributions when adding the 7 th electron onto the array and the system changes from a 6-electron ground state to a 7electron ground state. At the avoided crossings when adding charges #2, #4, and #11, the involved charge addition sites are single lattice sites, and the avoided crossing separation is large (small) when the two resonant sites are close to (far apart from) each other. At the avoided crossing for adding charge #7, due to finite hopping amplitude, the addition of a single electron leads to changes in occupation at multiple lattice sites. The resonance between the two coupled many-body charge configurations is a manifestation of the complex many-body interactions, such as the mutual charging energy (long range Coulomb interactions) and tunnel coupling, between many-body states in the array system. Supplementary  Fig. 7 Supplementary Fig. 7 (c). In analogy to single-electron resonant tunneling in a double-dot system, the addition energy separation at an avoided crossing in a 3x3 array with two in-plane gates is characterized by ( + 2 ′ ), where is the mutual charging energy (long-distance e-e repulsion) between the two charge addition distributions (see the red-blue charge distribution plots in Supplementary Fig. 7(d)) that are on resonance at the avoided crossing, and ′ is the resonant tunnel coupling between them. The listed addition energies (2 nd column) at the avoided crossings are extracted from the numerically simulated addition energy spectrum with hopping t = 0.5 meV, as shown in Supplementary Fig. 7(c). The at an avoided crossing is illustrated in the subplot in Supplementary Fig. 7(c). The values (6 th column) at the avoided crossings are extracted from numerically simulated addition energy spectrum with hopping t = 0 meV. The uncertainty in the least significant digit is limited by the x-axis (gate voltage) resolution in the numerical simulation. The ′ values in the last column are estimated as the difference between the addition energy and at the avoided crossings. The numbers support the observation that the mutual charging energy is inversely proportional to the spatial distance between the charge occupation configurations, while the resonant tunneling rate is exponentially dependent on this spatial distance.

Avoidedcrossings in
region reproduced from Fig. 3a in the main text. (b) Simulated eigenenergy spectrum for occupation = 8 eigenstates along the dashed-line detuning axis in (a). (c) Close up eigenenergy spectrum at the red circle in (b). The eigenenergy lines are plotted with 80% transparency, and the color intensity reflects the level of degeneracy. The eigenenergy levels bundle into minibands of highly degenerate or neardegenerate eigenstates. (d) The charge distribution of selected eigenstates at different locations along the detuning axis as highlighted in (c). Charge distribution configurations are nearly identical for eigenstates within the same miniband. Moving along the detuning axis, the ground state alters its charge distribution. The energy separation at the avoided crossing shown in (c) is determined by the tunnel coupling t ′ between the two distributions whose eigenstates are hybridizing at the avoided crossing.

Supplementary Note 8. Characterizing addition energy spectrum in the second and third arrays
Supplementary Fig. 9. Extracting addition energy levels in the third array using the measured Coulomb blockade diagrams. (a) Coulomb blockade conductance diagram taken along the 45-degree diagonal axis that is marked by dashed lines in (b). Horizontal dashed lines in (a) correspond to the bias voltage level of (b). The vertical dashed lines in (a) mark charge addition positions on the gate-voltage axis where conductance via the addition electrons becomes visible at finite biases. (b) Conductance gate-gate maps that are measured at representative bias voltages, +10 mV, +4 mV, -10 mV, and -4 mV, respectively. Maps are plotted in log-color scale to show conductance features of small amplitudes. (c) Charge addition energy spectra extracted by scaling the measured charge addition levels using properly chosen gate lever arms. The slopes of the charge addition lines in the gate-gate maps allow us to determine whether the conducting electron is added onto sites in the upper, middle, or lower row of the array, and to choose proper gate lever arms for converting charge addition positions on the gate voltage axis to energy levels in the addition energy spectrum. N 0 represents the initial number of excess electrons in the array.
Experimental determination of the charge addition spectrum in a quantum dot array is a prerequisite for more complex measurements and control of many-body states in the array. We use the positions of Coulomb oscillation peaks in the measured Coulomb blockade diagram, as marked by vertical dotted lines in Supplementary Fig. 9 (a), to extract the addition-energy spectrum along the 45-degree diagonal axis in the gate-gate voltage parameter space. To choose the proper gate lever arms to convert measured charge addition levels from the gate voltage axis to energy space, we determine the average gate lever arm ( = (∑ 1, + 2, ∈ )/3 for adding an electron to sites in row ) according to the gate-gate slope of its charge addition boundary. We plot the charge addition energy spectra in Supplementary Fig. 9 (c), for both positive bias and negative bias conditions. The color of the energy levels corresponds to the gate-gate slope of the conduction lines, which also corresponds to the row of sites where the conducting electron is added. The differences in the extracted addition energy levels between a positive bias condition and a negative bias condition result from the differences in specific atomic configurations along the source/drain direction. We note that the absolute numbers of excess electrons in this device are unknown; however, the assigned charge numbers in Supplementary Fig. 9 (c) do not affect the underlying physical phenomena in this work. Supplementary Fig. 10. Characterizing the addition energy spectrum in the second array. (a) Gate-gate map of conductance through the array measured at V bias = 6 mV and at the base temperature of the dilution refrigerator (T = 10 mK). (b) (c) Finite bias differential conductance spectroscopy taken from the cut along the dashed line in (a), and measured at T = 10 mK (b) and T = 12 K (c), respectively. When doing a transport measurement, a large positive bias at the drain lead (with source lead grounded) pulls down the chemical potential in the drain lead, and individual electrons propagate through the lower Hubbard band within the bias window. Therefore, the differential conductance spectrum at positive (negative) bias in (b) represents the addition energy spectrum of the lower (upper) Hubbard band. The addition-energy levels in the measured lower Hubbard band section are detuned by the differential gate-gate voltages and exhibit crossings of addition energy levels (solid circle) in (b) and (c), which corresponds to the avoided crossing region (solid circle) in (a).
As discussed in the main text, avoided crossing separations in charge stability diagrams are determined by the effective mutual charging energies and tunnel coupling ′ between relevant many-body states. Comparing the simulated diagrams for the second array (left panels in Supplementary Figs. 5a and 5b) with those for the third array (main text Figs. 3a and 3b) we observe the expected larger avoided crossing separations.