Rydberg Entangling Gates in Silicon

In this paper, we propose a new Rydberg entangling gate scheme which we demonstrate theoretically to have an order of magnitude improvement in fidelities and speed over existing protocols. We find that applying this gate to donors in silicon would help overcome the strenuous requirements on atomic precision donor placement and substantial gate tuning, which so far has hampered scaling. We calculate multivalley Rydberg interactions for several donor species using the Finite Element Method, and show that induced electric dipole and Van der Waals interactions, calculated here for the first time, are important even for low-lying excited states. We show that Rydberg gate operation is possible within the lifetime of donor excited states with 99.9% fidelity for the creation of a Bell state in the presence of decoherence.

Donor electron spin qubits in silicon are competitive in most of Di Vincenzo's criteria for a quantum computer implementation: scalability (and CMOS compatibility), high-fidelity initialisation and readout as well as extremely long coherence times [1][2][3][4][5]. However, published experimental results of entangling gates in silicon report only up to 86% [6] or 90% [7] fidelity. Moreover, entangling exchange interactions can change by orders of magnitude if the distance between orbital ground-state spins is varied by one lattice site [8]. Due to the impact on gate durations, this renders scaling challenging as even hydrogen lithography cannot to date perform donor placement accurate to within a single silicon lattice site [9,10].
In this paper we draw on the analogy between donors trapped in the silicon lattice and alkali Rydberg atoms trapped in light fields [12] to propose a new Rydberg entangling gate which is robust to variations in the interqubit distance, with theoretical fidelities up to 99.9% in silicon. Our gate relies on the enhancement of the interactions u between atoms in an orbital excited (Rydberg) state |r which leads to rapid gate operations on the order of hundreds of picoseconds. We show that the theoretical fidelities achievable by our new gate protocol are an order of magnitude higher than existing Rydberg atom schemes. This is because our gate allows for the excited state interaction u to be on the same order of magnitude as the Rabi frequency Ω coupling one of the qubit levels |1 spin-selectively (see Fig. 1) to |r .
Specifically, we show here that while in cold atoms excitation to high-lying excited states is needed to induce strong dipole interactions, for donors in silicon the interactions generated by the low-lying excited states are large enough due to the large dielectric screening and small effective mass in silicon. We calculate for the first time * e.crane@ucl.ac.uk . Global gates (blue) enable use of induced dipolar interactions. Readout devices can be placed on the layer below the one shown here. This scheme can be extended for a surface code quantum computer similarly to [11]. b) The ground state hyperfine levels are used as a qubit basis, with a higher lying orbital excited state as |r . c) The interactions between two atoms in the excited state barely oscillate, leading to a distance-robust two-qubit phase acquired after a Rabi cycle (residual oscillations are due to multivalley interference in the exchange interaction).
entangling gate operation in the presence of a finite excited state lifetime T 1 , assuming the decoherence time to be twice the lifetime, and show that our proposal outperforms existing Rydberg schemes. By using the experimentally determined T 1 and our numerically determined u, we find that high fidelities independent of placement errors can be reached in silicon. Lastly, we show that this gate is robust to variations in donor energy levels (inhomogeneous broadening) and that only tuning of a single independent parameter is needed. Rydberg entangling gates.-Whereas previous gate schemes [21][22][23] required /T 1 Ω < u, our proposal only requires /T 1 Ω, u, allowing for much smaller gate durations and hence less probability of decay from |r , leading to higher gate fidelities. This is because previous schemes were based on the blockade, in which the doubly excited state |rr is tuned out off resonance due to strong interactions induced by Ω u. Therefore, the gate duration τ ∼ 1/Ω is limited to τ 1/u. In our protocol, due to Ω being on the same order of magnitude as u, the gate duration can be τ ∼ 1/u.
The pulse sequence (see supplement for details [24]) consists of acting twice with a pulse with Rabi frequency |Ω| detuned from the |r transition on both donors simultaneously for a duration chosen such that the first pulse returns |01 and |10 to themselves. The phase ξ of Ω is then chosen such that |11 returns to itself after a second pulse of the same duration and Rabi frequency |Ω|. During the two pulses, |11 picks up a two-qubit phase due to the interactions in |rr . Lastly, the detuning ∆ is chosen such that the final phases of all two-qubit states implement a controlled phase gate. Global addressing enables parallel entangling gate operations as demonstrated in Ref. [23]. Conditions for a successful gate implementation are: spin-selective excitation to |r (either leaving the |0 − |r transition off resonance or forbidden by a selection rule) and u, Ω > /T 1 . In the following, we show that these two conditions can be fulfilled in both shallow and deep donors.
Spin-selective excitation to |r .-The singly ionised deep double donor selenium ( 77 Se + ) has a large hyperfine interaction in the ground state, giving rise to singlet S0 and triplet T0 states as the qubit basis with T 1 beyond 6 minutes and T 2 beyond 2 s [5]. For |r both the ± valley states of 1sT 2 Γ 7 [5] and 2p0 [14] can be used, as they are dipole-allowed spin-selective transitions with easily accessible excitation energies.
For the shallow donors phosphorus (P), arsenic (As), antimony and bismuth, the qubit can be encoded in the hyperfine S0 and T0 states with a T2 time exceeding 7s [3]. For |r we take the 2p0 state, which can be accessed from the D0X level [3,16] or using micromagnets to emulate spin-orbit coupling, as the larger extent of the 2p0 wavefunction would lead to a larger Zeeman splitting in |r than in |g [25]. Higher lying states would give rise to larger u, at the cost of a smaller Ω or the need for a two-photon excitation, as is the case for cold atoms.
For all proposals, we have checked that Ω > h/T 1 is feasible with laser intensities below 1 W/µm 2 , and we list their parameters in table I.
Having discussed that spin-selective excitation of orbital-excited states is indeed possible, we now calculate the interactions between donors in |r to find u.
Rydberg interaction u.-Entangling gate schemes so far have focussed on harnessing the highly oscillatory and exponentially decaying ground state ferromagnetic exchange (J) [6,26]. This interaction is negligible in high-lying excited states of Rydberg atoms, where u is dominated by the electric dipole interaction [21] where p i is the dipole operator of atom i, R is the interatomic distance, and n is the unit vector pointing along the inter-donor axis. In the absence of an electric field, neutral atoms and donors do not have a permanent dipole moment so V dd only contributes at second order in perturbation theory: Van der Waals (V VdW ) interactions falling off as ∼ 1/R 6 , see Supplemental Material [24] whose strength is determined by the polarisability which scales with the principal quantum number as n 11 . Contrarily, J is determined by the size of the wavefunctions which scale as n 2 and decay exponentially with interdonor distance. Therefore, dipolar interactions dominate for large principal quantum number and distance. In silicon, the six conduction band minima lead to donor electron multivalley wavefunctions. Here, we calculate multivalley results for J, electron-electron repulsion (W ) and for the first time also V VdW and the induced dipole interaction V dd in the presence of an electric field. We employ the FEM to solve the Schrödinger equation exactly within effective mass theory: where U cc is a central-cell potential for the S states which have a non-negligible probability of being at the core, U Ef is the electric field, i is the binding energy of the donor and m * is the effective mass of the electron in the silicon lattice, see Supplemental Material [24]. Then we couple the result into the six valleys. The FEM solves for any number of eigenstates, enabling a convergence check of perturbative calculations. For V VdW interactions, we typically include ∼ 40 states. Similarly, we obtain V dd from the dipole moment which we get from the slope of the multi-valley Stark shift. This perturbatively calculated Stark shift does not contain contributions from the continuum, which we quantify by comparing to the single valley energies directly generated by the FEM while varying U Ef (see Supplemental Material [24]). As the outermost donor electron is very weakly bound to the core, we set the maximally applicable electric field (which wouldn't ionise |r during its lifetime) by generalizing the method in ref. [27] to donors in silicon, see Supplemental Material [24], consistent with experimental results. Fig. 2 shows all interactions between two 2p0 state phosphorus donors and the relevant Rydberg interaction u which corresponds to the additional energy experienced by the atoms when they are doubly excited as opposed to singly excited: where the subscript |g denotes the ground state. J |gg , V VdW |gg , V VdW |rg , J |rg are negligible for the range of donor separations considered in this work. Contrary to the expectation from the ground orbital state and despite the small principal quantum number, both V VdW and V dd interactions are non-negligible. When two donors are aligned along the polarisation and electric field axis as is shown in figure Ia), V VdW and V dd dominate, with an oscillatory contribution from J. When aligned along an axis perpendicular to the polarisation, J and V VdW dominate, with no oscillations due to the 2p0 state not showing multivalley oscillations in the plane perpendicular to the polarization axis (see Fig. 2c)). In the Supplemental Material [24], we show that for the 1sT 2 Γ 7 state in Se+ the dipole interactions are less dominant and u is mainly given by J.
For a successful gate operation in a quantum computer, the interactions u need to be switchable witout affecting the surrounding qubits. Selecting the qubit pair to be excited to |r requires Stark shifting the donors into resonance with the laser using electric fields [6,28], analogously to proposals for Rydberg atoms in optical lattices [29]. Single qubit operations can be done in the same way, using global microwave illumination [28]. We performed an electrostatic simulation of the gate configuration sketched in Fig. 1a), showing in Fig. 2 that the electric gate configuration sufficient for shifting donors off-resonance would not create stray electric fields harmful to device operation. It would also be possible to apply an electric field roughly parallel to the inter-donor axis, In all donors and |r studied here, u is on the order of 1-10meV, fullfilling u > /T 1 . Given a large enough laser power, the interaction scale sets the ultimate limit to gate operation times, enabling entangling gate durations of a fraction of a nanosecond, nine to ten orders of magnitude faster than the qubit lifetime.
In order to show that high-fidelity entangling gates are possible within the short donor lifetime, we calculate the fidelity of Bell state creation for various donors in the presence of decoherence by simulating the whole pulse sequence within a Lindblad Master equation approach in the following.
Fidelity of Bell-state creation in presence of decoherence.-We model a two donor system with the Hamiltonian [22] and simulate the full pulse sequence for creating the Bellstate |Φ + = 1 √ 2 |00 + |11 within a Markovian Lindblad Master equation given by where ρ is the two-donor density matrix (with each donor restricted to the states |0 , |1 and |r ), L j are the jump operators describing decoherence and H eff = We include dephasing between |1 and |r with rate 1/(2T 1 ) and spontaneous emission from the Rydberg state, see Supplemental Material [24], with rate 1/T 1 . Decoherence processes between the qubit levels are negligible on the time scale of a single two-qubit gate.
We compare the Bell-state creation fidelity F = Φ + | ρ |Φ + between the previous schemes and our proposal, see Supplemental Material [24] and Fig. 3. According to the hierarchy of scales /T 1 Ω < u present in the previous schemes discussed above, there is a clear optimum for Ω [30]. In the gate proposed in this paper, there is also an optimal Ω which we determine numerically along with the other pulse parameters (see supplement [24]). We find that the only parameter which requires tuning is Ω, which depends linearly on u: For a given u/γ, the optimal Rabi frequency for our proposal is an order of magnitude larger than in the previous schemes, leading to much shorter gate durations and higher fidelities. We also find the fidelity to be robust to small variations in Ω and u. It also remains unchanged under realistic inhomogeneous broadening, which would lead to a variation of the excited state energies and therefore the detuning, see Supplemental Material [24]. We use these results to predict the Bell-state creation fidelities possible in silicon by entering the interactions u calculated in the previous section and the lifetimes summarized in Table I. In Fig. 3b) we show the gate fidelity as a function of the distance between donors. Importantly, fidelities of approximately 99.9% can be reached. Moreover, the fidelity is only weakly dependent on distance such that placement errors as large as several nanometers make little difference and it does not make a difference if both donors reside on different sublattices of the silicon lattice (see supplement [24]) To conclude, we have proposed a new Rydberg gate, potentially increasing fidelities in cold atom experiments by an order of magnitude. Furthermore, we have shown that they can be implemented in donors in silicon, reaching fidelities of 99.9% and having a high tolerance for donor placement error, opening the possibility to implement entangling gates in solid state quantum computers with ion implantation [33]. The insensitivity to variations in the interaction strength and the inhomogeneous broadening not only brings advantages for scalability, but also reduces fine-tuning in device operation and may therefore be advantageous for other solid-state platforms. We emphasize that the physical nature of the interactions is irrelevant, in particular both exchange and dipolar interactions may be used. Acceptors do not have multivalley oscillations [34] and would therefore be ideal candidates for this robust gate implementation as their dipolar interactions may be large enough to obtain high fidelities Predicted Bell-state creation fidelities for general Rydberg gates (a) and for Si:P or Si:As (b) a) Results from the Lindblad simulation of the pulse protocols from Jaksch et al. [22] and Levine et al. [23] along with the protocol proposed in this paper and the experimental data on alkali atoms in optical tweezers [23,31]. b) Fidelity map (plotted for 95% and above) as a function of distance between donors shows that Fidelities are insensitive to donor placement. Experimental lifetime T 1 = 235ps as shown in Tab.I. |r is the 2p0 state polarised along {0,0,1}. The same fidelity region appears, diametrically opposed with regards to donor 1. 8nm minimum separation for avoiding the molecular limit in the 2p0 state [32]. Fidelity robust to donors belonging to different silicon sublattices (see Supplemental Material [24]) or diffusing out of plane (Si:As diffusion <1nm [10] First, we present the Lindblad master equation approach to simulating the whole pulse sequence on two qubits in the presence of dephasing and spontaneous decay. Second, we present the original resonant gate and its optimisation [1]. Third, we present the recently proposed ultrafast off-resonant blockade gate [2] discussed in the main text. Finally, we present our gate which shows much higher fidelities than the two previous proposals.
To benchmark the gates, we calculate the fidelity of creating the Bell-state |Φ + = 1 √ 2 |00 +|11 by simulating the whole pulse sequence in the presence of decoherence. We model a two-donor system with Hamiltonian where Ω i , ∆, u are Rabi frequency of the laser acting on donor i, detuning between the laser frequency and the transition frequency between |1 and |r , and interaction strength between the two donors, respectively. We restrict the donor energy levels to the two qubit levels |0 , |1 and the Rydberg state |r . In order to take into account dephasing and spontaneous emission from the short-lived Rydberg state, we solve a Markovian Lindblad Master equation for the density matrix ρ of the form where H eff = H − i 2 j L † j L j and L j are the jump operators given by where γ se , γ de denote the rates of spontaneous emission of phonons and dephasing between the excited state |r and the qubit level |1 . We assumed spontaneous emission does not bring |r to |0 and we neglect decoherence processes between the qubit levels as these are negligible on the time scale of a single two-qubit gate. We solve the evolution for arbitrary decoherence rates and interactions and express all results in units of the decoherence rate γ se , assuming that γ de = 0.5γ se . In order to deduce the necessary interaction strength and Rabi frequency for a successful gate operation, the corresponding value for γ se given the donor/atom species has to be inserted as discussed in the main text.
We start in the initial state |Ψ 0 = (|1 + |0 )/ √ 2 ⊗ (|1 +|0 )/ √ 2, and entangle the two qubits by numerically solving the time evolution of the pulse sequence. Finally, we apply a single qubit phase gate to rotate to |Φ + [2]. We assume the latter to be perfect, i.e. we act with the appropriate unitary operator on the density matrix instead of simulating the pulse. We use the fidelity of creating |Φ + , defined as as a measure how much decoherence alters the success of the gate.

arXiv:2008.11736v2 [quant-ph] 24 Nov 2020
Resonant blockade gate [1] Pulse sequence |11 |10 |01 |00 1. π-pulse on atom 1 i|r1 i|r0 |01 |00 2. 2π-pulse on atom 2 i|r1 i|r0 i|0r |00 3. π-pulse on atom 1 -|11 -|10 -|01 |00 The resonant Rydberg blockade gate requires the qubits to be individually addressable and the laser to be on resonance with the orbital transition. The pulse sequence (c.f. Fig.1) applies a π pulse to the first atom, a 2π pulse to the second atom and a π pulse to the first atom again. In the initial state |00 the pulse sequence has no effect on the qubits because the |0 state is offresonant with the laser. Due to the Rydberg blockade, |11 acquires the same phase as |10 , as the second atom cannot be excited to |r , in total implementing the truth table of a controlled-Z gate.
In the presence of a non-zero decoherence rate γ, there is an optimum for the Rabi frequency Ω: If Ω ≈ γ, the fidelity is low because the pulses take too long, however if Ω ≈ u, the fidelity is also low as the blockade condition is not well fulfilled. We hence optimize the Rabi frequency Ω to maximize the fidelity F.
In Fig. 1 we show that this intuitive picture is correct, i.e. there is a clear optimum value for the Rabi frequency. Interestingly, the optimal Rabi frequencies are rather large showing that the gate operation time should be as small as possible to reduce loss from the excited state. This means that the hierarchy of scales is given by γ Ω < u. Importantly for the purposes of the donor implementation, we also find that for a given value of the Rabi frequency, the fidelity is only weakly dependent on the exact value of u. For large u, a plateau is reached beyond which the fidelity does not increase due to decoherence processes then being the limiting factor. At this point, the Rabi frequency should be increased to increase fidelities further.  [2]. a) Global pulse sequence acting with |Ω| on both qubits simultaneously: t1 pulse duration chosen such that the first pulse returns |11 to itself (whereas |01 and |10 are left in an arbitrary location on the Bloch sphere), ξ phase of the second pulse (evolve with Ω exp(iξ)) chosen such that |01 and |10 return to themselves after the pulse sequence, ∆ detuning of the |1 to |r transition chosen such that the phases highlighted in b) are equal. Finally, the third step is to apply single qubit phase gates to both qubits, with phase φ, which corrects for global single qubit phases built up in the dynamics. b) Mapping of the qubit states due to the Cz gate. c) The dynamics of the states |11 (with Rabi frequency √ 2Ω in the case of a perfect blockade) and |01 (with Rabi frequency Ω, leading to a different path over the Bloch sphere) in terms of two level systems. d) Optimal Rabi frequency in the off-resonant blockade gate as a function of interaction strength and decoherence time.
In the improved blockade gate as presented and implemented in Ref. [2], only two global pulses with fixed detuning ∆ = 0 are needed. The pulse sequence proceeds as follows: 1. Evolve for time τ with Ω 1 = Ω 2 = Ω.
2. Evolve for time τ with Ω 1 = Ω 2 = Ω exp(iξ), again followed by a single qubit phase gate on both qubits (c.f. Fig. 2a)). In the above sequence, the gate time τ is chosen such that the state |11 returns to itself after the first pulse, i.e. to the time corresponding to a blockaded 2π pulse, τ = 2π/ √ 2Ω 2 + ∆ 2 . Similarly, the phase ξ of the laser in the second pulse is chosen such that the state |01 returns to itself after both pulses. Finally, the detuning ∆ is chosen such that the phase acquired by both states differs by π (in the sense that this phase difference remains after application of the global single qubit phase gate, c.f. Fig. 2b).
In the strongly interacting regime, the parameters were analytically calculated to be given by [2] ∆/Ω = 0.377371 This leaves the Rabi frequency as a free parameter, which we optimize for a given interaction strength u and decoherence rates γ se and γ de = 0.5γ se . Similarly to the resonant gate, we find a clear optimal Rabi frequency for a given interaction strength, which we show in Fig. 2 and find slightly higher fidelities than for the resonant gate (see main text). A major limitation for both blockade gate proposals is the hierarchy of scales γ Ω < u, given by the requirement for the blockade condition to be fulfilled. This limits the maximum Rabi frequency, setting a lower bound for the gate duration and leaving more time for decoherence processes to kick in. Our calculation implies that the experiment in Ref. [2] operated very close to the optimal Rabi frequency, in part explaining their improvements over previous results.

