Optimisation of electrically-driven multi-donor quantum dot qubits

Multi-donor quantum dots have been at the forefront of recent progress in Si-based quantum computation. Among them, $2P:1P$ qubits have a built-in dipole moment, enabling all-electrical spin operation via hyperfine mediated electron dipole spin resonance (EDSR). The development of all-electrical multi-donor dot qubits requires a full understanding of their EDSR and coherence properties, incorporating multi-valley nature of their ground state. Here, by introducing a variational effective mass wave-function, we examine the impact of qubit geometry and nearby charge defects on the electrical operation and coherence of $2P:1P$ qubits. We report four outcomes: (i) The difference in the hyperfine interaction between the $2P$ and $1P$ sites enables fast EDSR, with $T_\pi \sim 10-50$ ns and a Rabi ratio $ (T_1/T_\pi) \sim 10^6$. We analyse qubits with the $2P:1P$ axis aligned along the [100], [110] and [111] crystal axes, finding that the fastest EDSR time $T_\pi$ occurs when the $2P:1P$ axis is $\parallel$[111], while the best Rabi ratio occurs when it is $\parallel$ [100]. This difference is attributed to the difference in the wave function overlap between $2P$ and $1P$ for different geometries. In contrast, the choice of $2P$ axis has no visible impact on qubit operation. (ii) Sensitivity to random telegraph noise due to nearby charge defects depends strongly on the location of the nearby defects with respect to the qubit. For certain orientations of defects random telegraph noise has an appreciable effect both on detuning and $2P-1P$ tunneling, with the latter inducing gate errors. (iii) The qubit is robust against $1/f$ noise provided it is operated away from the charge anticrossing. (iv) Entanglement via exchange is several orders of magnitude faster than dipole-dipole coupling. These findings pave the way towards fast, low-power, coherent and scalable donor dot-based quantum computing.


I. INTRODUCTION
Quantum computation architectures require long coherence times and a clear route towards scaling up. 1,2 Solid state spin qubits [3][4][5][6] are excellent candidates for large scale quantum computation 7 with outstanding coherence and fidelity, with Si:P donors 5,8-10 having a number of advantageous features. The strong Coulomb confinement potential of the donor atom comes for free and is reproducible, which, when coupled with extensive materials knowledge from the Si microfabrication industry, constitutes a viable avenue towards scalability. Thanks to the weak spin-orbit coupling 11 of Si electrons, the presence of hyperfine-free isotopes, 12,13 and absence of piezoelectric coupling to phonons, 14,15 the coherence time of Si:P donor electron spins is the longest among solid-state qubits. [16][17][18] Exceptional experimental progress has been registered in the past decade. 8,11,[18][19][20][21][22][23][24][25][26][27] Optimisation of the speed and scalability of multidonor quantum dot qubits will improve with our ability to operate them using only electric fields. Electric fields are much easier to apply and localise than magnetic fields, can be switched much faster, and consume less power. Spin qubits in single donors have exceptionally long coherence times but, lacking an intrinsic dipole moment, are challenging to operate electrically. Moreover, the multi-valley nature of the ground state (GS) [28][29][30] makes the exchange coupling 31,32 dependent on individual donor positions within the crystal. These oscillations have been mitigated by placing the donors along certain crystallographic directions 33 or by the presence of strain, 34 while the use of multi-donor quantum dots provides another approach to reduce this atomistic scale sensitivity by introducing valley-weight anisotropy. 33,35,36 Multi-donor quantum dot qubits can also be designed with a built-in dipole moment by making the charge densities different in adjacent dots, which enables electron dipole spin resonance (EDSR). 24,37,38 At the same time, dipole moments are known to expose the qubit to charge noise. The main questions at this stage in the development of multi-donor quantum dot qubits are: What is the largest achievable EDSR Rabi ratio in a donordot qubit, and can such an electrically operated qubit be robust against charge noise? These questions are theoretically challenging, since multi-donor dots are more difficult to model analytically and very expensive to treat computationally. [39][40][41] In addition, an understanding of EDSR and coherence properties in these multi-donor systems requires one to treat spatially and temporally random functions, which further increase the complexity and computational cost of the problem.
With these observations in mind, in this work we develop a variational effective mass wave function (EMA) wave function for 2P 22,42 and 2P : 1P 36 multi donor quantum dots in Si, and use it to study quantitatively the properties of a 2P : 1P double donor quantum dot spin qubit. It was recently shown that strong EDSR can be achieved in donor quantum dots using the hyperfine interaction, which plays the role of an effective spin-orbit field built into the qubit. 37, 38 We focus on the following properties: (i) the impact of qubit geometry on the EDSR gate time and Rabi ratio, and (ii) the impact of nearby charge defects, inducing random telegraph noise, on the coherence properties of the 2P : 1P spin qubit.
To determine the role of geometry in the electrical operation of the qubit, we consider three different orientations of the 2P:1P axis (the qubit axis): [100], [110] and [111], and compare the time scales relevant to qubit operation for these orientations. We find that the fastest EDSR time T π is achieved when the qubit axis is [111], while the largest Rabi ratio is found when the qubit axis is [100]. This can be explained by the difference in the wave function overlap between the 2P and 1P sites: both T π and T 1 /T π decrease linearly with t, while the tunneling t is highest for the 2P : 1P axis [100]. This is also the direction along which the 2P and 1P wave functions have the largest overlap, resulting in the slowest T π and largest Rabi Ratio. The tunnel coupling t is smallest for [111], hence T π is the fastest. In contrast, the orientation of the 2P axis only changes T π by 1-2%, and we conclude it has no visible effect on qubit operation. Since aligning the qubit axis [100] yields the largest Rabi ratio, in the latter part of the paper, devoted to coherence, we focus on this particular geometry. For the spin relaxation time T 1 , our theory yields 1/T 1 ∝ B 5 , consistent with an acoustic phonon mediated valley population mechanism. 43 This result in the low temperature limit is in accordance with earlier singleshot readout experiments. 8,44,45 Next we provide a quantitative analysis of the effect of charge defects on qubit operation. Importantly, a nearby charge defect can have a significant effect on both the detuning and 2P : 1P tunnelling, an aspect for which no quantitative studies exist in donor systems, although the issue has been known to exist in quantum dot qubits. 46,47 The contributions of the defect potential to detuning and tunnelling depend strongly on the defect location and orientation with respect to the qubit axis. As a result of the change in the 2P : 1P tunnel coupling, random telegraph noise can affect the EDSR gate fidelity in the vicinity of the anti-crossing. On the other hand, for 1/f noise the effect on the 2P : 1P tunnelling due to an ensemble of nearby defects is expected to be washed out, and operating the qubit away from the anti-crossing will drastically reduce its sensitivity to 1/f noise. In addition, our calculation shows that the anisotropic hyperfine interaction due to dipolar coupling between the electron and nuclear spins causes negligible decoherence. Finally, we examine entanglement via exchange 3,5 and dipole-dipole coupling 48,49 and find the former to be several orders of magnitude more efficient than the latter. Taken together, these findings pave the way towards high-performance all-electrical, coherent and scalable quantum computation using multi-donor quantum dots.
The outline of this paper is as follows. In section II, we present the variational wave function for a 2P donor-dot. Section III outlines the 2P : 1P qubit architecture, allelectrical qubit control, fidelity and coherence by studying EDSR, phonon relaxation and dephasing, and also explores the long-distance qubit couplings. We end with a brief summary.

