Quantum simulation of the Hubbard model with dopant atoms in silicon

In quantum simulation, many-body phenomena are probed in controllable quantum systems. Recently, simulation of Bose–Hubbard Hamiltonians using cold atoms revealed previously hidden local correlations. However, fermionic many-body Hubbard phenomena such as unconventional superconductivity and spin liquids are more difficult to simulate using cold atoms. To date the required single-site measurements and cooling remain problematic, while only ensemble measurements have been achieved. Here we simulate a two-site Hubbard Hamiltonian at low effective temperatures with single-site resolution using subsurface dopants in silicon. We measure quasi-particle tunnelling maps of spin-resolved states with atomic resolution, finding interference processes from which the entanglement entropy and Hubbard interactions are quantified. Entanglement, determined by spin and orbital degrees of freedom, increases with increasing valence bond length. We find separation-tunable Hubbard interaction strengths that are suitable for simulating strongly correlated phenomena in larger arrays of dopants, establishing dopants as a platform for quantum simulation of the Hubbard model.

In this section we discuss the extraction of the energies of the states observed in dI/dU tunneling spectra. We find that the the number of observed states, the state's relative energies to eachother, and their energies relative to the sample's Fermi energy, are only compatible with two interacting holes on pairs of acceptors. We further estimate the total tunnel coupling Γ Σ to the tip and sample reservoirs for sequential hole transport, based on the dI/dU lineshape, and show that hΓ Σ is much less than the inter-acceptor tunnel coupling t for separations d considered in the main text.
For the acceptor pair with d/a B = 2.2 ( Figure 3A, main text), measured dI/dU is plotted in linear scale (Supplementary Figure 1A) and log scale (Supplementary Figure 1B). Away from the acceptor pair (Supplementary Figure 1A/B, orange squares), direct tunneling of electrons from the valence band to the tip (holes from the tip to the valence band) is obtained for U 0 V as illustrated in Supplementary Figure 1C, while direct tunneling of electrons from the tip to the conduction band occurs for U 1.15 V as illustrated in Supplementary Figure 1D. Above the pair (Supplementary Figure 1A/B, black squares), we observe a well-isolated peak at U ≈ 0.3 V, a collection of two strong (closely spaced) peaks at 0.6 V and 0.75 V, and fourth peak at U ≈ 1.0 V.
Plotted on a logarithmic scale, both the current (not shown) and numerically differentiated dI/dU data (Supplementary Figure 1B), black squares) grow exponentially with increasing U for the first peak before reaching a local maximum. As suggested by this observation 2,3 , we fit the first peak of the differential conductance to a state probed by sequential tunneling with broadening determined by the T = 4.2 K temperature of the reservoir, cosh −2 (eα(U − U 1 )/2k B T ), where e = 1.602 × 10 −19 C is the elementary charge, k B = 1.381 × 10 −23 J/K is the Boltzmann constant, U 1 is the peak voltage, and α is the well-known lever arm from sequential transport 2 describing the variation in the energy level ∆E = eα∆U due to a change in applied bias ∆U , as schematically illustrated in Supplementary Figure 1E.
The least-squares fit is in good agreement with the measured dI/dU for both the full-width at half maximum ∆U = 3.6k B T /eα and the exponential decay of the tail exp(eα(U − U 1 )/k B T ). A peak voltage U 1 = 334.0 ± 0.9 meV and a lever arm α = 0.0144 ± 0.0006 were obtained for tunneling of holes from the reservoir, to the acceptor pair, to the tip. The value for α is similar to those found for single-hole transport through individual acceptors in flashed p-type silicon 3 and single arsenic donors in flashed n-type silicon 4 . This excellent agreement shows that broadening of the single-hole transport peaks is due to the temperature associated with the Fermi-Dirac distribution in the hole reservoir.
The small value of α means that the "broadened" lineshape in Supplementary Figure 1A and Supplementary Figure 1B is due to weak electrostatic control of the sample potential by the STM tip, as expected due to strong screening by the carrier reservoir. This is further confirmed by the negligible "spectral shift" 5 , that is, the similarity of the apparent gap ≈ 1.15 eV in dI/dU measurements of Supplementary Figure 1A/B to the actual 1.14 eV band gap for heavily-doped p-type silicon at low temperature 6 .

