Adequacy of Si:P chains as Fermi–Hubbard simulators

Abstract

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 R₀, 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 R₀. 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.

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⁵ 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,⁸ 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_{\rm {tot}}^z = 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₀) substitutional atomic sites in Si along a [110] crystalline direction. The many-body 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.

Fig. 1

Representation of a 1D chain of N _S dopants in Si aligned along the [110] direction and with interdonor separation R₀. We adopt periodic boundary conditions, connecting site j = N _S with site 1

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

$${\cal F}_{1,j}(A) = \mathop {\sum}\limits_{n = 1}^{N_\ell } w_n\left( {\left\langle {A_1A_j} \right\rangle _n - \left\langle {A_1} \right\rangle _n\left\langle {A_j} \right\rangle _n} \right)$$

(1)

where the average is taken over the thermally excited equilibrium occupations, $N_\ell$ is the total number of states, w _n is the Boltzmann weight of a level n at a given temperature T and $\left\langle \cdots \right\rangle _n = \left\langle {{\mathrm{\Phi }}_n\left| \cdots \right|{\mathrm{\Phi }}_n} \right\rangle$ is the expectation value of the operator for $\left| {{\mathrm{\Phi }}_n} \right\rangle$, the n^th eigenstate of H.

We define the dimensionless correlation function ${\cal A}_j(A) = {\cal F}_{1,j}(A){\mathrm{/}}{\cal F}_{1,1}(A)$, so that for any T the self-correlation is ${\cal A}_1 = 1$, while ${\cal A}_j = 0$ when the values of A at sites 1 and j are completely uncorrelated. It is now straightforward to define charge-charge $\left[ {{\cal C}_j = {\cal A}_j(\rho )} \right]$ and spin-spin $\left[ {{\cal S}_j = {\cal A}_j\left( {S^z} \right)} \right]$ electronic correlations with the total spin at site j component along the quantization axis defined as $S_j^z = \frac{1}{2}\left( {c_{j, \uparrow }^ + c_{j, \uparrow } - c_{j, \downarrow }^ + c_{j, \downarrow }} \right)$.

Mermin–Wagner theorem²¹ 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₀.