II. VARIATIONAL 2P WAVE FUNCTION
The effective Hamiltonian for one electron in a double donor quantum dot is written as 50,51 where m * denotes the effective mass, including anisotropy within the quantum dot (the longitudinal effective mass m l =0.916m e and the transverse effective mass m t =0.191m e , with m e the bare electron mass), R L and R R stand for the position of the left and right donors in the 2P donor quantum dot respectively, and H v D represents the short-range part of the Coulomb potential of each donor, which gives rise to the valley-orbit coupling, discussed below.
The EMA wave functions for individual donors at position r are given by 16,52 where ξ is the index for the six valleys in Si, and |k ξ | = ±k 0 ≡ ±0.85 (2π/a Si ), with a Si = 5.43Å the lattice constant of Si. The latticeperiodic function is u ξ (r) = K c ξ K e iK·r , with K reciprocal lattice vectors. Here φ D,ξ (r − R D ) denotes the hydrogenic part of the wave function. The resultant states 53 contributing to the envelope of the donor wave function 54 take the form of a deformed hydrogenic 1s orbital with two radii a and b, φ ±y (r) and φ ±z (r) can be obtained from Eqn. 2 by interchanging y ↔ x and z ↔ x repectively. The singledonor variational wave function is given by ξ φ D,ξ (r)e ik ξ ·r u ξ (r). The valley-orbit coupling 55-59 is added to a single donor through the term H v D = U 0 δ(r − R D ). In the basis spanned by the six valley wave functions {φ D,ξ (r)e ik ξ ·r u ξ (r)} this will yield the standard Kohn-Luttinger matrix elements ∆ , ∆ ⊥ and δ, and we choose the value of U 0 to reproduce the correct energy splitting between the ground and first excited states. 60 This enables us to circumvent complications associated with the central cell correction r cc : a single r cc cannot be chosen to produce both the donor binding energy and the valley-orbit coupling correctly. 30  (1,0) in the x-y plane, with secondary maximum due to short-range traveling waves at (−1, −2), (1,2) and so on.  For a 2P quantum dot, the wave functions for the left and right donors are given by φ L,ξ (r) = φ ξ (r − R L ) and φ R,ξ (r) = φ ξ (r−R R ) respectively. We determine the dot wave function based on the variational method developed in Ref. 62 for the hydrogen molecule, then add the valleyorbit coupling due to the two donors as above. The variational method has been shown to reproduce the spectra of H 2 and He with great accuracy. 62 We treat a and b as variational parameters 62 which are a function of donor distance, and define |F D,ξ = φ D,ξ e ik ξ ·(r−R D ) . The matrix elements of the Hamiltonian H ±ξ = ε ξtξ t ξ ε ξ between the wave functions for the same valley {F L,ξ , F R,ξ } are the on-site energies ε ξ = φ D,ξ |H 2d |φ D,ξ and the tunneling energiest ξ = φ L/R,ξ |H 2d |φ R/L,ξ , with ε ξ ,t ξ real. These matrices are diagonalised by symmetric and anti-symmetric combinations |S ξ , |A ξ , which provide an orthonormal basis of valley wave functions. Since thẽ t ξ 's are negative, the symmetric functions |S ξ have a lower eigen-energy E ξ = ε ξ −|t ξ | than the anti-symmetric functions |A ξ with E ξ = ε ξ + |t ξ |. We can ignore the anti-symmetric wave function as the inter-valley matrix elements are small compared to the energy difference of 2t ξ between the symmetric and anti-symmetric functions. We are left with a 6 × 6 matrix in the manifold spanned by |S ξ . Minimization of S ξ |H|S ξ yields the variational parameters a and b (Table I) 63,64 Comparison of the ground state energy, valley composition and exchange oscillations between our method and much more advanced computational approaches such as finite-element method 61 , Nano-Electronic MOdelling 65 is surprisingly encouraging. Figure 1a depicts the color coded amplitude of the 2P quantum dot wavefunction 61 in the x-y plane with the single donors situated at the closest lattice spacing of a Si = 0.54 nm along the [100] direction, using the envelope function and incorporating the phase factor e iK·r that determines the short-range structure of Ψ 2D . In contrast the 2P [110] wavefunction is shown in Fig. 1b, where the closest spacing between the two donors is √ 2a Si = 0.76 nm. Fig. 1c shows the diminished secondary maximum of the 2P [111] wavefunction in the x-y plane, with both donors outof-plane. The ground-state wave function simplifies to Ψ 2D (r) = 1 √ ξ w 2 ξ ξ w ξ S ξ u ξ (r), with w ξ representing the valley weight. The ground state of the double-donor quantum dot is symmetric under inversion yet, unlike the case of a single-donor quantum dot, the wave function is no longer spherically symmetric. Importantly this spatial anisotropy has a significant effect on the valley composition of the ground state. The large effective mass anisotropy (EMA) results in a large enhancement of the kinetic energies of the valley states lying perpendicular to the orientation of the 2P quantum dot. When the 2P axis is [100] the contribution to the kinetic energy stemming from the x-valleys will be smaller by a factor of 4.8 than that due to the y and z valleys perpendicular to the 2P axis. In this case the ground state would comprise only the y and z valleys. The valley-orbit coupling (VOC) arising from short-range Coulomb potential H v D (Eqn. 1) on the other hand is responsible for matrix elements of approximately the same magnitude connecting all valleys, and therefore favours a ground-state superposition of all valleys with comparable weight. The competition between the effects of effective mass anisotropy and valley-orbit coupling ultimately determine the ground state valley weight ratio.
• When the 2P axis [100] the valley weight ratio w ±x : w ±y : w ±z = 0.88 : 1 : 1, indicating a higher contribution to Ψ 2D from the y and z valley states. • When 2P axis is [110] the valley weight ratio w ±x : w ±y : w ±z = 0.91 : 0.91 : 1 indicating equal and lower contribution from x and y valleys.
• When the 2P axis is [111] all valley weights are equal, as the effects of EMA and VOC are symmetric in all valleys. To summarise this section, we have presented a simple and practical variational wave function describing the ground state of a 2P donor dot, whose structure reflects the interplay between the effective mass anisotropy and valley-orbit coupling, to yield an accurate value for the ground state energy. This wave function is very convenient in formulating a simple analytical model for the operations of a 2P : 1P qubit, as the next section shows.

