Superexchange coupling of donor qubits in silicon

Full configuration interaction (FCI) is a method to obtain the exact numerical solution of a many-body Schrödinger equation, limited by the number and quality of single electron basis states. The single electron basis states used here are calculated using a 10-band $sp^{3}d^{5}s^{*}$ atomistic tight-binding (TB) method in NEMO3D [41, 42]. This approach uses a localized atomic orbital-based method with nearest-neighbour interactions. The TB parameters are optimized to reproduce the bulk silicon band structure [43]. The phosphorus donors are represented using a Coulomb potential with a central cell correction at the donor site which can successfully determine the experimentally measured energy spectra of donors in silicon [44, 45].

The schematic representation of the simulations is illustrated in Figure 1(a). There are four phosphorus donors placed in a chain, each with an electron. The end electron spins, $Q_{1}$ and $Q_{2}$ are the qubits separated by R. The middle spins, $M_{1}$ and $M_{2}$ work as mediators with corresponding donors being separated by $r_{M}$. The separation between $Q_{1}$ and $M_{1}$ is $r_{1}$ and the separation between $M_{2}$ and $Q_{2}$ is $r_{2}$. The simulation domains used in the calculations entail $\sim 1.14$ million atoms which account for $\sim 60$ nm of silicon in the direction of separation and $\sim 20$ nm in the other two directions. The single electron basis states, solved from a parallel Block Lanczos eigensolver, are used to calculate the anti-symmetric Slater determinants of the multi-electron problem. The lowest two single electron valley-orbital states of the system are shown in Figure 1(b) and (c). Here we see bonding and anti-bonding state formation in the middle two donors $M_{1}$ and $M_{2}$ of a four-donor chain. The rest of the single electron molecular basis states are shown in Figure 8 in the Supplementary Information where we see the next two valley-orbital states are mainly localized in the outer donor dots.

The four-electron wavefunction is a superposition of various symmetry-permitted configurations of the Slater determinants. All possible integrals between the Slater Determinants with pairwise electronic interaction operators are computed to capture Coulomb, exchange, and higher-order correlations. The four electron Hamiltonian constructed from the Slater determinants is solved using the block Krylov-Schur algorithm within the Trilinos framework [46]. The number of single electron basis states in FCI calculations is increased until the eigenvalues of the four-electron system converge within a chosen tolerance – see Supplementary Information. On average, 56 single-electron basis states are sufficient to reach convergence. The eigenvalue of the ground state is separated from the triply degenerate excited states by the indirect exchange coupling. We place the donor atoms in different separations and orientations in our simulations and calculate this exchange coupling.

To analyze and interpret our FCI results, we compare them with an effective Hamiltonian model [30]. We can represent the four-spin system with a spin Hamiltonian such as -

where $\sigma_{i}$ is a Pauli matrix corresponding to an electron spin located on $i^{th}$ donor in the chain. $j_{M}$ is the exchange coupling between the middle two spins, $M_{1}$ and $M_{2}$. $j_{1}$ is the exchange coupling between $Q_{1}$ and $M_{1}$ and $j_{2}$ is the exchange coupling between $Q_{2}$ and $M_{2}$ – see Figure 1(a) for the schematic of the system.

We only consider the subspace with spin-zero states ($S_{z}=0$) to construct the Hamiltonian as we are interested in the singlet-triplet oscillations in the qubits $Q_{1}$ and $Q_{2}$. When $j_{1},j_{2}\ll j_{M}$, we can isolate the low energy states of the four-electron spin Hamiltonian using a Schrieffer-Wolff transformation. In that case, the effective Hamiltonian in a Heisenberg exchange form, $H_{SW}$ [30] now becomes -

In the regime where $j_{1},j_{2}\ll j_{M}$, the lowest energy eigenstates are characterized by the singlet state formed within the two middle dots and singlet or triplet states formed within the outer dots [30]. The energy difference, $\Delta E$ between the lowest long-distance singlet and triplet states, i.e.

is what we call superexchange ($\Delta E$) in this paper – see Figure 3(a) for the energy level diagram of the system. For $j_{1},j_{2}\ll j_{M}$ the middle-dot singlet manifold is well separated from all the middle-dot triplet states desirable for coherent coupling of the outer spins only, without any interference from the middle-dot spins.