Blockade-inspired off-resonant gate
We propose a "phase accumulation" variant of the above gate, similar to the "Gate A" discussed in Jaksch et al. [1]. To do so, we use the same pulse sequence as in the off-resonant blockade gate but reverse the reasoning -instead of tuning the gate time to the blockaded Rabi frequency, we use the unblockaded one, given by τ = 2π/ √ Ω 2 + ∆ 2 . Hence, the states |01 and |10 return to themselves after a single Rabi pulse. Accordingly, we tune the phase ξ such that the state |11 returns to itself after the second pulse, accumulating a two-qubit phase due to the non-vanishing probability to populate Pulse sequence of proposed gate a) Global pulse sequence acting with |Ω| on both qubits simultaneously: pulse duration τ chosen such that the first pulse returns |01 to itself (whereas |11 is left in an arbitrary location on the Bloch sphere), ξ is the phase of the second pulse chosen such that |11 returns to itself after the second pulse, ∆ detuning of the |1 to |r transition chosen such that the phases highlighted in b) are equal. c) Bloch-sphere depiction of the dynamics of the states |11 and |01 . The |11 state touches the Bloch sphere spanned by |11 and |rr at times (n/2)τ with n = 1, 2, 3, 4. This is more clearly shown in d), which depicts the state probabilities when starting in the initial state |ψ(0) = (|0 + |1 ) ⊗ (|0 + |1 )/2. e) Optimized pulse parameters in the blockade-inspired off-resonant gate. Note the very weak dependence on the interaction strength.
|rr . That this proposal is advantageous compared to the blockade gates may seem counterintuitive at first due to the high loss rate from |r . However, in this gate the condition Ω < u is relaxed and hence the gate duration can be much smaller compared to the blockade gates. We optimize the pulse sequence and find that the optimal Rabi frequency is proportional to the interaction strength. Furthermore, all other parameters to be almost independent of the interaction strength. We obtain these results numerically by fixing τ = 2π/ √ Ω 2 + ∆ 2 and optimizing Ω, ∆, ξ for a given u/γ, obtaining the optimal pulse parameters shown in Fig. 3. In the large interaction limit, we find Ω/u = 1.45747, Note that the optimal Rabi frequency is about an order of magnitude larger than in the blockade gate schemes (given a fixed value of the interaction strength), leading to an order of magnitude smaller gate durations and to the much larger fidelities shown in the main text.
Robustness to placement error We checked that a displacement within the 2x3 lattice site physical limit on the precision placement of hydrogen lithography of phosphorous in silicon [3] would still be within the width of the "fidelity resonance", as would be a small diffusion out of plane, see Fig. 4. This result implies that a sample made with hydrogen ligthography would not need tuning -assuming the interaction between donors as a function of distance is known, and for an ion implaneted sample, a simple procedure to finding good pulse parameters could be a rough scan of the sample with an scanning tunnelling microscope (STM) after fabrication to determine the locations of the donors. Note that while placement-insensitivity is stronger in the blockade gates, there is still an error introduced by Ω/u differing from the optimal value if u varies. Negligible effect of inhomogeneous broadening.-Infidelity as a function of interaction strength scaled by |r lifetime, for a fixed Ω and interaction strength u. Even a relatively large inhomogeneous broadening of 0.1Ω corresponding to roughly 0.2meV for Si:P, which is a factor of two larger than found in experiment [4], has very little effect on the fidelity.
Inhomogeneous broadening refers to a change in the excited state energies due to a differing electrostatic environment of the donor, e.g. caused by randomly distributed silicon isotopes. For the Si:P 1s-2p+-transition, this effect can result in a variation of excited state energies by about 0.1meV as measured in Ref. [4]. In our gate scheme, not taking into account this effect leads to an error in the detuning chosen for the pulse. In order check that this leads to insubstantial changes in the gate fidelity, we plot the infidelity for a large range of detunings around the optimal value in Fig. 5 for a fixed Rabi frequency Ω = 1000γ (corresponding to 2.8meV for the Si:P lifetime, such that a variation of 0.1 meV in the detuning corresponds to about 0.04Ω). The fidelities are not heavily affected, and we conclude that the gate does not need tuning to account for inhomogeneous broadening.