Results for the extreme cases T = 0 K and T = 300 K, are shown in the Supplemental Material. At absolute zero, ${\cal 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₀ ~ 3 nm there is an alternation in ${\cal C}_j$ among successive j’s, i.e., a temperature-activated delocalization. This indicates that, for this particular R₀, states with metallic character within k_BT ~ 0.3 meV of the ground state dominate the Boltzmann average.

Fig. 2

Electron-electron correlations with respect to site j = 1 (therefore ${\cal C}_1 = {\cal S}_1 = 1$, see text). a Calculated charge-charge correlations for the indicated temperatures. Each box shows ${\cal C}_j$ (see color code) as a function of interdonor distance (R₀) and position along the chain (d). Note that for T = 4 K and R₀ ~ 3 nm the charge-charge correlations are temperature activated. b Same as a for calculated spin–spin correlations. An anti-ferromagnetic (AF) behavior up to fourth neighbors is clearly observed—the fifth, sixth and seventh neighbors pairs are equivalent to third, second and first pairs, respectively, due to periodic boundary condition. Note that the magnetic behavior is not monotonic with R₀

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₀—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₀ 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₀ = 6.53 nm the chain is fully correlated, while for R₀ = 6.91 nm the correlations vanish, reappearing for R₀ = 7.30 nm at this temperature.

Fig. 3

a Absolute values of the calculated Hamiltonian parameters Hubbard (U), direct (V), multi-valley hopping (t) and onsite energy (ε). Spheres mark the values of R₀ allowed for donor pairs in Si along [110]. b Relative values of the Hubbard energy and k_BT at temperature T = 100 mK with respect to the hopping absolute value. The sharp oscillations are due to |t| alone, as the others are constants. The temperature of 100 mK is attainable for experiments performed under dilution refrigerators. Small $k_BT{\mathrm{/}}\left| t \right|$ and large $U{\mathrm{/}}\left| t \right|$ favor the experimental implementation of the proposed simulator (see text)

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

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 $\left( {{\cal C}_j} \right)$ correlations at the indicated temperature, δ, R₀ 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 ${\cal S}_j$ correlations. In all cases we compare disordered (δ = 0.4 nm) and ordered (δ = 0 nm—circles) cases. Note that ${\cal S}_j$ is more susceptible to disorder than ${\cal C}_j$

These nanowires are quasi-1D chains where electrons can follow only one path.²² To investigate effects of disorder and temperature we compare results for perfect and disordered chains for two nominal distances between target sites (R₀). Both interdonor distances are chosen to be significantly larger than the range of disorder $R_0 \gg \delta = 0.4$ nm. We choose the two distances sustaining (R₀ = 5.37 nm) and loosing (R₀ = 7.68 nm) the AF correlations at low T > 0 (See Fig. 2b). Results for ${\cal C}_j$ shown in Fig. 4b show that in all cases effects of the 2D disorder are mild, and do not affect significatively ${\cal 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.²⁵

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 on-demand 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.⁵ 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).²² Each basis orbital is an effective mass Kohn–Luttinger (KL) variational wavefunction for the ground state (A₁ symmetry),²⁸

$${\mathrm{\Psi }}_{{\mathrm{KL}}}({\bf{r}}) = \frac{1}{{\sqrt 6 }}\mathop {\sum}\limits_{\mu = 1}^6 F_\mu \left( {{\bf{r}} - {\bf{R}}_{\bf{i}}} \right)\phi _\mu \left( {{\bf{r}} - {\bf{R}}_{\bf{i}}} \right),$$

(2)

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 $\left\langle {100} \right\rangle$ directions, i.e. ±x, ±y, ±z, $\left| {{\bf{k}}_\mu } \right| = k_0 = 0.85\left( {\frac{{2\pi }}{{a_{Si}}}} \right)$, and a _Si is the conventional Si lattice parameter.²⁹ 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,³⁰ single impurity wavefunction³¹ and two impurities spectra, both in the ionized states³² and neutral excited states.³³

We use isotropic hydrogen-like envelopes, F _μ (r) = F(r) = (πa^*3)^−1/2exp(−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. ³⁰. Screening is treated by considering a functional form for the Coulomb potential that interpolates the expected asymptotic behaviors V(r → 0) = ±e²/4π$\epsilon _{\mathrm{0}}$r and V(r → ∞) = ±e²/4π$\epsilon _{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

$$V(r) = \pm \frac{{e^2}}{{4\pi r}}\left[ {\frac{1}{{\epsilon _{Si}}} + \left( {\frac{1}{{\epsilon _0}} - \frac{1}{{\epsilon _{Si}}}} \right)e^{ - r/r^ \ast }} \right],$$

(3)

where $\epsilon _{Si}$ is the Si static dielectric constant, and $\epsilon _{\mathrm{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.³⁰

We study a first nearest neighbors Hamiltonian written in the LCDO basis.²² Defining the creation and annihilation operators $c_{i,\sigma }^ +$ 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,\sigma } = c_{i,\sigma }^ + c_{i,\sigma }$ and ϱ _i  = n_i,↑ + n_i,↓, the Hamiltonian reads

$$H = \mathop {\sum}\limits_{i,\sigma } \varepsilon _in_{i,\sigma } + \mathop {\sum}\limits_{\langle i,j\rangle ,\sigma } t_{ij}c_{i,\sigma }^ + c_{j,\sigma } + \mathop {\sum}\limits_i U_in_{i, \uparrow }n_{i, \downarrow } + \mathop {\sum}\limits_{\langle i,j\rangle } V_{ij}\varrho _i\varrho _j.$$

(4)

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_BT/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_BT/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 non-monotonic behavior with R₀, as a general trend small distances favor a CDW phase, while increasing R₀ 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.