We first analyze the equidistant case, where the separation between each neighbouring pair of donors is equal, $r_{1}=r_{2}=r_{M}=R/3$ – see Figure 1. Here we have approximately equal values of exchange coupling between each pair of donors, i.e. $j_{1}\approx j_{2}\approx j_{M}$. In Figure 2 (with red dots), we show the values of the indirect exchange coupling for these equispaced donors (4P chain) with constant but increasing separation between the dots along the $[100]$ (left) and $[110]$ directions (right) in the silicon crystal, calculated using FCI. For comparison, we also plot (with blue dots) the direct exchange values for two donors, i.e. 1P-1P system, separated by R, as previously calculated in [35]. Looking at the fitted dashed lines between the 4P chain and 1P-1P cases, we can see that the presence of the mediator donors dramatically increases the exchange coupling between the two outer spins. The indirect exchange coupling for a donor chain of $R$ of about 25 nm reaches $\sim 10^{2}$ GHz ($10^{-1}$ meV), while for direct exchange for the same separation without the presence of mediators would fall below $10^{-4}$ GHz ($10^{-7}$ meV) (extrapolated from the dataset in [35]).

We also see that the direct exchange decays much faster than the indirect exchange by comparing the slope of these two plots, true for both [100] and [110] crystal directions. For the [100] case in Figure 2(a), the slope of the fitted dashed line for the direct exchange as a function of qubit separation is $\sim-0.38$/nm whereas the slope for the indirect exchange is $\sim-0.12$/nm. Similarly in the [110] case in Figure 2(b), the slope of the fitted line for the direct exchange is $\sim-0.36$/nm whereas the slope for the indirect exchange is $\sim-0.11$/nm. In general, the indirect exchange decays almost three times less with donor separation than the direct exchange. This dependency with distance lets the qubits be far separated in the device and allows less correlated noise between them while maintaining strong spin coupling. Since the nearest neighbour exchange coupling has an exponential dependence with increasing separation of the donors, one might expect superexchange to change as an exponential function of $r_{1}+r_{2}-r_{M}$ – see Equation 2. On the contrary, the nearest neighbour exchange coupling changes as an exponential function of the qubit separation $R=r_{1}+r_{2}+r_{M}$ which results in a faster decay of direct exchange than the indirect one.

We observe no oscillations in the exchange coupling when the donors are separated in the [100] direction whilst we see oscillations in the exchange energy for separation in the [110] orientation, a signature of inter-valley quantum interference arising from the band structure of silicon [33, 35, 47].

To obtain superexchange and coherent control between the electron spin qubits at the end of the chains $Q_{1}$ and $Q_{2}$, it is essential that the middle spins $M_{1}$ and $M_{2}$ form a singlet-like state. Otherwise, there will be additional oscillations originating from the middle spins impeding coherent manipulation of the qubits. Thus we are looking for an operational regime where the ground state and the first excited state are well separated from the higher energy states in which the indirect coupling can be termed superexchange. In Figure 3(a) we plot the energies of all the eigenstates of the effective Hamiltonian $H_{eff}$ as a function of $j_{1,2}/j_{M}$ varying from 0 to 1. When the outer spins are weakly coupled with the middle spins ($j_{1,2}/j_{M}\sim 0$), we see two clear branches of energy states separated by $d$ but the separation between the branches decreases as $j_{1,2}/j_{M}\sim 1$.

In Figure 3(b), we plot the energy difference between the two lowest states, as a function of $j_{1,2}/j_{M}$ as calculated from $H_{eff}$ and $H_{SW}$. These two lowest states are the long-distance singlet and triplet states, namely $S^{l}$ and $T^{l}_{0}$. Here we see that, in the regime where $j_{1,2}/j_{M}\sim 0$, the solutions from the effective spin and the Schrieffer-Wolff Hamiltonian are the same since the assumption of the transformation is valid here. As $j_{1,2}/j_{M}$ increases, the two solutions start to diverge suggesting the Schrieffer-Wolff transformation no longer holds. The lowest two energy states (spanned by $\ket{\uparrow S\downarrow}$ and $\ket{\downarrow S\uparrow}$) are no longer well separated from the excited states and there are now triplet-like admixtures in the singlet-like state of the middle two spins.

In Figure 3(c), we can see the contributions in the ground state $S^{l}$ of the inner-dot singlet states (i.e. $\ket{\uparrow S\downarrow}$ and $\ket{\downarrow S\uparrow}$) and inner-dot triplet states (all the remaining basis states). When $j_{1,2}/j_{M}=1$ (i.e. the equidistant case where all donors are equally separated), the inner-dot singlet and triplet contributions in $S_{l}$ are both approximately 50% and we no longer have a well-defined two-level system of $S^{l}$ and $T^{l}_{0}$, but now have contributions from the four higher states shown in 3(a). As $j_{1,2}/j_{M}$ decreases, the contribution of the inner-dot singlet reaches more than $90\%$ where $j_{1,2}<\sim 0.3j_{M}$ (shaded region). This regime can be considered as a threshold for coherent operation of the qubits. However, we note that the threshold we mention here ($j_{1,2}/j_{M}\sim 0.3$) is not the upper limit of the coherent regime. The ratio should ideally be very close to $0$.

