Adequacy of Si:P chains as Fermi–Hubbard simulators

The challenge of simulating many-body models with analogue physical systems requires both experimental precision and very low operational temperatures. Atomically precise placement of dopants in Si permits the construction of nanowires by design. We investigate the suitability of these interacting electron systems as simulators of a fermionic extended Hubbard model on demand. We describe the single-particle wavefunctions as a linear combination of dopant orbitals (LCDO). The electronic states are calculated within configuration interaction (CI). Due to the peculiar oscillatory behavior of each basis orbital, properties of these chains are strongly affected by the interdonor distance R0, in a non-monotonic way. Ground state (T = 0 K) properties such as charge and spin correlations are shown to remain robust under temperatures up to 4 K for specific values of R0. The robustness of the model against disorder is also tested, allowing some fluctuation of the placement site around the target position. We suggest that finite donor chains in Si may serve as an analog simulator for strongly correlated model Hamiltonians. This simulator is, in many ways, complementary to those based on cold atoms in optical lattices—the trade-off between the tunability achievable in the latter and the survival of correlation at higher operation temperatures for the former suggests that both technologies are applicable for different regimes. Precisely positioned impurities in silicon could perform quantum simulations of models that are intractable for classical computers. Amintor Dusko and co-workers from Brazil’s Universidade Federal do Rio de Janeiro and the University of Ottawa in Canada have shown numerically that if phosphorus atoms in silicon devices are arranged in a line, their electrons can provide an experimental realization of the many-body Fermi-Hubbard model. This model, first presented in 1963, provides a heavily simplified description of high-temperature superconductors. However, a theoretical solution cannot be obtained for large systems. By simulating the model experimentally, Dusko et al.’s proposed system circumvents these computational limitations, and could reach temperature regimes that are not accessible with state-of-the-art cold atom simulators. Successful experimental realization could therefore provide new insights into a model that has been intensively studied since the 1960s.


INTRODUCTION
Strongly interacting fermions are the main ingredient for some phenomena in the forefront of physics, such as high-T c superconductivity and topological phase transitions. [1][2][3] In one dimension, correlations may be identified through collective electronic behavior, such as charge density wave (CDW), bond-order wave (BOW) and spin density wave (SDW). 4,5 The complexity in describing correlated particles led to the idea of simulating these systems with artificial architectures mimicking many-body models, such as the Hubbard Hamiltonian. Control over this Hamiltonian parameters requires a level of experimental precision achieved only recently with cold atoms in optical lattices. This experimental platform presents low tunneling probability for atoms in neighboring optical lattice sites, which sets a low energy scale for quantum effects. Optical lattice-based experiments obtained many-body manifestations as spin and charge correlations for 1D 5 and 2D 6,7 lattices. In the 1D case, related to the present study, the required temperature is in the nanokelvin range. Other promising proposals as the quantum simulation using a semiconductor quantum dot array were presented recently, 8 and will become a competitive architecture once long quantum dots arrays become feasible.
In the last few years, the expertise in positioning dopants nanostructures in Si has evolved. [9][10][11][12][13] The precision necessary for quantum applications with regard to impurity placement has pushed the development of these techniques to the point that atomic scale certainty is now a reality. [14][15][16] Precise atomic chains of these impurities are fabricated and, as demonstrated theoretically here, they may constitute convenient simulators for the extended Hubbard model. Multi-valley effects are ubiquitous in Si and valley interference impacts the tunnel coupling. We investigate-within a realistic model-how electronic properties of donor nanowires in Si can be controlled by design to emulate Hubbard systems, even allowing for some disorder effects.