INTERACTIONS
With regards to the interactions between donors, we need to consider the tight binding model, with additional dipolar forces due to a donor being made up of a core plus  (13) where t is the single hopping matrix element between neighboring sites and where V is the Coulomb interaction, given by V (r) = V 0 /r where V 0 = e 2 4π 0 S , e is the electron charge and S = 11.4 is the dielectric constant of silicon.
The contribution from U ii ii = W ii is i =i W ii n ini wheren i = σ a † iσ a iσ which corresponds to the essentially classical Coulomb interactions between donors on neighbouring sites. From this we get: • Intersite Coulomb electron -electron repulsion where R = R 2 − R 1 is the separation vector between the two donors.
The contribution from U iiii corresponds to the on-site interaction, which can lead to anti-ferromagnetic superexchange which is negligible in our case.
The contributions from U ijji can be re-written, making use of Pauli matrix identities, into: (16) It corresponds to the ferromagnetic exchange coupling. Note that for calculating the total interaction for the gate, we are only interested in the case where both spins are on resonance with the laser, so both spins are aligned. This leads to − → S i · − → S j = 1 4n inj . • Ferromagnetic exchange coupling We refer to U ijji as J ij and we have: where R = R 2 − R 1 is the separation vector between the two donors.
The leading correction to this electron-only model of donor interactions is given by dipole interactions.
• Dipole-dipole interactions take the form: where p i is the dipole moment of atom i, R is the inter-atomic distance, and n is the unit vector between the two donors.
V Vdd is the dominant interaction in cold Rydberg atoms. At first order, the dipole-dipole interaction cancels, due to the vanishing dipole moment of donors. However, when an electric field is applied, dipole moments p i are induced and can be calculated as the derivative of the Stark shift of the state i which we study in App. .
A dipole can furthermore be induced by the vacuum, leading to Van der Waals forces, which are calculated in second order perturbation theory of V Vdd .
• Van der Waals forces: V VdW ij = − n ,l ,j ,m n ,l ,j ,m | n , l , j , m | ⊗ n , l , j , m | V Vdd |n, l, j, m ⊗ |n, l, j, m | 2 E |n ,l ,j ,m + E |n ,l ,j ,m − 2E |n,l,j,m To conclude, the Hamiltonian we use in this problem, in which we have checked that we can neglect the tunnelling at the distances we consider and we do not include electron-core interactions [6,7] as they cancel in the expression of the total interaction in the main text, corresponds to: In the presence of an electric field, W is dropped in favor of V Vdd which becomes the leading contribution.
In the case of the entangling gate, we are interested in the energy difference occurring when both donors are excited to the Rydberg state, compared to when none or only a single of the two is excited, as mentioned in the main text, which corresponds to E |rr , in which both spins are excited, E |rg in which one donor is excited and the other is in the ground state and E |gg in which both donors are in the ground state: where, again, V Vdd dominates W in the presence of an external electric field and we neglected the extremely small contributions coming from V VdW and J for |rg and |gg . FIG. 7. In the case of Si:P with |r chosen as 2p0, as is the case in the main text, the fidelity is robust to donors belonging to different silicon sublattices.
From these expressions, we are able to obtain the expression for the total interaction: J |gg , V VdW |gg , V VdW |rg , J |rg are negligible for the range of donor separations considered in this work. The interaction energy is plotted as a map for the 2p0 and 1sT 2 Γ 7 states in Se+ in Fig. 6.
For donor spin qubits in Si:P, most schemes make use of the ground state ferromagnetic exchange interactions, which oscillates wildly depending on the sublattice on which the donors reside. In Fig. 7 we show that in the gate proposed in this paper, the fidelity is not heavily affected by donor qubits residing on different sublattices.