A. Energy spectrum
We extracted the lever arm α and energies of states E i − E 1 = eα(U i − U 1 ), by least-squares fitting of dI/dU to a sum of single-hole tunneling peaks 1,3,4,7,8 . Fits Figure 1A and 1B (solid green lines) along with individual peaks (dashed green lines), on linear and log scales, respectively. Extracted energies E i − E 1 = eα(U i − U 1 ) of the excited states are plotted in Supplementary Figure 2, relative to the energy E 1 of the two-hole ground state, for d/a B = 2.2, d/a B = 2.7 and d/a B = 3.5.
For all d the extracted energy splittings between the states (Supplementary Figure 2) are too small to be associated with single-hole occupied states (A − 2 states) of coupled acceptors. For d/a B = 2.2, two energetically nearby states were found at 5 and 6.5 meV, while a higher energy state was found at ≈ 3.5 meV above those. Here, the observed splittings are incompatible with the "1s" hybridization energy 2t ≈ 25 meV between tunnelhybridized single-hole states estimated by numerics (Supplementary Note 3). Similar arguments hold for the d/a B = 2.7 (d/a B = 3.5), where the excited state splittings ≈ 2 meV (≈ 1 meV) are incompatible with the single-hole coupling energy scale of 2t ≈ 14 meV (2t ≈ 7 meV).
Even more importantly, it follows from the electrostatic arguments given below that the coupled-acceptor peaks observed in experiments are too close to the hole reservoir's chemical potential to correspond to ionizations of the very "deep" single-hole (A −1 2 ) states into zero-hole (A −2 2 ) states. Rather, they are only compatible with ionizations of the two-hole (A 0 2 ) ground and excited states of acceptor pairs into one-hole (A −1 2 ) states, which are energetically very similar to the ionization energy of a neutral (A 0 ) acceptor. For d/a B = 2.2 → 3.5 we found ionization transitions A −1 2 →A −2 2 for the last hole (Supplementary Note 3) require more than 100 meV energy. This exceeds the Fermi energy E F − E V ≈ 50 meV of the reservoir (E V is the valence band edge) by more than 50 meV. Now, for U = 0 V, the tip's contact potential bends states by an amount −eαU FB = −eα(W − Φ tip ) relative to the sample reservoir, where W = 5.2 eV is the work function of degenerately doped p-type reservoir, and Φ tip is tungsten tip workfunction. Let us assume a work function Φ tip = 4.8 ± 0.3 eV typical of a tungsten tip 3 . Then, the states bend down by 5.8±4.2 meV at zero bias, much smaller than the value > 50 meV required to depopulate a single hole transition). Therefore at zero bias, the observed acceptor pairs are in an A − 2 charge state, and a bias U = α −1 × 50 mV much less than −1 V would be required to remove the last hole by tip-induced band bending. This is incompatible with U ∼ 0.2 to 0.3 V of peaks observed in experiments. On the other hand, ionization energies of the two-hole states, A 0 2 →A − 2 of coupled acceptors are calculated (Supplementary Note 3) to be very energetically similar to A 0 ionization energies. Consequently, transitions populating ground and excited two-hole will occur for U > 0, explaining how they (much like single-acceptor ionizing resonances 3 ) are observed in our experiments.

B. Transport coupling to reservoir
We have used the dI/dU lineshape to obtain an upper bound on the tunnel coupling Γ Σ = Γ in + Γ out to the environmental reservoirs -the largest lifetime broadening energy scale 2 hΓ Σ that can be included in the lineshape analysis of the data 3 . As necessary to observe molecular states, we find Γ Σ for data in Figure 3A, 3B, and 3C, to be considerably smaller than the inter-acceptor tunnel coupling t estimated from the hybridization energies (Table 2). Tip height-dependent measurements (Supplementary Note 5) confirmed the exponential dependence of the tunnel current on tip height, indicating the rate to the tip is the slower (dominant) rate. The dwell-time (and coupling Γ Σ ), determined by the fast rate, is therefore controlled by the coupling to the hole reservoir in the substrate, Γ in in Supplementary Figure 1E.