RESULTS
We focus on the Hilbert subspace of neutral (half-filling) chains with N S = 8 sites, periodic boundary conditions and zero total spin projection S z tot ¼ 0 (sketched in Fig. 1). We arbitrarily set the quantization axis to z, without any regard to its significance with relation to the crystallographic directions. In the absence of external magnetic fields, this choice is arbitrary and the solutions to this Hamiltonian is invariant under a rotation of spin states. Donor positions are assigned at evenly spaced (R 0 ) substitutional atomic sites in Si along a [110] crystalline direction. The manybody state is described within the configuration interaction (CI) framework [17][18][19][20] and diagonalized exactly. Since Si is a material with very low spin-orbit coupling and no piezoelectric phonons, it is reasonable to assume that spin relaxation times are much longer than all other time scales involved in the experiment, so that thermalization does not remove the system from the S z = 0 subspace. Given a set of local operators {A i } acting on site i, a pair correlation function may be defined with i = 1 taken as the reference site where the average is taken over the thermally excited equilibrium occupations, N ' is the total number of states, w n is the Boltzmann weight of a level n at a given temperature T and Á Á Á iis the expectation value of the operator for Φ n j i, the n th eigenstate of H.
We define the dimensionless correlation function A j ðAÞ ¼ F 1;j ðAÞ=F 1;1 ðAÞ, so that for any T the self-correlation is A 1 ¼ 1, while A j ¼ 0 when the values of A at sites 1 and j are completely uncorrelated. It is now straightforward to define charge-charge C j ¼ A j ðρÞ Â Ã and spin-spin S j ¼ A j S z ð Þ Â Ã electronic correlations with the total spin at site j component along the quantization axis defined as S z j ¼ 1 2 c þ j;" c j;" À c þ j;# c j;# . Mermin-Wagner theorem 21 states that one-dimensional chains at finite temperatures can not sustain a long range order that breaks continuous symmetries. Still, the range of the pair correlations is a valuable figure of merit for the appearance of collective behavior in low-dimensional systems. We initially discuss charge correlations as a function of temperature T and interdonor spacing R 0 .
Results for the extreme cases T = 0 K and T = 300 K, are shown in the Supplemental Material. At absolute zero, C j is restricted to holes in the nearest neighbors of the reference site, namely j = 2 and j = 8, with strong localization due to the Mott mechanism. Any other pair {1, j} remains essentially uncorrelated. Results for 300 K are presented merely as an illustration of high temperatures limit, when even the nearest neighbors' correlations are lost. Neither the proposed device nor the model developed here are suitable for this temperature range. At room temperature, even the nearest neighbors correlations are lost. Figure 2a shows results at two experimentally achievable low temperatures. For T = 0.1 K the 0 K results are essentially reproduced while at T = 4 K and R 0~3 nm there is an alternation in C j among successive j's, i.e., a temperature-activated delocalization. This indicates that, for this particular R 0 , states with metallic character within k B T~0.3 meV of the ground state dominate the Boltzmann average. Spin correlations propagate further into the chain. The antiferromagnetic correlations among nearest neighbors hint at the establishment of a SDW phase, as expected for this range of parameters in the U vs V phase diagram. For T = 0 K, the antiferromagnetic-like behavior is observed along all j (see Supplemental Material), still with stronger correlations with neighboring sites (j = 2, 8). Results at room temperature, as for the charge, show no indication of correlation between sites. The T = 0.1 and 4 K results show a strong sensitivity of spin correlations with R 0 -some specific distances sustain the ground state antiferromagnetic tendency while others show very weak correlation signatures. This is a consequence of the oscillatory behavior of the hopping t with R 0 as can be observed in Fig. 3a. If the hopping is small (large), correlations will be weaker (stronger). Thus, for T = 0.1 K and R 0 = 6.53 nm the chain is fully correlated, while for R 0 = 6.91 nm the correlations vanish, reappearing for R 0 = 7.30 nm at this temperature.
Unavoidable positional disorder impacts all Hamiltonian parameters. We estimate this effect in the electronic properties through a simple model for disorder where P donors can occupy any position inside a disk with radius δ = 0.4 nm around a target substitutional site (see Fig. 4a). With this uncertainty radius, each donor can occupy 5 positions. This level of uncertainty is realistically achieved for STM placement of donors. 10,13 These nanowires are quasi-1D chains where electrons can follow only one path. 22 To investigate effects of disorder and temperature we compare results for perfect and disordered chains for two nominal distances between target sites (R 0 ). Both interdonor distances are chosen to be significantly larger than the range of disorder R 0 ) δ ¼ 0:4 nm. We choose the two distances sustaining (R 0 = 5.37 nm) and loosing (R 0 = 7.68 nm) the AF correlations at low T 4 0 (See Fig. 2b). Results for C j shown in Fig. 4b show that in all cases effects of the 2D disorder are mild, and do not affect significatively C j general trends-even for temperatures up to 4 K. Magnetic correlations, on the other hand, are more clearly affected by disorder, as seen in Fig. 4c. On average, disorder leads to less than 20% reduction for non neglectable correlations. The dispersion among spin correlations for individual realizations, indicated by the error bars in Fig. 4, means that disordered samples may present sizeable correlations, eventually stronger than the ordered ones. In real experiments, small chains are susceptible to this dispersion and may result in enhanced magnetic behavior.

DISCUSSIONS AND CONCLUSIONS
Perhaps the most important challenge for the implementation of the simulator described here is the measurement of the correlation function. While a direct measurement of charge and spin is possible 23,24 -these are the basis of the Kane model of quantum computation-it might be easier to extract these correlations from charge transport measurements. 25 The natural electronic correlations that appear in these chains may constitute an important resource for the study of many-body physics. It displays peculiar properties, constituting a unique example of a strongly interacting system with disordered tunnel coupling due to valley interference. This kind of random phase of the tunnel coupling element is the main ingredient in models displaying critical unitary statistics. 26,27 Moreover, the ongoing development of nanofabrication capabilities suggest that ondemand models may be analogically simulated. For instance, the intricate phase diagram of the Fermi-Hubbard problem may be unveiled by spin-resolved or density-resolved microscopy measurements. 5 Such application is under intensive investigation within cold atoms in optical lattices, and the present technology may complement these efforts. While not as easily tunable, the mass fabrication of circuits of donors adopting the know-how from available semiconductor technology would allow to chart the behavior of electrons over a wide range of attributes. The resilience of correlations in Si:P chains under relatively high temperatures suggests an attractive avenue for future experimental investigation.
We have shown here that, up to currently accessible values of position disorder and temperature, dopant arrays in silicon preserve quantum correlations among atoms in a diluted chain.
Our key point is that this system constitutes a robust implementation of the Fermi-Hubbard model in a semiconductor system with on-demand Hamiltonian parameters.