To explore the coherent-control regime of the 4-donor chain using our FCI calculations, we switch from the equidistant case discussed in Figure 2 to a chain with varied $r_{1,2}/r_{M}$ ratios, see Figure 4. Here, we present values of the indirect exchange coupling for chains oriented along the [100] and [110] direction in 4 (a), (b) respectively, as a function of middle donor separation while keeping the outer donor separation $R$ constant. The lower middle donor separation corresponds to a higher exchange coupling $j_{M}$ and lower $j_{1}$ and $j_{2}$ couplings, moving the system closer to the regime well described by the Schrieffer-Wolff approximation. In general, we see from Figure 4 that the exchange coupling exponentially decreases as we decrease the middle donor separation. The grey dashed line represents the separation where $r_{1}=r_{2}=r_{M}$. The exchange coupling shows two different dependencies on either side of this point. The first is where the distance between the middle donors is smaller than the values of $r_{1}$ or $r_{2}$ and the second is where it is large. For donors separated along the [110] direction, we see clear oscillations in $\Delta E$ as a function of $r_{M}$.

As we have discussed in the previous section the equidistant donor chain is not suitable for coherent manipulation of distant spins due to the admixture of states. The right side of the grey dashed line where $r_{1}=r_{2}\leq r_{M}$ is therefore not valid for operating qubits. To specify the limits of these different regimes, we use the values of direct exchange from the 1P-1P results in Ref. [35] to find $r_{1,2}$ and $r_{M}$ and the corresponding threshold. Decreasing the separation of the middle dots gives rise to a value of the superexchange that decreases exponentially while increasing the inner singlet contributions to the ground state – see Table 1 in the supplementary information. So there is a trade-off between the strength of coupling and operation fidelity. However, we see even small changes in the middle donor separation, of the order of $\sim$ 2 nm, can increase the inner singlet contributions to the ground state from $\sim 49\%$ to $\sim 98\%$ while keeping the value of superexchange still high enough ($\sim 100$ MHz or $\sim 1\times 10^{-3}$ meV) for use in realistic devices.

To address the tunability of superexchange, we have applied an electric field along the direction of the donor chain in our simulations where all donors are separated by $9.7758$ nm ($r_{1}$ = $r_{2}$ = $r_{M}$ = 9.7758 nm). For a realistic applied electric field of $2$ MV/m, we observe an exchange coupling of 18.09 GHz ($74.83\,\mu$eV) compared to 11.01 GHz ($45.556\,\mu$eV) under no electric field. The electric field detunes the qubits and increases the virtual tunnelling between them which results in slightly increased (less than a factor of 2) superexchange. However, when the donors are all charge neutral with one electron on each phosphorus atom, it is difficult to reach the (2,1,1,0) charge regime from a (1,1,1,1) regime without applying a very large bias. A similar challenge has been observed in [35] where applying an electric field of 2 MV/m the exchange coupling was modulated by 5 times. In the case of superexchange, the sensitivity of the coupling with an applied electric field is even less, hence, the distanced qubits are less susceptible to charge noise caused by electric field fluctuations compared to the nearest neighbour ones. To explore different levels of control over superexchange, we apply a J-gate bias on a system with a
phosphorus-doped silicon top gate based on the original Kane architecture [14]. We place gates 30 nm above the donor plane to control the potential barrier between donor atoms and the nearest neighbour exchange and observe their effect on superexchange. Similar three-dimensional tuning of the potential barrier between in-plane phosphorus-doped gates has been recently experimentally demonstrated [48] in precision engineered STM tunnel junctions where the gates were degenerately phosphorus-doped silicon layers. The schematic of the system is illustrated in Figure 5(a). The device consists of five quasi-metallic gates (source S, drain D, and three top gates - $G_{1}$, $G_{2}$ and $G_{M}$). $G_{1}$ and $G_{2}$ are used to tune $j_{1}$ and $j_{2}$ and $G_{M}$ is used to tune $j_{M}$. The total electrostatic potential profile of the device in the x-z plane is shown in Figure 5(b) where a cut is taken at the donor location. Here we see that, in the donor plane, the top gates create a parabolic potential profile. This total electrostatic potential is added to the tight-binding Hamiltonian while solving the energy levels of the donor quantum dots. In Figure 5(c-h), 1D cuts are taken at the donor locations shown for different $V_{M}$ and $V_{1}$(=$V_{2}$), i.e. for different voltages applied at gates $G_{M}$ and $G_{1}(G_{2})$, respectively. Here S and D are grounded, i.e. $V_{S}=V_{D}=0$. We note that just the presence of the electrostatic leads around the donors also modifies the total electrostatic potential of the device, even without an applied voltage. This is due to the band structure mismatch between the leads and the surrounding material, and the band-bending caused by the leads [49]. We calculate the total electrostatic potential profile of the device using a multi-scale modelling technique that combines the atomistic calculation of the band structure of the gates with a non-linear Poisson solution of the entire device [50].