Supplementary Note 2. ACCEPTOR DEPTH ESTIMATE
We have estimated the depth below the silicon surface of the acceptor dopant pairs, presented in Figure 3A, 3B and 3C of the main text, to be approximately 0.5 nm, 0.9 nm, and 0.9 nm, respectively. The depth of isolated dopants in scanning tunneling spectroscopy is typically estimated fitting the spatial variation of a spectral feature such as a band edge, in a bias condition where the dopant is ionized 3,4,9,10 , to a dielectric screened Coulomb potential. However, as illustrated in Supplementary Figure  1E, tunneling from the reservoir to the acceptor pair is a transition from a one-hole state to a two-hole state. Consequently, for voltages U below the lowest voltage peak at U = U 1 in Supplementary Figure 1A and Supplementary Figure 1B, the two-acceptor system is in Coulomb blockade with single-hole occupation, rather than being doubly ionized, such that fitting to dielectric-screened singleion potentials is inappropriate 3,4,9,10 . For U above U 1 the two-acceptor system has a high probability of occupation by two holes since Γ out Γ in . Hence, the method of fitting the band profile to dielectric screened ion potentials used in the above references is problematic.
To estimate the depth of our acceptors, we use the wellknown result that the spatial extent of the bound state measured by STS increases with increasing dopant depth (see e.g., reference 11). We compare the full-width at half-maximum W of the spatial tunneling distribution fit in Supplementary Note 4 to the same quantity calculated for subsurface acceptor-bound holes in the 6 × 6 spinorbit coupled valence band. We provide this estimate for completeness, since the fitting procedure described in Supplementary Note 4 does not rely on prior knowledge of the depth or Bohr radius of the dopant.
The full-width at half maximum W for the 1s orbital density φ s (r) ∼ exp(−|r − r 0 |/a * ) at z = 0 can be di- Evaluating this quantity we obtain W = 1.48 ± 0.12 nm, W = 2.25 ± 0.39 nm, and W = 2.19 ± 0.25 nm for the measured molecules with d/a B = 2.2, 2.7, and 3.5, respectively. Comparing these results to the full-width at half maximum W predicted for the orbitals using the numerical solution to Kohn-Luttinger 6 × 6 Hamiltonian presented in Supplementary Figure 3, the W extracted from experiments translate to depth estimates of ≈ 0.53 nm, ≈ 0.9 nm, and ≈ 0.9 nm, respectively.

A. Correlation with lever arms
When the acceptors are located deeper in the sample, the lever arms are anticipated to be smaller 3 . The relative depths of the planes in which the acceptor pairs are estimated to be located agree with the extracted lever arms. We obtain a lever arm α = 0.0144 ± 0.0006 for the shallower d/a B = 2.2 pair. The other two d/a B = 2.7 and d/a B = 3.5 pairs which are slightly deeper have lever arms α = 0.0088 ± 0.0010 and α = 0.0055 ± 0.0003, respectively. We note that fluctuations in the depth of the reservoir can also influence the lever arm. We anticipate that varying the dopant depths has a small effect on correlations in our experiments where the overall system is neutral. We have observed that the neutral level of nearsurface acceptors (N = 1) remains bulk-like 3 for depths 0.5-2.0 nm, as expected, due to a competition between the confinement and the dielectric mismatch. Fixed ionization energy essentially fixes the Bohr radius, fixing t and U/t. We note however that the apparent spatial extent of the wavefunction increases with increasing depth. This is not because of a change in the physical Bohr radius. Rather, it occurs because STM probes the wavefunction near the surface.

Supplementary Note 3. THEORY OF INTERACTING ACCEPTORS
We employed a microscopic configuration interaction framework to calculate eigenstates of interacting holes bound to proximate acceptors, in a Kohn-Luttinger 6 × 6 spin-orbit coupled representation (heavy holes, light holes, and split-off holes), and in a simplified single-band representation for reference. In both representations, the two-hole ground state at relatively small d/a B ∼ 2 was found to be predominantly composed of the |e + e− singlet of two even orbitals, where + and − denote +"3/2" and −"3/2" respectively for the valence band holes, and denote ↑ and ↓ respectively in the single-band approximation. The probability amplitude of the |o + o− singlet of odd orbitals was found to increase with increasing d, signaling interaction-driven generation of quantum correlations. The corresponding spatial tunneling probability distribution was predicted using the quasi-particle wavefunction (QPWF) density appropriate for few-particle states with quantum correlations 12,13 , and are presented herein for the non-trivial 6 × 6 spin-orbit coupled case.