FINITE ELEMENT METHOD MULTIVALLEY WAVEFUNCTIONS
In the following we explain how to obtain multivalley wavefunctions of donor electrons in silicon, which result from coupling the hydrogenic wavefunctions (envelope functions) into the manifolds determined by the symmetry group of the system. The envelope wavefunctions can be obtained by solving the Schrödinger equation for hydrogen, with the appropriate effective mass. We do this using the Finite Element Method [8] (developed in engineering for numerically computing approximations to solutions to partial differential equations and recently applied to donors [8]).
The effective mass approximation states that close to a band extrema, such as the conduction band minimum, each electron can be described by a mean-field Hamiltonian (every electron experiences the same average periodic potential) which is that of a single free electron with a modified mass (the effective mass) in an impurity potential.
Silicon has a face centred cubic lattice which leads to its brouillin zone (its primitive cell in reciprocal space) being a truncated octohedron. The latter contains six square faces, the centers of which are called X points. The conduction band minima in silicon are positioned at the X-points and are all equivalent. Using the KP method, the effective mass at a band extrema -such as the conduction band minima -can be determined. The effective mass corresponds to the dispersion of the energy E k as a function of momentum k (E k = 2 k 2 2m ). In silicon, the conduction band minima are anisotropic: they are ellipsoids which lie along the axis linking the center of the brillouin zone (Γ point) to the valley's X point. In valley specific coordinates, we set this to be the z axis. The effective mass along z is referred to as longitudinal m * l = 0.191m e , where as along x and y it is transverse m * t = 0.916m e , in reference to the conduction band minimum ellipsoid, where m is the mass of a free electron. The kinetic energy term in valley-specific coordinates is then: where γ = m t /m l . The impurity potential is the Coulomb potential of the proton, which at large distances can be approximated as a point charge at the nucleus, felt by the valence electron: where e is the electron charge. However, the valley degree of freedom in multivalley semiconductors allows the wavefunctions to have a non-negligible probability of being at the core where the Coulomb approximation breaks down. To account for this, we add another potential to the Hamiltonian, widely referred to as the central cell potential U cc (r), which takes the form of a delta function at the core and which we model using a Heavyside step function [9]. The value of the central cell potential is determined using a bisection algorithm by requiring the energy of the solution of the FEM method to match the experimentally determined energies for the multivalley states. It is also possible to add to the Hamiltonian the confining potential of an external electric field applied along the real-space z axis (which corresponds to the z-axis in valley-specific coordinates), which is: The Schrödinger equation which we will simplify in the following, in the valley-specific coordinates where z is the valley axis, is: (28) The wavefunction is written It is composed of the envelope function F i (r), where the indices i runs over the hydrogenic states, the Bloch wavefunction φ(k j , r) = e −ikj .r u kj (r) which is the product of a plane wave and a lattice periodic function u kj (r), where j runs over the six valleys, and lastly a multivalley parameter α i which couples the function into various manifolds according to the T d point group mentioned above. We use the finite element method to obtain F i (r) and couple it to the Bloch wavefunction and into the multivalley manifolds in a second step. As all the conduction band minima are equivalent, the solution to this Schrödinger equation suffices for all the manifolds representing the T d point group.
We write the Hamiltonian in the atomic units a B = 4π 0 R /e 2 m t and E H = e 4 m t /( 4π 0 R ) 2 : We make four steps to simplify the Hamiltonian to improve FEM convergence: • Transform to an anisotropic frame where z = √ z to have a symmetrised Hamiltonian.
• Transform to spherical polar coordinates (the kinetic energy has spherical symmetry so this is allowed).
• Choose the wavefunction where m is the magnetic quantum number, to obtain a 2D Hamiltonian.
Additionally, for wavefunctions with parity and magnetic quantum numbers of opposite parity y(η, θ) = −y(η, π − θ) therefore at the boundary: The final symmetrised polar compressed 2D Hamiltonian reads: Where H(r) is the Heavyside Theta function and Us(r) is the Unitstep function and r c and r m cutoff values. The solutions obtained via FEM are then numerically normalised using the Cuba Vegas package [10]. We now describe the coupling to the multivalley manifolds. Firstly, we reduce the six valley problem (−x, x, −y, y, −z, z) to an effective three valley problem (x, y, z) with x,y and z ∈ [−∞, ∞]. We associate to each valley the envelope function in rotated cartesian coordinates to ensure consistency when all the valleys will be coupled, considering we took z as the valley axis we create a three part vector corresponding to the envelope wavefunction in the three effective valleys: (F i (z, x, y), F i (y, z, x), F i (x, y, z)).
We transform the six functions of the tetrahedral point group into the effective three valleys by taking the Bloch wavefunction into account and using trigonometric identities to couple two valleys of opposite sign, each with a plane wave factor e −ikix into one with sin 2k i x if the manifolds have an even total sign or cos 2k i x if it is odd. For the 1sT 2 Γ 7 state for example: We then normalise the wavefunctions.