Current STM lithography techniques allow us to place donor atoms in silicon with the precision of a single lattice constant [51, 52, 53, 54, 55]. As a consequence, we can investigate how small shifts in the exact P atom location within the four donor chain can affect the value of superexchange available. For that purpose, we analyze a chain oriented along the [100] direction with constant middle dot separation $r_{M}=5.431$ nm but with variable positions of the outer donors, described by $r_{1}$ and $r_{2}$ where $r_{1},r_{2}\geq$ 7 nm. We have deliberately chosen a smaller value of $r_{M}$ so that $r_{1},r_{2}>r_{M}$. We do not consider nearest neighbour separations less than 7 nm (qubit separation $R<$ 20 nm) since we are interested in the long-distance qubit coupling regime and want to exceed the distance achievable by direct exchange.

The results are shown in Figure 6. Here we can see the superexchange $\Delta E$ is maximum (11.53 GHz) for the smallest values of $r_{1}$ and ($r_{1}$ = $r_{2}$ = $7.6$ nm). This result arises both from minimizing the total chain separation, $R$ as well as maximizing the $r_{M}/r_{1,2}$ ratio – as observed in Figure 4. The value of $\Delta E$ decreases exponentially if either of the outer donors is moved outwards, i.e. if either $r_{1}$ or $r_{2}$ increases. However, at the same time, we can see $\Delta E$ does not change if both of the outer donors are shifted simultaneously in the same way – see the dashed lines where $r_{1}+r_{2}$ is constant. This result is consistent with the Schrieffer-Wolff approximation for superexchange in Eq. 2, where the dominating contribution to $\Delta E$ comes from $j_{1}j_{2}/2j_{M}$. While $j_{1}$ and $j_{2}$ change approximately exponentially with distance, their product will not change when $r_{1}+r_{2}$ is kept constant.

Thus we can conclude that any small shifts in the location of the donors within the [100]-oriented chain can result in large changes in total $\Delta E$. However, interestingly, a simultaneous shift of both outer (or both middle) donors in the same direction would not have any impact on the superexchange (indicated by the black dashed lines). The situation is completely different for the [110] crystalline direction due to the oscillations of two-donor exchange with distance which lifts the smooth $j_{1}j_{2}/2j_{M}$ dependency for varied $r_{1},r_{2}$ when $r_{1}+r_{2}$ is constant. Here $\Delta E$ will change in an oscillatory manner if one or both middle donors are shifted.

Finally, we comment on the presence of nuclear spins in the donor system and how they can also impact superexchange and coherent electron spin manipulation. The nuclear spins of phosphorus donor atoms couple with their electron spins through the hyperfine interaction. From the perspective of the electron spin, this can be treated as a small additional magnetic field dependent on the nuclear spin polarization. In the case of an electron singlet-triplet spin qubit, localized within the middle double donor dot, this hyperfine interaction can create a magnetic field gradient that mixes singlet and triplet states [17, 56]. For a 1P-1P system, the gradient is $\sim$ few Milli Tesla which gives rise to a Zeeman energy difference, $\Delta E_{z}$ $\sim$ $5\times 10^{-4}$ meV between the donors when the two nuclear spins are oriented antiparallel ($\Uparrow\Downarrow$) and $\Delta E_{z}$ $\sim 0$ when the two nuclear spins are aligned in the same direction ($\Uparrow\Uparrow$).

To create well-defined eigenstates as well as singlet and triplet states, the exchange coupling $\Delta E$ needs to dominate over $\Delta E_{z}$. We can consider the 4-electron system discussed in this paper essentially as two pairs of singlet-triplet qubits – one defined between the two inner dots and one between the two outer dots. To avoid mixing of the singlet and triplet states it is desirable for both $j_{M}$ and $\Delta E$ to be significantly greater than $\Delta E_{z}$. From Figure 2 this is satisfied when the middle donor separation $r_{M}$ is smaller than 10-15 nm, which includes most of the results presented in this work. The value of superexchange $\Delta E$ ultimately depends both on $r_{M}$ and $R$. From Figure 4 we can see it is possible to design 4P configurations which belong to the regime where $j_{1},j_{2}/j_{M}<0.3$ and at the same time satisfy $\Delta E>\Delta E_{z}$. These considerations are however no longer relevant if we deterministically initialize the nuclear spins to all-parallel before any operation using a local Nuclear Magnetic Resonance (NMR) [10]. In this case, there is no magnetic field gradient and thus no limitations on the minimum values of $j_{M}$ and $\Delta E$.