A. Overview
In the absence of spin-orbit coupling, wavefunctions of interacting spins can be classified into singlets and triplets 14 . For holes in the valence band, the atomic orbitals have p-like symmetry, and the angular momentum L = 1 couples to the spin (intrinsic) angular momentum S = 1/2, such that the resulting Bloch states have definite J and J z . While intrinsic angular momentum is not a good quantum number, spin can nevertheless be generalized to hole pseudospin, i.e., Kramers doublets of time reversal symmetric linear combinations of Bloch states 15 .
Coupling of hole spins is well understood in the Kohn-Luttinger framework [15][16][17] . It has been theoretically shown that coupled hole pseudospin interact to produce singlets and triplets 15 , as for coupled electron spins. In the Kohn-Luttinger approach generalized valence band holes spins as written as spinors are the envelope functions corresponding to the Bloch states |J, m J . Due to spin-orbit coupling such spinors are ∼ 90 % polarized (rather than 100 % polarized) into single Bloch components. This produces a small band-mixing effect when the pseudospins are coupled together 16,17 , that we take into account exactly in our numerics.

B. Full-Configuration Interaction
We employ the Kohn-Luttinger FCI approach developed for valence band holes in reference [17]. Twohole eigenstates |Ψ are expressed as superpositions of antisymmetrized two-hole Slater determinants |αβ = c † α c † β |0 = 2 −1/2 (|α 1 |β 2 − |β 1 |α 2 ) with probability amplitudes d αβ , In |αβ above, c † j creates a single-hole state |j = c † j |0 , where |0 is the vacuum state. As discussed in detail below, the single-particle states |j were chosen as the eigenstates of the non-interacting Hamiltonian H 0 of a pair of acceptor ion potentials and a hard-wall interface with dielectric mismatch. The Hamiltonian for the two-hole system was taken as H = H 0 (r 1 ) + H 0 (r 2 ) + V h−h (r 1 , r 2 ) where V h−h (r 1 , r 2 ) is the Coulomb repulsion term for the holes. Substituting |Ψ above into H|Ψ = E|Ψ and taking the inner product with νλ|, we obtain Γ νλ,αβ = dr 1 dr 2 for each |νλ , where Γ νλ,αβ = ν| 1 λ| 2 V h−h |α 1 |β 2 is an off-diagonal (interaction) term evaluated in the spinor representation, J ν,λ = Γ νλ,νλ is the direct Coulomb interaction, and K νλ = Γ νλ,λν is a direct exchange interaction. Finite Γ νλ,αβ terms can introduce correlations 17 by creating eigenstates that are admixtures of two-particle Slater determinants |αβ . The single particle Hamiltonian employed was H 0 = H(k) + V ion,A + V ion,B + V image,A + V image,B + V wall taking into account ion potentials V ion,A/B of both acceptors, image charges V image,A/B due to dielectric mismatch 18 , and an infinite (hard-wall) V wall potential a distance 1a B ∼ 1 nm away from the ion representing the semiconductor/vacuum interface. For H(k) we considered both the 6×6 spin-orbit coupled Kohn-Luttinger valence band 19,20 H KL (k), as well as a single-band, isotropic, parabolic dispersion relation H p (k) = 2 k 2 /2m * where m * is an effective mass chosen to reproduce boron's binding energy 21 . The latter allows us to calculate states of a simple "scaled hydrogen molecule" without complex spin-orbit coupling and band structure effects. The ion potentials were taken as dielectric-screened Coulomb potentials. For the 6 × 6 scheme we included additional charge δq = −0.1e at the ion site as a crude core-correction to obtain the correct spatial extent r and ionization energy of the states 22 . The interaction term V h−h (r 1 s, r 2 ) is a dielectric screened Coulomb interaction including an image charge distribution produced by dielectric mismatch 23 .
For the single band theory, the four lowest singleelectron eigenstates of H 0 were used to construct 6 possible two-particle configurations, all of which were included for a full configuration interaction (FCI) approach. We considered two even, spin-degenerate orbitals |eσ = c † eσ |0 , and two odd, spin-degenerate orbitals |oσ = c † oσ |0 . For the 6×6 Kohn-Luttinger calculations, we use a basis that reflects the four low-energy states of single acceptors with s-like envelopes -the predominantly |m J | = 3/2 ("3/2") state, and the predominantly |m J | = 1/2 ("1/2") state 24 . Near an interface, the "1/2" state experiences more confinement than the "3/2" state, making the "3/2" state the ground state 1,25,26 . For lone subsurface acceptors ∼ 1 nm from an Si[001]:H interface, we have typically measured ∼ 1 − 2 meV splittings (see e.g., lone acceptor data in Supplementary Figure 2), in good agreement with few meV values predicted by 6 × 6 Kohn-Luttinger. See reference [1] for further details. In the two-acceptor potential H 0 these hybridize into eight single-hole eigenstates, which were used to construct 28 possible two-hole configurations, all of which were included in the state vector for a FCI approach. For the acceptor distances in the main text, four were predominantly composed of heavy holes ("3/2" states) and four were predominantly composed of light holes ("1/2" states). In both groups of four states, two are found within a Kramers doublet with even parity symmetry for the majority spin, and two are found within a Kramers doublet having odd parity symmetry for the majority spin 16,17 .
A numerical representation of the FCI matrix (Equation 2) was directly diagonalized to obtain state vectors d and energies E. The FCI matrix was obtained by evaluating the Coulomb matrix elements J νλ , exchange matrix elements K νλ , and the interaction terms Γ νλ,αβ by Monte-Carlo integration, using the VEGAS method for adaptive sampling and 10 7 iterations per integral. The basis of single particle eigenstates required for the integrals was obtained numerically using the finite difference scheme, for both the single-band and 6 × 6 representations, on a real-space grid of 110 × 110 × 110 points with a discretization step size of 0.21 nm, using Luttinger parameters for silicon's valence band 22 . To normalize the measured and theoretical (6 × 6 and single-band) interacceptor distances, we choose the value a B = 1.3 nm reproducing the 44.4 meV binding energy of a Boron acceptor in silicon 21 in a single-band approximation.