III. 2P : 1P QUBIT
In order to operate a multi-donor dot qubit electrically we require a dipole moment, which a 2P donor dot does not have. On the other hand, a double donor dot qubit in a 2P:1P configuration has a difference in the charge densities between the 2P and 1P sites which gives rise to a dipole moment that can couple to an applied electric field. In what follows we concentrate on such a 2P : 1P qubit and discuss the electrical operation and EDSR using the hyperfine coupling to the nuclei. Next, considering the fact that this dipole moment also couples to phonons and charge fluctuations, we discuss the coherence properties of an electrically operated 2P : 1P qubit. Finally, we discuss briefly the prospects for entangling two 2P : 1P qubits using the exchange interaction as well as the dipole-dipole interaction.

A. Effective Hamiltonian for a 2P : 1P qubit
We first introduce the general formalism describing a 2P : 1P qubit. We adopt the Hund-Mulliken approach, which takes into account the overlap of the 2P and 1P wave functions. Fig. 2a represents a schematic of the STM device architecture showing the proposed qubit system, in which the 2P axis and the 2P : 1P axis are both [100]. The 2P dot is situated on the left and the 1P single donor quantum dot is on the right. The ground states are orthonormalised as: we include the on-site energies, inter-dot tunnelling, and Coulomb interaction effects as described in Ref. 67. The 2P : 1P exchange energy is given by the difference between the two lowest eigenvalues of the 4 × 4 matrix spanned by the two-particle basis mentioned above (See supplementary material).
A detuning dc electric field is applied between the inplane gates M , L and R to drive the 2P : 1P qubit to the charge anti-crossing. The detuning, denoted by δ, represents the corresponding energy difference between the 2P and 1P sites. The resultant ground and excited orbital state energies are e 1 = 1 2 (e l + e r + δ − δε), e 2 = 1 2 (e l + e r + δ + δε). The anti-crossing occurs at the point where e 2 − e 1 = δε = (e l − e r + δ) 2 + 4t 2 → 2t; e l , e r being the on-site energies of 2P (left) and 1P (right), t signifies the tunneling energy. The hyperfine interaction between the electron and effective nuclear field from the three donors is with A 0 = (2/3)µ 0 γ e γ n ; γ e , γ n being the gyromagnetic ratios of the electron and proton respectively; − → S denote nuclear and electron spin operators and i = 1P, 2P L , 2P R label the locations of the nuclear spins. 68 We assume isotopically purified 28 Si so that there are no hyperfine fields other than those due to the donors given above. An external magnetic field B applied along z resolves the orbital states by Zeeman energy ±e z to produce electron spin-up and spin-down states. 37, 69 We include a driving ac electric field E(t) applied between M , L and R (Fig. 2b), which gives rise to an oscillating electrical potential v ac (t) = e E(t)x. The total Hamiltonian H q in the 2P : 1P orbital+spin basis {G ↑, G ↓,E ↑, E ↓} reads: The matrix elements the expectation values of the nuclear spin angular momentum operators. k gg b z stems from the out-of-plane hyperfine interaction, and this is simply added to the total Zeeman splitting e z , modifying the spin states splitting as z = e z + k gg b z . In summary, we have determined a simple, exactly diagonalisable analytical Hamiltonian that describes the electrical operation of a 2P : 1P qubit. This will be the topic of the next section.