Ionisation rates
The ionisation rate for a hydrogen atom in an excited state is given in [12] to be: 3n 3 E f (35) with n, m, n 1 and n 2 the principal, magnetic, and parabolic quantum numbers respectively (in this paper we consider the 2p0 state which has n 1 = 1, n 2 = 0) and E f is the applied field strength in atomic units. We scale the results for the hydrogen atom by using experimental results acquired by [11] on Si:P in the ground state to determine the atomic units for the electric field used in Eq. 35. For an atom in the ground state Eq. 35 reduces originally given in [13], with α = 4 √ 2m t E 3/2 f /(3e ) the atomic field and ω = 12E b / , where E b is the binding energy, e is the electron mass and m t is the tunnelling mass. The tunnelling mass depends on the direction in which is applied the electric field, with the [111] direction yielding the slowest tunnelling times [11]. The 2D Hamiltonian which enables us to use the FEM for investigating induced dipolar interactions in this paper requires us to consider applying the electric field along [100]. The minimum ionisation rate is given by the Bohr period τ B = 2π /E b [11], on the order of 10 −13 seconds for Si:P and 10 −14 seconds for Si:Se+.
The results of these formula comply with experimental results for the Si:P 1sA state [11]. The classical ionisation threshold for Si:P 2p0 is at 0.28 V/um. In order to reduce the ionisation probability during the 2p0 lifetime of 235ps, we can see from Fig. 8 that the electric field should be on the order of 0.2 V/um.
The results for the theoretically calculated ionisation rates for Si:Se+ can be found in Fig. 8. The 2p0 result was obtained inserting the 1sA α and ω values into Eq. 35. The result for 1sT 2 (n 1 =0, n 2 =0) coincides with the result obtained when calculating α and ω values from the binding energy of 1sT 2 and directly using Eq. 36. In order to reduce the ionisation probability during the 1sT 2 experimentally measured lifetime of 7.7ns, we can see from Fig. 8 that the electric field should be on the order of 8 V/um for 1sT 2 and 1V/um for 2p0. In effect, these electric fields yeild negligible dipolar forces, as can be seen from