C. Correlations in ground state
In this section, the ground state of the FCI Hamiltonian in Supplementary Note 3 B is discussed. Our FCI predicts a singlet ground state for both the spinorbit coupled 6 × 6 representation and for the isotropic single-band case. In a single-band representation (ignoring spin-orbit coupling), the ground state and contributions for which for majority pseudospin have odd parity, For small d/a B ∼ 2, the dominant configuration (approx. 90 %) in the ground state singlet was found to be d (e,3/2)(e,−3/2) , while for increasing d, the probability amplitude d In this section, we discuss theoretical predictions of the correlation coefficient C in Figure 4A, as well as theoretical calculation of the quantity measured for spatially resolved single-hole tunneling into multi-particle states with correlations, the QPWF density 12,13,29 . The results of Supplementary Note 3 C are considered, for both the single-band and spin-orbit coupled 6 × 6 Kohn-Luttinger representations.
As described in Supplementary Note 1, the tunneling probability density measured in experiments is determined by tunneling of a hole from the acceptor molecule to the tip, leaving behind a single-hole eigenstate |i on the acceptor molecule. The corresponding QPWF density Γ i out (r) ∝ |M iS (r)| 2 is obtained 12 , where M iS (r) = i|Ψ(r)|Ψ S , Ψ(r) = j c j φ j (r) is the field operator, r = (x, y, z 0 ) is position of the tip apex orbital, and |Ψ S is the ground state. For the single-hole transport processes probed, it is necessary to sum over the final (single-hole) states |i of the molecule after hole tunneling to the tip. The total tunneling probability density is then Γ out (r) = i Γ i out (r). For the single band case, evaluating the QPWF using the field operator gives Γ out (r) ∝ |γ ee | 2 |φ e (r)| 2 + |γ oo | 2 |φ o (r)| 2 . Results for C = 2|γ oo | 2 are plotted in the main text ( Figure 4A, dashed line). For the 6 × 6 spin-orbit coupled representation, the quasi-particle wavefunction densities are also readily obtained, such that Γ out (r) Evaluating this expression results in a QPWF density that is a sum of even and odd QPWF contributions for all d and both orientations 100 and 110 , as expected considering the generalized parity of holes 16 . Results for C, again defined as twice the probability density of odd orbitals, are plotted in the main text ( Figure 4A, coloured lines) for both 100 and 110 orientations.
The predicted spatial tunneling probability density is shown 1 nm above the ion in Supplementary Figure 4 for an acceptor molecule with d = 3.8a B . The predicted correlation coefficient is C = 0.6, similar to C = 0.78 ± 0.16 extracted from measurements of the d = 3.5a B acceptor molecule, in Figure 4A of the main text. Contributions to the spatial maps (line profiles) from even and odd QPWF densities are separately plotted in Supplementary Figure  4A (Supplementary Figure 4D) and Supplementary Figure 4B (Supplementary Figure 4E), respectively. Their sum, the total tunneling probability density (line profile), is plotted in Supplementary Figure 4C (Supplementary Figure 4F). Upper and lower grey lines in Supplementary Figure 4F denote normalized QPWF densities for the even QPWF and odd QPWF components, respectively, which in this case directly correspond to the profiles in Supplementary Figure 4D and 4E, respectively.