B. Hyperfine mediated EDSR and relaxation
In the vicinity of the charge anti-crossing, the hyperfine interaction enables an electron spin-flip transition (1, 0) ↑, ↓↔ (0, 1) ↓, ↑ between the left and right dots. If the ac electric field is in resonance with the qubit Zeeman splitting, an electron spin flip takes place owing to the time dependent modulation of the difference between the in-plane hyperfine interactions on the left and right dots, a mechanism analogous to that observed using a micromagnet, 15,25,[70][71][72] termed hyperfine mediated electron dipole spin resonance (EDSR). 23 To describe the effect of a driving AC electric field in the plane of the qubit, the Schrieffer-Wolff (SW) transformation is applied to Eqn. 4 after rotating the qubit basis, which enables us to obtain a low-energy, 2×2 effective qubit Hamiltonian (See supplementary material). The effective Zeeman splitting ζ due to an applied external magnetic field and out-ofplane hyperfine coupling is much smaller than the orbital energy splitting δε = √ 2 + 4t 2 , where = e l − e r + δ and t denotes 2P − 1P tunneling. Using 2ζ δε, a condition which boils down to ζ t at the charge anticrossing, we expand up to second order in small terms. The EDSR rate for a π-rotation is given by the off-diagonal termH 12 between G ⇑ and G ⇓ as follows, The same matrix elementH 12 enables relaxation via coupling to phonons, 11,73 as described in the supplementary material. The EDSR rate is independent of the qubit Zeeman splitting, a fundamental difference between EDSR mechanisms based on the hyperfine and intrinsic spin-orbit interactions. 74,75 For the results in Fig. 3 the 2P : 1P qubit is operated at the charge anticrossing where the EDSR time is T π = 2t vacbkge . Following Eqn. 5. T π ∝ t implies that at larger 2P-1P separations EDSR is faster, as 2P-1P tunneling decreases. Nevertheless our SW approximation breaks down beyond a certain value of the 2P:1P separation where 2P-1P tunneling is very small and the condition ζ t is no longer satisfied.
The nuclear spins are polarised in the direction of the external magnetic field i.e. the z-direction; and in an ensemble-averaged experiment the in-plane b + , b − would average to zero leading to a washing out of the Rabi oscillations. 23 To have non-zero transverse components one needs to apply a π 2 pulse to the 1P nuclear spin. We consider three possible orientations of the 2P:1P qubit axis,-along [110], [110] and [111] respectively. We have found that changing the direction of the 2P axis alters the results by 1−2%, which in practice has no visible effect on the qubit EDSR and relaxation. As the valley compositions of 2P 100, 2P 110 and 2P 111 are not drastically different, for fixed orientation of the 1P donor dot the overall 2P-1P tunneling matrix does not change much. On the other hand, rotating the 2P:1P qubit axis means the overlap between same valleys (e.g. 2P x to 1P x , 2P y to 1P y etc.) changes considerably, which produces significant changes to the overall tunneling matrix elements. Therefore, for simplicity, in the three 2P-1P configurations considered here we take the 2P axis to be oriented along the same direction as the 2P : 1P qubit axis (see inset to Fig. 3a). Fig. 3 shows the effect of varying the orientation of the 2P : 1P qubit axis, as well as the effect of changing 2P − 1P separation on the EDSR time T π and Rabi ratio of the relaxation time T 1 to the gate time T π , all evaluated at the charge anticrossing. The number of gate operations per relaxation time decreases as the dots are further apart. We consider first having the 2P : 1P axis [100]. For a constant 2P-1P separation of 13 nm, the 2P-1P tunnel coupling takes the value t = 1.7 meV, with valley contributions t ±x = 0.31 meV, t ±y = 0.26 meV, t ±z = 0.26meV . The magnitude of the 2P − 1P tunnel coupling is at its highest for 2P:1P [100]. We find that the x-valleys make the dominant contribution to t (Fig. 3a). Along this direction the wave function overlap between 2P and 1P , and with it the tunnel coupling, is maximal. Since T π and T 1 /T π are both ∝ t, this means the EDSR gate time is longest along this direction, but the Rabi ratio is also at its largest. Next, having the 2P : 1P axis [110] results in the x-and y-valleys making a higher contribution to tunneling than the z-valleys, namely t ±x = 0.20 meV, t ±y = 0.20 meV, t ±z =0.11meV , yet the total value of the tunnel coupling decreases as compared to [100], t = 1 meV. This is because the overlap of the 2P and 1P wave functions is smaller than it is along [100]. Hence the Rabi ratio T 1 /T π as well as the EDSR time T π have smaller values along [110] than along [100]. When the qubit axis is [111] all valleys contribute equally: t ±x = 0.08 meV, t ±y = 0.08 meV, t ±z =0.08meV , with the total 2P-1P tunneling t = 0.5 meV, lowest among the three orientations; so the Rabi ratio further decreases ( Fig. 3b)  We have thus established that a 2P : 1P qubit can exhibit fast hyperfine-mediated EDSR and over a million electrical operations within one relaxation time, and that, owing to the difference in interdot tunnelling for the three orientations investigated, the performance is best when the qubit axis is [100]. 2, 1/f ) due to 1/f noise, as a function of detuning (δ). Formally T * 2,1/f tends to infinity here since we have not included higher order terms in the defect potential.