Stark shifts
In this paper, we propose a simple way to calculate the exact single valley Stark shift of donors using the FEM, by including a slope corresponding to the electric field in the Schrödinger equation in 28. This single valley Stark shift enables us to check our multivalley perturbatively calculated Stark shift because it includes contributions from the continuum. This is only possible for a field applied along the polarisation axis as, when applied perpendicular the Schrödinger equation no longer reduces to a 2D Hamiltonian.
Multivalley oscillations are on the same scale as the single valley wavefunction for the singly ionised donor Se+, in which the valence electron is tightly pulled in towards the core. All Stark shifts are plotted in Fig. 9, and are very close in magnitude for both methods -single valley or multivalley, for electric fields smaller than the ionising electric fields for the donor excited state lifetimes previously determined.
The perturbatively calculated single valley calculations include single s states with no central cell corrections (ccc) in the 2p0 calculations, and with a ccc for the 1sT 2 Γ 7 calculation in order to be able to compare directly to the FEM results. We can see both results are close, despite the perturbative results not including contributions from the continuum. The multivalley calculations include, for the ground as well as for higher lying s states, A, E and T manifolds, each with their own ccc.
For the 2p0 states of both donors, the main contributing state is 2s in the single valley case, and 2sE and 2sA in the multivalley case. For Se+ 1sT 2 Γ 7 in both single and multivalley cases, the main contributions to the Stark shift parallel to the polarisation of the state come from the p0 states and perpendicular to the polarisation comes from the p+ states.
Considering these Stark shifts, with laser linewidths below 0.1meV and Rabi frequencies on the order of 1meV, small electric fields would suffice to shift the donors on and off resonance. Large fields could even be applied to all the ground state donors not contributing to the gate to fully ensure they are off resonance with the laser exciting to the Rydberg state, because ground state donors can resist far larger fields than the excited states, which have a short lifetime and have higher ionisation probabilities.