Supplementary Note 4. FITTING OF MEASURED SPATIAL TUNNELLING PROBABILITIES
The experimentally obtained spatial tunneling probabilities for the two-hole ground state of acceptors, in Figure 3A, Figure 3B, and Figure 3C of the main text, were fit assuming even and odd linear combinations of atomic orbitals with s-like envelopes for φ e (r) and φ o (r). This data was obtained using a tip height established by constant current imaging conditions, at a bias where the topography was nominally flat apart from dimer row corrugations.
In this section we describe the fits and show that the probability density of an s-like wavefunction describes the measured probability density of the ground state and first excited state of an isolated acceptor, in agreement with what is expected for the s-like "3/2" pseudospin ground state and s-like "1/2" pseudospin excited state.
Spatially resolved tunneling spectra were measured for isolated acceptors, consistently revealing two dI/dU peaks associated with the s-like J = 3/2 ground state manifold. The ∼ 2 meV splitting of the isolated acceptor in Supplementary Figure 2 was found to be in good agreement with single-particle eigenstates calculated from our numerical Kohn Luttinger solver for depths determined as in previous work 3 , that is, by a fitting the profile to the ionized acceptor's perturbation to the band edge 9,10 . The spatial tunneling probability distributions for the lowest energy state and first excited state of an isolated acceptor were obtained by integrating the two low energy dI/dU peaks. Results are shown in Supplementary Figure 5A and Supplementary Figure 5B for the single acceptor in Supplementary Figure 2.
The dimer lattice structure of the hydrogen-terminated 2 × 1 surface is evident in the states Supplementary Figure 5A and Supplementary Figure 5B, running along the [110] direction indicated. The dimer and lattice frequencies were Fourier filtered out of the image 30 , after which only the characteristic envelope of the probability density of the states remain. Results for the ground state and first excited state are presented in Supplementary Figure  5C and 5D, respectively. The tunneling probability density reflects the hole density, and the appearance of two s-like envelopes is in agreement with expectations of the s-like envelopes for the "3/2" and "1/2" states 1, 24 . We observe that the measured low-frequency envelope for the probability density is closer to isotropic for the ground state, and spatially elongated along the dimer direction for the excited state.
The probability density was fit along a line in the plane (z = 0) of the surface of the form |φ s (r)| 2 , where φ s (r) ∼ exp(−|r − r 0 |/a * ) is a 1s orbital envelope function. In this expression, r = (x, y, z) is the tip position, r 0 = (x 0 , y 0 , z 0 * ), x 0 and y 0 are centre coordinates of the orbital, z * 0 is an effective depth, and a * is an effective Bohr radius. The profile of the measured ground state (solid squares) and least-square fit (solid line) are shown for a profile along the direction parallel (black) and perpendicular (green) to the dimer row in Supplementary Figure 5E, demonstrating excellent agreement. The same quantities are plotted for the first excited state in Supplementary Figure 5F. A more noticeable anisotropy of the excited "1/2" state was found.
Motivated by the theoretical observation that the QPWF density partitions into contributions from even and odd states, we categorize the measured QPWF density into contributions from even and odd QPWFs, that is, linear combinations of atomic orbitals with even and odd parity, respectively. For a distance R between acceptors, we use φ e (r) = (1/I e (R))[φ s (r−R/2)+φ s (r+R/2)] and φ o (r) = (1/I o (R))[φ s (r − R/2) − φ s (r + R/2)], where I e (r) = 1 + (1 + r + r 2 /3) exp(−r) and I o (r) = 1 − (1 + r + r 2 /3) exp(−r) in atomic units. In the fits presented in the main text, a * , z * 0 , R, and |γ ee | 2 are free parameters, while 1 = |γ oo | 2 + |γ ee | 2 is simply a consequence of normalization. Effective Bohr radii a * ∼ 1 nm and depths z * 0 ∼ 1 nm were obtained from the least squares fits. In our fits the values for a * and z * 0 define a full-width at half maximum W for the spatial tunneling probability of each acceptor. In Supplementary Note 2, the actual depth z 0 of the acceptors was estimated using W derived from experimentally fit values a * , z * 0 .