METHODS
These atomistic wires may extend throughout several nanometers, and a full description of the Si atoms would not be feasible. Instead, we describe the wire electronic states as a linear combination of donor orbitals (LCDO). 22 Each basis orbital is an effective mass Kohn-Luttinger (KL) where R i is the position vector of the substitutional donor at site i. The index μ = 1…6 labels the degenerate minima of the Si conduction band at k μ along the six equivalent 100 h i directions, i.e. ±x, ±y, ±z, , and a Si is the conventional Si lattice parameter. 29 In this approach, the influence of the Si host is explicitly included in the orbitals, allowing the investigation of longer chains than the conventional fully atomistic approach. When directly compared to experiments, this multivalley central cell corrected KL wavefunction leads to the correct single impurity spectrum, 30 single impurity wavefunction 31 and two impurities spectra, both in the ionized states 32 and neutral excited states. 33 We use isotropic hydrogen-like envelopes, F μ (r) = F(r) = (πa *3 ) −1/2 exp(−r/ a * ). These isotropic envelopes include a species-dependent Bohr radius a * obtained by considering a screened potential, affected by the Si host and the donor singular potential. Isotropic envelopes do not incorporate the effective mass anisotropy around the conduction band minima in Si, which is not relevant here. Its validity, including for the current system, is discussed in ref. 30 . Screening is treated by considering a functional form for the Coulomb potential that interpolates the expected asymptotic behaviors V(r → 0) = ±e 2 /4πϵ 0 r and V(r → ∞) = ±e 2 /4πϵ Si r, where the + (−) signal stands for electron-electron (electron-proton) potential. This transition between bare and screened point charge potentials occurs at a phenomenologically determined screening length r * , such that the full potential reads where ϵ Si is the Si static dielectric constant, and ϵ 0 the vacuum susceptibility. In what follows, we consider this screening both on single and two particle Hamiltonian terms. In the electron-electron Coulomb potential, we take r * = 0.1 nm, a typical value for a Si environment. 30 We study a first nearest neighbors Hamiltonian written in the LCDO basis. 22 Defining the creation and annihilation operators c þ i;σ and c i,σ for an electron at the orbital centered in R i with a spin projection σ along a quantization axis, and the corresponding number and charge density operators n i;σ ¼ c þ i;σ c i;σ and ϱ i = n i,↑ + n i,↓ , the Hamiltonian reads H ¼ This is readily recognizable as the extended Hubbard Hamiltonian, 34,35 with parameters ε i (onsite), t ij (hopping), U i (Hubbard) and V ij (direct). Analytic expressions for theses parameters are given in Supplemental Material, calculated values for a set of interdonor distances are presented in Fig. 3a. These parameters, which are consistent with typical orders of magnitude obtained experimentally, 9,11,13 support the idea that chains of dopants in Si constitute a strongly correlated system. Figure 3b shows that the ratio between the onsite Coulomb repulsion and the tunnel coupling is U/t ≈ 1 to 100, which ranges from the metallic regime, through the Mott insulator transition up to the strong localization driven by interactions. Still, the tunnel coupling t is strong enough that quantum fluctuations are dominant over thermal excitations even at dilution fridge temperatures T ≈ 100 mK. At this temperature, we have k B T/t ≈ 10 −4 to 10 −2 . Even at liquid He temperatures, this ratio is lower than 0.1 for all ranges of interdonor separation suggested here. In comparison, state-of-art cooling techniques applied to cold atoms still are not able to achieve ratios lower than k B T/t ≈ 0.2.
There is strictly no long range order in a one dimensional chain. Still, a rich variety of low-temperature electronic ordering tendencies appears at the range of parameters discussed here. Regardless of the strongly nonmonotonic behavior with R 0 , as a general trend small distances favor a CDW phase, while increasing R 0 we pass through a BOW phase and a SDW phase is favored at larger distances (see Supplemental Material). We investigate signatures of these many-body effects from charge and spin correlations.

Data availability
Supplementary information is available at npj quantum information website.

AUTHOR CONTRIBUTIONS
Numerical and analytical calculations were performed by A. Dusko., A. Delgado supervised the many-body methodology. A. Dusko, A. Saraiva and B. Koiller discussed the work together and wrote the manuscript. (c) Disordered Chain (δ=0.4nm) Perfect Chain (δ=0nm) Fig. 4 a Model of 2D disorder: the j th donor occupies the target (blue) or additional Si structure sites (green) within a distance δ from it. The figure corresponds to δ = 0.4 nm. b Disorder effects on charge C j À Á correlations at the indicated temperature, δ, R 0 and site. Square data points and error bars give averages and standard error calculated over an ensemble of 1000 8-site disordered chain realizations. c Same as b for S j correlations. In all cases we compare disordered (δ = 0.4 nm) and ordered (δ = 0 nm-circles) cases. Note that S j is more susceptible to disorder than C j Si:P chains as Fermi-Hubbard simulators A Dusko et al.