C. Nuclear flip-flops and the anisotropic hyperfine interaction
We now consider the effect of the nuclear field of the three donor nuclear spins (2P+1P) on the electrical operation of the qubit. To begin with, the nuclear dipole-dipole interaction 71 could induce nuclear-nuclear flip-flops between the three donors independently of the applied electric field. Although the nuclei are relatively close (the two nuclei on the 2P dot are 0.5 nm apart, whereas the nuclei on 1P is situated 12 nm apart from nuclei on 2P ), the time scale for nuclear flip-flops is expected to be 10 ms, 71 which is very slow compared to the electric field frequency since the rate depends on the product of two nuclear g-factors.
We address potential feedback effects on the nucleus in the course of EDSR. One possibility is for the nuclear spin to experience electrically driven spin resonance, in which the nuclear spin is flipped by the combination of the ac electric field and inhomogeneous Zeeman field between 2P -1P due to the electron spin, which enters in the hyperfine interaction. This possibility can be safely discounted, since for it to occur the hyperfine interaction with the electron spin would have to couple the nuclear ground state to orbital excited states of the nucleus. These excited states are extremely high in energy, of the order of MeV. 76 Moreover, the ac electric field frequency will match the Zeeman splitting of the electron spin states, and will be detuned by three orders of magnitude from the nuclear spin state splitting due to the difference in g-factors.
The electron spin may also have a back-action on the nuclear spin, since the electron spin is rotated by EDSR, and this rotation could in principle affect the hyperfine interaction. Flipping the electron spin could result in the hyperfine interaction being also flipped. Nevertheless, as long as the external magnetic field is large enough, this sign flip should have a minimal influence on the nuclear spin. It is worth emphasising that the electron-nuclear system is not allowed to evolve at its natural frequency. Rather, the coefficients A(R i ) are modified by the electric field fast enough that the electron effectively senses an inhomogeneous magnetic field that depends on its position, analogous to a micromagnet. 77 The anisotropic hyperfine (AHF) interaction between the electron and nuclear spin arising from magnetic dipolar coupling 78 can lead to fluctuations in T π . Compared to the contact hyperfine Hamiltonian (Eqn. 3), the AHF Hamiltonian 79 3 γ e γ n 2 3xixj −r 2 δij r 5 . The contribution of AHF to the EDSR matrix element is 0.007 times the contact hyperfine at the anticrossing, 72 which will have a negligible effect on qubit operation.
We conclude that nuclear flip-flops and the anisotropic hyperfine interaction do not play an important part in the electrical operation of 2P : 1P qubits.