A. Error analysis in fit
The sum of square errors was used to determine the importance of the Coulomb correlation parameter γ oo in the fits. The first column (Fig #1) of Supplementary Table 1 gives the results for the SSe when |γ oo is a free parameter, while the second column (Fig # 2) of Supplementary Table 1 gives the results for SSe when |γ oo is artificially forced to zero.
Here we see the fit is considerably improved by including the Coulomb correlation parameter in the fit.

B. Independence of fit results on tip height
Due to the very small lever arm α ∼ 10 −2 , the data is acquired very near flat-band 3 , and as such we do not observe the so-called "ionization parabolas" (observed elsewhere and associated with tip-induced band bending 9,31 ) for single B:Si acceptors 1 , or for B:Si acceptor pairs (Fig.  1b, main text).
To provide an independent check that the tip does not influence the results of the QPWF measurement, we present measurements on the d/a B = 2.2 and d/a B = 3.5 pairs from the main text at a tip height +60 pm above the position where the data was taken in the main text.  Supplementary Fig 7D), are in excellent agreement with the same QPWF model (solid red line, Supplementary  Fig 7C and Supplementary Fig 7D) used in the main text. Reference curves are shown for the totally uncorrelated state (γ oo = 0) and maximally correlated state (γ oo /γ ee = 1).
The extracted results for the correlations are |γ oo | 2 = 0.15 ± 0.06 for d/a B = 2.2 and |γ oo | 2 = 0.40 ± 0.09 for d/a B = 3.5. These results for +60 pm tip heights are within the experimental errors of the results for the data in the main text, which are |γ oo | 2 = 0.12±0.06 for d/a B = 2.2 and |γ oo | 2 = 0.39 ± 0.08 for d/a B = 3.5.

C. Large distance limit of fitting model
For any given experimental signal-to-noise ratio, there always exists a large enough d/a B where it is impossible to distinguish the existence of correlations in the ground state QPWF. This follows from the fact that our model determines correlations based on the contrast between orbitals |φ e (r)| 2 and |φ o (r)| 2 . However, since |φ e (r) are atomic orbitals at sites A and B, and the product 2φ A (r)φ B (r) is exponentially suppressed with increasing d/a B . The results of Table 1 show that the Coulomb correlations embodied by |γ oo | can be reliably detected, for the signal-to-noise ratio of our experiment.  Figure 6D) are slightly spatially elongated along the 110 direction. This is consistent with the expected QPWF density in the interpretation of Figure 5 (main text), that the excited state orbitals have one hole occupying a "1/2" pseudospin state and the other occupying a "3/2" pseudospin state. These states display slightly different spatial anisotropies (Supplementary Figure 5C and Supplementary Figure 5D).
While we have measured the spatial tunneling probabilities for the excited states, measuring excited-state quasi-particle wavefunctions requires slow tunnel-in from the tip into the dopant system (not slow tunnel-out to the tip). Slow tunnel-in from the tip is required to make the current limiting slow tunnel rates additive for excited states (rather than competitive with the ground state as in the case herein), in the presence of Coulomb blockade. Some details are given in reference 8.