Stray Electric Fields
In this section, we show that the stray electric fields from the gates local to the surrounding donors (Stark shifting them off resonance with the laser) do not interfere with the gate operation.
In order to reduce the probability of the laser exciting a transition, the detuning must be increased proportional to the Rabi frequency, for example for a 0.1 probability, ∆ = 4Ω. For realistic Rabi frequencies between 10 8 Hz and 10 11 Hz, depending on the experimental laser power, this sets an upper bound on ∆ of around 0.4 meV.
The 2p0 state in Si:P ionises in relatively small fields as can be seen in the previous section entitled 'Ionisation rates and Stark shifts'. A 0.4meV Stark shift would require an electric field of 0.5 V/µm, which we calculated using the multivalley perturbation theory method described in that section.
As is shown in the figure in the main text, we use the Finite Element Method to solve for the electric fields in the proposed geometry. We find that if the left local gates of donors not participating in the entangling operation are held at -0.0025V and the right ones at 0.0025V, with the gates surrounding the central donors all held at 0V, and the top global gate at -0.015V and the bottom one at 0.015 V, then the donors not participating in the entangling operation will get shifted off resonance with the laser with a field of 0.5V/µm, and the two donors participating in the readout, in the case of Si:P 2p0 being chosen as |r , will get a small field parallel to the interdonor axis due to the global gates, of around 0.18 V/µm, which will not ionise them within the excited state lifetime, but will induce dipolar interactions. We conclude that with this proposed geometry, there are no unwanted stray electric fields.