D. Charge noise effects on qubit operation: decoherence and gate errors
We assume the sample is isotopically purified, so that there is no random hyperfine field contributing to qubit decoherence. Nevertheless, as compared with singledonor quantum dot qubits, the 2P : 1P system itself has a dipole moment and is therefore exposed to charge noise. We now discuss qubit decoherence in the presence of this charge noise, which causes certain quantities to fluctuate in the effective 2×2 qubit Hamiltonian. 80 The presence of charge defects in the qubit environment can affect both the 2P -1P tunnelling energy, producing a fluctuation ∆t in the tunnelling, as well as in the charge distribution between 2P and 1P leading to a fluctuation ∆v in the detuning. 8,26,81 We focus first on random telegraph noise (RTN), considering a charge defect in the vicinity of the physical qubit, and varying its location. 82 The setup is illustrated in Fig. 4a. The defect is characterised by one switching time τ , which we take to be τ = 1µs, since experimentally the effect of fluctuators with switching times longer than this can be eliminated by means of dynamical decoupling. For the range of defect distances studied here, the fluctuations ∆t, ∆v satisfy ∆t, ∆v ( /τ ), hence, as described in Ref. 82, the dephasing rate can be expressed as where ζ is the effective Zeeman splitting, k ge is a matrix element of the hyperfine interaction introduced above, the orbital gap is δε, and = (e l − e r + δ). The RTN dephasing time is determined by fluctuations in both the tunnelling and the detuning, and in what follows we study their relative contributions for different defect locations. Fig. 4b shows that ∆t, the fluctuation in the tunnel coupling, is a strong function of the angle φ, which is the azimuthal angle characterising the position vector of the charge defect. It is symmetric about π, while its absolute value reaches a maximum at φ = π/4, 7π/4. Interestingly, the integral of ∆t over the angle φ is zero, which implies that, if a significant number of defects were randomly located in the vicinity of the qubit, their effect on tunnelling is expected to be washed out. On the other hand, the fluctuation ∆v in the detuning is nearly constant for all φ. Finally, the variation of ∆t and ∆v with charge defect distance R is shown in Fig. 4c, showing that, as expected, they both decrease with increasing defect-qubit separation at approximately the same rate. The dependence of ∆t and ∆v on the vertical position of the charge defect is shown in panel 4d. Whereas there is little overall variation in the range studied, the absolute values of both ∆t and ∆v reach their maxima at h = 0, R = 30 nm, φ = ±40 o , i.e. the worst position for a defect to impact T * 2 is in the same plane as the qubit, at an angle ±40 o to the inter-donor axis and at a distance 30 nm away. A charge defect is shown therefore to have a more detrimental effect when placed closer to the 1P dot, and, when it lies in the same plane as the qubit.
At the anticrossing, where = 0, the term ∝ ∆v does not contribute to T * 2,RT N , the only contribution stems from the tunneling fluctuation ∆t. Away from the anticrossing the RTN dephasing time increases because the orbital energy gap δε increases (Fig. 4e). Both the tunnelling fluctuation ∆t and detuning fluctuation ∆v contribute to dephasing, and T * 2,RT N ∝ δε 8 . The potential of a charge defect makes a substantial contribution to the off-diagonal terms in the qubit Hamiltonian, which are responsible for EDSR. This can lead to EDSR gate errors. To quantify this effect we consider  between the fluctuation-induced term in the Hamiltonian and the EDSR Rabi frequency. At the anticrossing, where only tunneling fluctuations ∆t contribute, the maximum value of this ratio in the range studied is 0.57 MHz with a driving electric field of 10 kV /m, which corresponds to 2% of the Rabi frequency for our qubit setup. This suggests gate errors could be an important effect of random telegraph noise. At the same time, the fractional change in the Rabi frequency could be reduced by simply increasing the driving electric field.
For 1/f noise, which arises from an incoherent superposition of a large number of random telegraph sources, we will assume that the effect of tunneling fluctuations ∆t is negligible, using the reasoning presented above. On the other hand the detuning is roughly constant as the charge defect location is rotated in the plane so the detuning fluctuations are expected to be important at all times. In light of this, the 1/f noise spectrum 27 S(ω) = Ak B T /ω, in the qubit subspace, takes the form where k ge is hyperfine energy, = e l − e r + δ; e l , e r are the on-site energies of 2P and 1P dots, δ is the detuning field, ζ signify qubit Zeeman splitting and δε = √ 2 + 4t 2 is orbital energy splitting respectively. A = 0.1 µeV is estimated from experiments. Using S x (t) = S 0x e −χ(t) (See supplementary material), we obtain: At the anti-crossing the numerator in Eqn. 7 is zero in this approximation as the detuning fluctuation does not contribute, leading to a minimum in dephasing rate. Clearly, when higher-order terms in the defect potential are considered, the value at the anti-crossing will be a small but finite number. Close to the anti-crossing the numerator in Eqn. 7 is dominant resulting in two peaks in the dephasing rate (T * 2 ) −1 symmetric about the anticrossing, while further away from the anti-crossing the denominator, which is of higher order in the detuning, dominates the decoherence (Fig. 4f). In other words, (T * 2 ) 1/f decreases close to the 2P-1P charge degeneracy point, while increasing further away from anti-crossing. A direct comparison of the coherence properties with single-spin qubits is impossible in the absence of noise data for 2P : 1P (even for 1P ). Nevertheless, putting together all the above findings suggests the safest operational regime for the 2P : 1P qubit is away from the anti-crossing, where sensitivity to RTN dephasing and gate errors, as well as to 1/f noise, is minimised.
E. Noise sensitivity in 2P : 1P donor dot qubits vs electrically operated quantum dot qubits We discuss briefly the sensitivity to noise in 2P : 1P donor dot qubits as compared to electrically operated quantum dot qubits. A dipole moment can be induced in quantum dots either by the spin-orbit interaction 85 or by a nearby micromagnet. 86,87 Dephasing in these architectures was considered in Refs. 83 and 84 respectively. The sensitivity to noise is determined by two factors: (i) the magnitude of the spin-mixing interaction, which for quantum dots is either the spin-orbit coupling of the magnetic field gradient, while for a 2P : 1P qubit it is the difference in the hyperfine interaction between the 2P and 1P sites; and (ii) the asymmetry of the charge distribution, which in a quantum dot is the asymmetry between the ground state and first excited state wave functions, whereas for a 2P : 1P qubit it is the asymmetry between the 2P and 1P wave functions.
For single-spin qubits in QDs noise introduces a fluctuation between the charge distributions in the ground state and first excited state respectively, leading to a fluctuation δv using the notation in Ref. 83. A complete comparison is difficult to make given that dephasing in the quantum dot depends on a number of parameters including the dot radius and aspect ratio. Nevertheless, one comparison can be made straightforwardly: setting the Rashba spin-orbit interaction in a quantum dot qubit equal to the to the hyperfine interaction term k ge in the 2P : 1P qubit (which is a very realistic assumption), setting the Zeeman splittings to be equal in both qubits, and considering a defect a fixed distance away from either qubit and with the same switching time, we evaluate the resulting T * 2 for RTN and 1/f noise. Table. II establishes that noise properties are better for 2P:1P donor configuration than single-spin qubit in QDs, given that the multi-donor dot qubit is operated away from the anticrossing. The qubit orbital energy gap δε increases in the operational regime away from anticrossing, also the spin-orbit interaction coming from hyperfine decreases, resulting in longer dephasing times.
The orbital splitting dependence for 2P:1P is given by T * 2,1/f ∝ δε 4 similar to the micromagnet mediated QD qubit in Ref. 84. A somewhat different trend is observed for single-spin QD qubits in Ref. 83, where electrical operation relies on the spin-orbit interaction, resulting in T * 2,1/f ∝ δε 3 . This discrepancy is accounted for by the different energy dependencies between the matrix elements of the spin-orbit interaction and those of the inhomogeneous magnetic field of the micromagnet.
In summary, 2P : 1P qubits are more robust against noise than the equivalent electron quantum dot architectures.