Supplementary Note 5. TIP-HEIGHT DEPENDENT SPECTRA
In this section we discuss our confirmation of our model for the thermally broadened sequential tunneling in Supplementary Fig 1E, by detailed investigation of the tipheight dependence of the current, in agreement with our results on Arsenic donors in Si 4,8 . This is shown in Supplementary Figure 8A for tip heights that sequentially increase by 15 pm per curve, above a single B:Si acceptor presenting the unique topographical signature discussed elsewhere 1,3 . The data shows broadened current steps centred at U ≈ 0.25 V and U ≈ 0.5 V that have an exponential dependence on the tip height, and which fit very well the lineshape of purely thermally broadened current resonance of a dopant 1,3,4,7,8 . The peaks shown here correspond to the ±3/2 ground state and ±1/2 first excited state discussed in Supplementary Note 4 and reference 1.
Qualitatively, the exponential dependence on tip height verifies the assumption that Γ in Γ out , as expected for the large vacuum tunneling barrier, and as required to obtain the results in Figure 3 of the main text. We performed least-squares fitting to extract the height I 0 of the thermally broadened current steps. These I 0 are plotted as a function of tip height z − z 0 in Supplementary Figure 8b for the ±3/2 ground state (red curve) and the ±1/2 excited state (blue curve). The data are in excellent agreement with a z-dependence exp(−2κ(z − z 0 )) where κ = (1.15±0.05)×10 10 m −1 . This value is in good agreement with the expected value for κ = √ 2m 0 Φ B / for the ≈ 5 eV barrier expected for shallow valence band acceptor in silicon. Moreover, we find that the lever arm (not shown) and energy splitting (Supplementary Figure  8C) between the ground and excited state does is essentially independent of tip height. This provides further verification of the sequential tunneling model, and that the tip does not strongly influence the acceptor-bound states.
The dependence of γ ee and γ oo on U/t in Figure 2b can be found by a simple transformation. The orthonormal localized states of the Hubbard dimer model created by operators c † Aσ and c † Bσ are c † Aσ = w −1/2 (c † aσ − gc † bσ ) and (8) respectively 32 . Here c † aσ and c † bσ create (non-orthogonal) atomic orbitals with spin σ, w = 1 − 2Sg + g 2 , g = (1 − √ 1 − S 2 )/S, and S = a|b is the overlap between the atomic orbitals.
Re-writing the even and odd eigenstates of the singleparticle potential and equation the two singlet ground states we find γ ee = γ c + γ i and γ oo = γ c − γ i . These relationships are used to map the Hubbard model solution to the molecular orbital model in Figure 2B of the main text. This mapping is independent of normalization constants w, g, and S.

Supplementary Note 7. DEGREE OF ENTANGLEMENT
For fermions, entanglement is defined in terms of Slater decompositions 33 rather than the Schmidt decompositions of distinguishable particles, and can be expressed independent of basis by a degree of entanglement. Here we follow derivation in reference [28] for entanglement to describe correlations in a two-site Fermi-Hubbard system. For two particles that obey Fermi statistics, is characterized by an anti-symmetric matrix ω a,b , that is, a Slater decomposition. Transformed into a block diagonalized form diag[Z 0 , Z 1 , ..., Z N ] through a unitary rotation of the single particle states, where the number of nonzero z i is called the Slater rank, and if the Slater rank is 1, the quantum correlation of the state is zero 34 . It has been shown 28 that the z 2 i are the eigenvalues of the basis independent quantity ω † ω, and based on this observation, the degree of entanglement S of the two particles was defined as For a superposition of "even/even" and "odd/odd" sin-glets |S = (γ ee c † e,↑ c † e,↓ − γ oo c † o,↑ c † o,↓ )|0 , It directly follows that the non-zero eigenvalues z 2 i of ω † ω are |γ ee | 2 and |γ oo | 2 . Then, the degree of entanglement for the state |S is S = −|γ ee | 2 log 2 |γ ee | 2 − |γ oo | 2 log 2 |γ oo | 2 .