F. Entanglement: comparison of exchange and dipole-dipole coupling
To conclude our discussion, we focus briefly on two possibilities for entangling 2P : 1P qubits: exchange coupling 3,5 and dipole-dipole coupling. 48,49 We consider first the exchange energy in the Hund-Mulliken (HM) approximation for two 2P:1P multi-donor quantum dot qubits. Each qubit is oriented along [110], and the qubits are arranged head-to-tail, i.e. perpendicular to [110], as illustrated in Fig. 5. The ground state for each qubit is denoted by G cf. Eqn. 4. Variation of the exchange energy J 2Q with inter-qubit distance d reveals the absence of exchange oscillations when the qubits are operated near the anticrossing (Fig. 5). Next we compare the two-qubit exchange coupling with the two-qubit dipoledipole interaction. Exchange is negligible compared to the direct Coulomb interaction when the distance between the two 2P : 1P qubits is ∼ 20 nm. Using the basis of two-particle product states, we can calculate the 16 × 16 matrix consisting of the direct Coulomb interaction terms 48,49,[88][89][90] Then we use the Schrieffer-Wolff (SW) perturbative approach to transform the interaction into the 4 × 4 two-qubit spin states manifold. Keeping only the dipole terms, the third order SW produces the Ising interaction H   model to earlier results in terms of operation speed T π , relaxation time T 1 ; and also presents our two-qubit gate speeds. In summary, entangling 2P : 1P qubits using Coulomb exchange is several orders of magnitude faster than using the dipole-dipole interaction.

IV. SUMMARY
We have developed an analytical wave function for a 2P donor based quantum dot to construct a model of a 2P : 1P qubit, which is found to be a highly suitable candidate for all-electrical spin quantum computing. Fast EDSR can be achieved in a 2P:1P multi-donor quantum dot with gate times of 10−50ns at electric fields of 10000 − 50000 V/m at the charge anti-crossing. The spin relaxation time due to phonons satisifes 1/T 1 ∝ B 5 and, at the anti-crossing, allows in excess of a million qubit operations. Random telegraph noise can lead to sizable fluctuations in tunneling between the 2P and 1P sites as well as fluctuations in the detuning, which can cause both decoherence and gate errors leading to loss of fidelity. Nevertheless, we have shown that qubits can be immune to RTN and 1/f noise some distance away from the charge anti-crossing. Efficient entanglement can be achieved using exchange, which does not exhibit oscillations as a function of qubit separation, while the dipoledipole interaction is considerably slower. In the future our theoretical method could be extended to 3P multidonor quantum dot configurations and to several qubits in order to examine cross-talk and further issues related to scaling up multi-donor quantum dot qubits.

Acknowledgments
The authors would like to thank Andre Saraiva for many enlightening discussions.
where the 1P donor dot is situated at X and the 2P multi donor dot is situated at −X. S is the overlap integral and is given by: The two particle singlet-triplet states are: where Ψ s− Ψ s+ Ψ s0 are singlet states, symmetric in r 1 ↔ r 2 ; Ψ t is the asymmetric triplet state. The orbital Hamiltonian in the basis with, h tot =h l + h r + C, where e l is the on-site energy of Ψ − (r). U 1 denotes on-site repulsion.
The exchange energy is given by the the difference between the two lowest eigenvalues of the 4x4 orbital Hamiltonian in Eqn. 16.

II. EFFECTIVE EDSR HAMILTONIAN
stems from the out-of-plane hyperfine interaction, and this is simply added to the total Zeeman splitting e z , modifying the spin states splitting as z = e z + k gg b z . First, the upper diagonal 2 × 2 block of the H hf is considered.
Diagonalisation of H 2×2 hf gives us a new set of basis states with spin up and down combination, i.e. a rotated basis: The hyperfine strength A 0 depends on the detuning δ, which can be checked by varying δ in the qubit ground state G and excited state E wavefunctions. Near the anticrossing, A 0 has a linear dependence on detuning and allows the charge transition (1, 0) ↔ (0, 1) between left and right dot within 10kV /m electric field range.An electric dipole transition can be driven by an ac electric field E(t) applied between G1 and G2, modifying A 0 with time. If the ac electric field is in resonance with the qubit Zeeman splitting, a spin flip takes place due to the difference between the hyperfine energy of the left and right dot, a mechanism that has also been observed using a micromagnet. This is known as the hyperfine mediated electron dipole spin resonance (EDSR). Turning on the ac field, the EDSR Hamiltonian in the rotated basis is given by: where, ζ = 2 z + (b 2 x + b 2 y )k 2 gg . The matrix elements are calculated below: (α 2 + γ 2 ) = 1 can be deduced straightforward. Similarly, Evidently, H 23 = Q 2 = Q 1 . k ge b z /ζ 1 400 , so the terms P 1 , P 2 , Q 1 , Q 2 , to very good approximation, effectively become: The off-diagonal spin-flip term after applying the Schrieffer-Wolff transformation, under the assumption that the effective qubit Zeeman splitting is small compared to the orbital energy gap, ζ t: The diagonal elements from Schrieffer-Wolff transformation are given by, The 2x2 effective qubit Hamiltonian is given by H ij ; i, j ∈ {1, 2}.

III. RELAXATION TIME
The relaxation time T 1 carries information about the coupling of the electron spin with the crystalline and electromagnetic environment. The anisotropy of the deformation potential is taken into account, and the Hamiltonian takes the form: H e−ph = λ q −iq(a † qs exp(iq · r) + a qs exp(−iq · r)) × (Ξ d e λ,xqx + Ξ d e λ,yqy + (Ξ d + Ξ u )e λ,zqz ) where a † qs and a qs are the phonon creation and annihilation operators, q is the phonon wavevector,ê λ is the phonon polarization vector, the deformation potential is given by the tensor Ξ ij : The initial and final phonon number states produces n q+1 |a † qs |n q ∝ n q + 1; in the low temperature region n q + 1 ∼ 1, as n q follows Boltzmann distribution. The electron-phonon Hamiltonian can be written as (Ξ d e λ,xqx + Ξ d e λ,yqy + (Ξ d + Ξ u )e λ,zqz ) (34) where N V denotes the volume of the crystal, ρ is the mass density and ω qλ is the acoustic phonon frequency. Using ω qλ = v λ q, we obtain: The perturbation in our basis states due to presence of vibration would modify the qubit ground state |G = |G + E|H e−ph |G |E e 2 − e 1 the qubit spin states couling term becomes: G ↑ |H hf |G ↓ = 2 E|H e−ph |G * e 2 − e 1 E ↑ |H hf |G ↓ 1 where Ω denotes the angular integral. G|iq·r|E G|1+iq·r|E =¨¨B 0 G|E + G|iq·r|E ; using E ↑ |H hf |G ↓ ≡ |b|k ge , and assuming |b| = 1: B = 1 T out-of-plane magnetic field is used. The dependence of T 1 /T π and T 1 (inset) on qubit geometry comes from k ge , t, E|q · r|G terms.
IV. 1/f NOISE 1/f noise is Gaussian for semiconductors and can be described by the correlation function S(t) = V (0)V (t) . The Fourier transform of this noise correlation function varies inversely proportional to frequency, S(ω) = Ak B T ω , with A = 0.1 µeV is estimated from experiments. [2] The time evolution of spin component is given by, Using the noise correlation function in the qubit subspace, we can write: The main contribution to 1/f noise comes from the low-frequencies, so a low frequency cut-off ω 0 is chosen as the inverse of measurement time. Assuming, ω 0 t 1, and approximating ∞ ∞ dω sin 2 (ωt/2)