Full configuration interaction simulations of exchange-coupled donors in silicon using multi-valley effective mass theory

Donor spin in silicon have achieved record values of coherence times and single-qubit gate fidelities. The next stage of development involves demonstrating high-fidelity two-qubit logic gates, where the most natural coupling is the exchange interaction. To aid the efficient design of scalable donor-based quantum processors, we model the two-electron wave function using a full configuration interaction method within a multi-valley effective mass theory. We exploit the high computational efficiency of our code to investigate the exchange interaction, valley population, and electron densities for two phosphorus donors in a wide range of lattice positions, orientations, and as a function of applied electric fields. The outcomes are visualized with interactive images where donor positions can be swept while watching the valley and orbital components evolve accordingly. Our results provide a physically intuitive and quantitatively accurate understanding of the placement and tuning criteria necessary to achieve high-fidelity two-qubit gates with donors in silicon.


I. INTRODUCTION
Donor spins in silicon are prominent candidates for the construction of a large-scale quantum computer. The spin-1/2 degree of freedom of an electron bound to a donor atom is a natural qubit. Isolated in the environment of isotopically purified 28 Si, it has demonstrated exceptional coherence times. 1 High fidelity initialization, readout, 2 and single-qubit operations [3][4][5] are well established. Several proposals for multi-qubit gates rely on the exchange interaction, using either strong tunable exchange to perform SWAP oscillations, [6][7][8] or weak exchange and microwave pulses to perform a CROT gate. 8,9 Experimental work by He et al. 10 demonstrated fast SWAP oscillations in the strong exchange regime, using two donor clusters fabricated with a scanning tunnel microscope. Madzik et al. 11 showed an embryonic CROT gate as the coherent rotation of a target qubit, conditional an a specific state of the control, using a pair of ion-implanted 31 P donors, in the weak exchange regime.
Moving from these proof-of-principle demonstrations to systematic achievement of high-fidelity two-qubit gates will be greatly assisted by having an accurate and efficient tool to simulate the electronic structure of coupled donors. Here, we achieve this goal by developing and applying a full configuration interaction implementation of effective mass theory that avoids unjustified approximations made in some earlier theories. Our code can be made computationally very efficient, which allows us to study detailed properties of coupled donors as a function of a very wide range of parameters. This also allows us to pedagogically unpack the microscopic physics that determines the results. We trust that the physically transparent nature of this work will help future researchers gaining better insights into the subtle physics of donor spin qubits.
Exchange-coupled donor molecules (Fig. 1a) are many-electron systems and computationally very complex. Their study requires the application of quantum chemical methods and approximations 30,31 to the multi-valley effective mass theory frame work. Widely used is the Heitler-London approach. [32][33][34][35] Here, the two-electron wave functions are simply the symmetrized and anti-symmetrized products of the donor ground state orbitals. Early work by Koiller et al. 32 and Wellard et al. 33 approximated the donor potential by a simple Coulomb potential and ignored the valley-orbit coupling (VOC) terms. Extending on this work, Wellard et al. 34 and Pica et al. 35 included VOC and an isotropic central cell correction (CCC) to correct for the short range limit where the dielectric constant becomes meaningless and the impurity deforms the charge density that the electron interacts with. However, they failed to precisely reproduce the excited states within the 1s subspace. In the Hartree-Fock (HF) method the many-electron wave function is constructed from an antisymmetrized product of optimized molecular orbitals. Wu et al. 36 applied HF to an EMT implementation that includes a CCC and VOC. Their effective mass tensor is anisotropic, while their Gaussian basis is not.
To go beyond the HF, the Configuration Interaction (CI) method systematically adds more configurations, i.e. the symmetrized and antisymmetrized products of orbitals. The Hund-Mulliken approach adds ionic configurations to the two-electron basis. 37 Including all possible configurations is called Full CI and it is the optimal solution in a given valley-orbit basis. 18,30,31 Gonzalez-Zalba et al. 38 and Saraiva et al. 18 also implemented Full CI to model exchangecoupled donors. However, they used a hydrogenic model, i.e. an isotropic basis and effective mass, and neglected VOC. Additionally, they only applied a CCC to the donor ground state and assumed the excited states to be degenerate.
Our theory builds upon the multi-valley EMT implementation by Gamble et al.. 23,26 They solved the coupled six-valley Shindo-Nara equations, including silicon's full Bloch functions. 16 They introduced a new tetrahedral CCC and went beyond a minimal basis set. As a result, they were able to precisely reproduce the entire spectrum of valley-orbit split 1s states. The efficient use of Gaussian type orbitals allowed them to study ∼ 1.3 × 10 6 relative lattice placements of a single-electron two-phosphorus donor molecule.
Here, we add an efficient HF and Full CI solver to study exchange-coupled phosphorus donors. We introduce a efficient basis set of contracted Gaussian orbitals, to reduce computational costs. We perform a basis set analysis and find a fast and accurate choice. This allows us to perform exhaustive 3D iterations over all physically relevant relative donor placements for a given distance. We study the exchange interaction and valley configurations in detail for the [100], [110] and [111] crystal orientations and visualize the electron density. We find that a CROT gates are viable over a range of distances between 10 nm and 24 nm depending on the crystal direction. Along [100] and [110] the donor placement may straggle by 5 nm. Finally, we study the electric tunability of exchange-coupled donors. Using a simple noise model, we find that high-quality SWAP gates are viable.
The paper is organized as follows. In Sec. II, we review the multi-valley effective mass theory as used by Gamble et al.. 23,26 Sec. III introduces the HF and CI methods, applied to donors in silicon. In Sec. IV, we discuss our choice of basis functions. In Sec. V, we use our theory to calculate the key properties of exchange-coupled donors. Sec. VI compares our results to earlier work. Sec. VII discusses the exchange tunability with electric fields. In Sec. VIII, we summarize our main conclusions. This section gives a brief review of the underlying effective mass theory established by Gamble et al.. 23 The cornerstone of EMT is that low-energy conduction electron states ψ( r) in silicon only have support around the six valley minima (Fig. 1c): Here, we sum over the six valley minima at k µ with the corresponding Bloch functions ϕ µ ( r) = u kµ ( r)e i kµ r . The envelope functions F µ ( r) are subject to the multivalley EMT equations, a system of coupled Schrödinger-like equations: The kinetic energy operator of the µth valley T µ gives a different effective mass to motions along (m ) and perpendicular to (m ⊥ ) k µ , e.g.
is the potential of the impurity. The EMT equations are coupled via the valley-orbit potential which lifts the degeneracy of the valley states and results in the A 1 , T 2 and E valley-orbit states depicted in Fig. 1b. The relative amplitudes of F µ ( r) for each of these states are governed by the the symmetry of the crystal, namely T d . We give the relative amplitudes in vector form in natural ordering (+x, −x, +y, ...), Neglecting the electronic spin degeneracy, there is a non-degenerate ground state belonging to the A 1 irreducible representation of T d below a triplet of states belonging to T 2 and a doublet of states belonging to E. In practice, we need a suitable approximation for the product of Bloch functions ϕ * µ ( r)ϕ ν ( r) and to find the right binding potential U ( r) for our donor. The latter consists of the bulk-screened Coulomb potential U c ( r) = −e 2 /(4π 0 r r), where e is the electron charge, 0 is the permittivity of free space, r the dielectric constant of the host material, and r = | r| the distance from the nucleus of the substitutional donor.
In the short-range limit, i.e. near the unit cell containing the donor atom, the approximation of a bulk-screened Coulomb potential described by r breaks down. The impurity deforms the charge density the valence electrons it interacts with. To obtain the energies in Fig. 1b we need to include corrections close to the impurity.
Here, we use a central cell correction (CCC) in the style proposed by Gamble et al. 23 where we sum over all tetrahedral bond directions t i ∈ {(1, 1, 1), (−1, 1, −1), (1, −1, −1), (−1, −1, 1)}. All CCC parameters need to be found via optimization (see Sec. IV). Fig. 1d shows the constant energy surfaces of the CCC given in Tab. I. The Bloch functions can be expanded in a Fourier series, where the sum goes over reciprocal lattice vectors G. The Fourier coefficients {A µ G } were calculated by Gamble et al. 23 using density functional theory. Additionally, they established that the series of products of Bloch functions is well-converged at | G − G | ≤ 4.4 × 2π/0.543 nm and can be truncated at this point. To solve Eq. (2), we expand the envelope functions F µ ( r) over a finite set of orbitals F µ a ( r), = µ a C µ a χ µ a ( r).
We also define χ µ a as a valley orbit basis state. In this basis, the EMT Hamiltonian has the following form: The valley-diagonal blocks in the first line contain the kinetic energy and the impurity potential. They are coupled via the valley-orbit elements in the second row. Finally, we can express Eq. (2) as a generalized eigenvalue problem, where S µν a,b = χ µ a |χ ν b = δ µ,ν F µ a |F ν b is the block diagonal overlap matrix.

III. QUANTUM CHEMISTRY
Quantum chemistry is the field of science that seeks to describe the electronic structure of molecules from first principles using quantum mechanics. This section will present the Hamiltonian of a two-donor molecule. Then, we introduce computational chemistry methods, viz. Hartree-Fock (HF) and Configuration Interaction (CI), applied to donors in silicon.

A. The donor molecule
The Hamiltonian of the donor molecule is given by Here, h(i) the one-electron effective mass Hamiltonian of the ith electron in the field of both donors as presented in Sec. II. It is identical for both electrons and the 1 indicates the dependence on the space and spin coordinates of electron number 1. The electron-electron repulsion operator is given by where r i is the position operator of the ith electron and all other quantities are the same as for the donor Coulomb potential in Sec. II. Despite the superficial simplicity of the Hamiltonian 14, the multi-electron nature of the problem renders it very complex, often resisting semi-analytic solution and typically requiring sophisticated computational methods like the ones outlined below.

B. Hartree-Fock Approximation
In the Hartree-Fock (HF) model, we treat each electron as an independent particle and only consider interactions to the mean field of the other electrons.

Slater Determinants
We describe each electron by a spin orbital φ i , a one-electron wave function that is the product of a spatial orbital ψ i and the spin, which can either be up v ↑ or down v ↓ . In the simplest case we can have only one spin-up φ ↑ (1) = ψ ↑ (1)v ↑ (1) and one spin-down φ ↓ (1) = ψ ↓ (1)v ↓ (1) orbital. The many electron wave-function Φ(1, 2, ...) is a product of spin orbitals. It must not only satisfy the Schrödinger equation, but also be antisymmetric with respect to the interchange of the space and spin coordinates of any two electrons: This condition can conveniently be satisfied by arranging the spin orbitals in a Slater determinant (SD). In our two-electron example this is simply given by: Allowing each spin to have a different spatial orbital, as above, results in the Unrestricted Hartree-Fock (UHF) approximation. Restricting the spatial orbitals for each spin to be identical, i.e. ψ ↑ = ψ ↓ , is called Restricted Hartree-Fock (RHF).

Hartree-Fock Equations
To find the molecular orbitals (MOs) φ i that minimize the energy of a single Slater Determinant, we need to solve a set of effective one-electron equations, the HF Equations: This is the pseudo-eigenvalue problem, as the Fock operator itself depends on the MOs as outlined below. In this case it is given by Here, we introduced the Coulomb J i and exchange K i operators, which we define via the following integrals The integral over the Coulomb operator is just the classical repulsion between the charge distributions φ 2 j and φ 2 i . Note how the K operator exchanges the two orbitals on the right-hand side. The exchange integral has no classical analogy. Intuitively, the Fock operator can be interpreted as follows: The up spin electron sees the Coulomb potential from electrons in spin up and down orbitals plus the exchange potential from the spin up orbital. There is no exchange interaction between electrons of different spin. In RHF, the Fock operator is equal for both spin orbitals and Eq. (19) simplifies to We will focus on this case in the following.

Self-Consistent Field Iterations
The Fock operator itself depends on the orbitals φ i and Eq. (18) can only be solved iteratively.
In practice, the MOs are expanded over a finite orbital basis set, just like the EMT envelope functions in Eq. (11). The matrix representation of Eq. (18) in the valley-orbit basis yields the Roothaan-Hall equations the generalized eigenvalue problem, where S µ,ν a,b = χ µ a |χ ν b and F µ,ν a,b = χ µ a |F|χ ν b are µ-ν-valley matrix blocks of the Overlap and Fock matrix respectively. The first step to self-consistently solve Eq. (23) is to precompute the two-electron integrals. They are approximately valley local, 25,32,33 as the fast oscillatory part of the Bloch function, suppresses the contribution of mixed valley charge distributions, i.e. for k µ = k ν . We also introduced the four-index tensor (ab|cd) µκ , a convenient chemistry notation. Given an initial orbital ψ i (1) = µ a C µ a,i χ µ a (1), we find Note that the J matrix is valley diagonal. It can be interpreted as the mutual electrostatic repulsion between the classical charge densities in each valley.

C. Configuration Interaction
The HF method finds the energetically best solution a single SD can offer. However, this underestimates how the motion of the electrons influence each other, i.e. the electron correlation (EC). The Configuration Interaction (CI) method accounts for the EC energy by expanding the full multi-electron wave function Ψ over a set of SDs Each of these SDs represents an electron configuration. The full multi-electron Hamiltonian (Eq. (14)) is expressed in this basis. By diagonalization, we determine how these configurations interact. Including all possible configurations is called Full CI and it is the optimal solution in a given basis set.

Spin-Adapted Configurations
Our Hamiltonian in Eq. (14) does not contain any spin operators and S 2 and S z commute with H. As a result, the exact eigenfunctions of H are also eigenfunctions of the spin operators. In the two-electron case, these are the well-known singlet S (S 2 = 0) and triplet T (S 2 = 2) spin functions: In the previous section, we solved the RHF Roothaan's equations and found a set of MOs ψ i . The RHF ground state configuration was simply given by ψ ↑ = ψ ↓ = ψ 0 in Eq. (17), which corresponds to a spin singlet as expected. In principle, we can build a set of excited SDs Φ ij by promoting the up and down spin orbitals to excited MOs, ψ i and ψ j . However, configurations i = j are not eigenfunctions of the total spin operator S 2 . We use spin-adapted configurations (SACs) instead, which conveniently split H into singlet and triplet blocks. The spatial part of the singlet SACs are and for the triplet  where the latter requires i = j. Together with the spin functions in Eq. (29) and Eq. (30) they form the orthonormal basis in which we express our two-electron Hamiltonian: Note that using Full CI, it does not matter if we employ solutions of Eq. (13) or Eq. (23).

CI matrix elements
All we have to do do now is to express Eq. (14) in the singlet and triplet SAC basis. The first step is to express the one-and two-electron integrals in the MO basis Performing the transformation of the four-index tensor one index at a time reduces the complexity to O(N 5 basis ). With this, we can easily compute the the SAC matrix elements, which in the triplet case they are given by The singlet elements H S are a little more complex due to the cases i = j in Eq. (31). By solving the eigenvalue problem we find the singlet/triplet E S/T n energies as well as the corresponding coefficients a S/T ij (Eq. (33) /Eq. (34)) of the nth singlet/triplet state.
In the lower panel of Fig. 2 we perform the CI calculations including configurations in Eq. (31) up to and including the MO in the top panel. The resulting electron density is shown for z = 0. Including only a single MO is the RHF solution. Having only one orbital to work with, the RHF solution spreads it over both donors and introduces unnaturally high contributions from the E valley-orbit states, which are wide along the axis joining the donor. The resulting energy of E RHF = −98.50 meV is consequently quite high. Adding only the first excited MO shifts the electron density onto the donors, reducing the E-states contributions and lowering the energy to E CI = −107.03 meV. Adding more MOs increases this effect. The electron density of the Full CI solution, i.e. including all 36 MOs, looks much more like the expected solution of two A 1 -like states. 23 The resulting energy is lower at E CI = −109.64 meV. The energy difference between RHF and Full CI is called correlation energy and notably accounts for 10 % of the total energy within this basis set.

IV. BASIS SETS
To solve the EMT equations and to find our MOs, we expanded our envelope functions over a set of valley-orbit basis functions F µ a ( r) (Eq. (11)). In practice, it is impossible to use a complete basis, as we can only work with a finite number of functions. The size and the type of the basis we chose will determine the level of accuracy we can achieve.
Ideally, every single basis function already reproduces the shape of the wave function we would like to model. In our case, these are Slater type orbitals (STOs), which mimic the exact orbitals for the hydrogen atom. In the silicon crystal, the STOs are anisotropic; 41 a +z-valley orbital centered at the origin has the form: where N is a normalization factor, the n i are the exponents of polynomial prefactors and the α the anisotropic decay constants. Gaussian type orbitals (GTOs) allow for efficient evaluation of the one and two-electron integrals. The the one-and two-centered Coulomb kernels can be accurately decomposed into a compact Gaussian quadrature. In a GTO basis, the resulting matrix elements are then a sum of separable integrals over Gaussian functions, which can be computed analytically. 42 The GTO analogue to the STO above has the following shape The GTOs have a zero derivative at the center and are missing the cusp of a STO, making it hard to predict the proper behaviour close to the nucleus. Additionally, GTOs have a steeper flank. As a rule of thumb one needs three times as many GTOs as STOs to achieve the same accuracy. 31 This makes the transformation in Eq. (36), which scales as O(N 5 basis ), particularly costly in a GTO basis. The standard solution to this problem is to fit a number of n GTOs to a STO and contract them to represent one orbital. A single STO-nG basis function then has the following form where the N i and β i are the fit parameters (see Fig. 14). We choose n = 3 as adding more GTOs to the fit gives little improvement. The basis functions in the 5 other valleys are found via permutation of the exponents and polynomials. aving decided on the style of basis, we need to find a set of STO-nGs that models all features of the two-donor system with sufficient accuracy. We follow Gamble et al. 23 and variationally fix a CCC and basis that reproduces the energy levels of a single phosphorus donor. We try a small basis (SB) containing only 1s-orbitals (n x = n y = n z = 0). The CCC in Tab. I together with the orbitals #1 to #2 produce the D 0 energies on the right. The SB shows good agreement in the A 1 ground state. However, it was not possible to match the E and T 2 energies at the same time.
We proceed to test a larger basis (LB) that additionally contains 2s-orbitals (n x = n y = n z = 2). The CCC and orbitals #1 to #5 in Tab. II achieve good accuracy for all D 0 valley configurations.
To be able to model configurations with two electrons loaded on one donor, we optimize an additional orbital to produce the D − energy. We find that it is sufficient to add orbital #3 to the SB and #6 to the LB to achieve the D − energies on the right.  approach twice the energy of an isolated donor for large distances, while the singlet-triplet splitting decays. The second row shows the difference ∆E between the singlet/triplet energies calculated using the SB and the LB, while the third row compares the exchange interaction strength where F µ a is the ath basis state in the µth valley, which differs from the ath basis state in another valley through the permutation of the decay constants and polynomial prefactors. The SB 1s-orbitals have faster decaying flanks, i.e. the LB orbitals are flatter. As a result, the LB states put more support on the faster decaying red orbital to offset the steeper flank of the primary (blue lines in Fig. 3) orbital. Note that although the 2s-orbital contributions are small, they have high amplitudes at the CCC and are crucial to get the right valley energies. For d 5 nm the donor molecule is not dissociated yet. The narrow (red lines in Fig. 3) orbitals dominate the wave function. This can be understood considering the hydrogen atom. The ground state electronic wavefunction has a radius of r = a B = 0.53 Å. When we bring two hydrogen atoms together to zero distance, we get the limiting case of a helium atom, which has a narrower radius of r = 0.31 Å. Equally the short distance limit of a donor molecule is narrower, leading to the greater contributions of red orbitals. Both electrons are spread over the increasingly helium-like donor cluster. Neither our CCC nor our basis is designed for this limit, resulting in larger ∆E. For larger distances, the system acts more like two separate donors with one A 1 electron on each (see Sec. V). Their energies have been successfully optimized for both basis sets, resulting in ∆E nicely converging to zero.  The LB exchange interaction strength decays slightly faster due to the flatter nature of the relevant basis states (red and blue) that makes the difference between the bonding and antibonding states less pronounced. We conclude that for the study of two exchange-coupled donors with separation d > 5 nm both basis sets are equally viable. Fig. 4b shows the difference ∆J as a function of J. The diagonal dashed lines mark the ∆J/J = 100% and 10% regimes. We observe that it is very difficult to precisely compute very small values of J. This is a natural effect and will happen in all types of calculations.
Finally, let us consider the computational cost. The simulations were performed on the Gadi supercomputer. 43 A LB iteration typically costs 120 SU ('Service Units' or computer core hours), while a SB step costs 4.7 SU. We here showed that the SB is sufficient for d 5 nm and we can use it to perform exhaustive 3D iterations (see Fig. 6) including N ≈ 8500 iterations run with approximately 40 kSU. The carbon footprint of all simulations run in this project is calculated in B.

V. STUDY OF EXCHANGE-COUPLED DONORS
We now have an appropriate CCC and basis set and are ready use our Full CI code to study two exchange-coupled donors. Fig. 5a shows the exchange interaction strength as a function of donor separation along different crystal axes sketched in Fig. 5b.
Notably, in the large-distance regime we find that J [100] (d) > J [110] (d) > J [111] (d). This is due to the anisotropy of the envelope functions F µ ( r) (Eq. (1)) 41 as outlined in the following. Each F µ ( r) has an ellipsoidal shape squeezed along the valley orientation. The +x-valley combination of the A 1 ground state of the Bloch and envelope function |F +x ( r)ϕ +x ( r)| 2 on the surface of a sphere with r = 10 nm around the donor is shown in Fig. 5c. Fig. 5d shows the full A 1 ground state. Each valley contributes a belt perpendicular to its axis. Coming from large distances, we have two A 1 states approaching each other and the exchange interaction strength will scale as the overlap of their wave functions. Along the [100] crystal axis (x-axis) two of these belts overlap and result in a high J [100] . Along   [110] (bisecting x and y-axis) only the z-valley ellipsoids contribute fully giving an intermediate J [110] , while going diagonally along [111] corresponds to the dark region in Fig. 5d and a small J [111] .
Another interesting feature is that the exchange interaction strength has a dent at roughly d ≈ 4 nm, especially for [100]. At this point the two lowest spin triplet states have an avoided crossing. 22 This occurs when the spin triplet arising from an antibonding combination of A 1 orbitals (whose energy increases with decreasing distance) crosses the energy of the bonding combination of T 2 orbitals (whose energy decreases with decreasing distance). The consequence of this hierarchy inversion has been observed in an experiment with closely-spaced implanted donors. 44 The first row of Fig. 5e shows the three lowest energy states E S 0 (d) and E T 0/1 (d) as a function of donor separation along the three crystal directions. The splitting between the A 1 -like singlet (blue) and triplet (yellow) increases until the triplet meets the A CROT operation relies on the magnetic drive of a single electron spin resonance line. On the one hand, J needs to be large enough to allow the individual drive on a specific spectral line giving us a lower bound to J. A large J will however induce deviations from the logical qubit basis given by the simple spin product states providing an upper boud to J. A high-fidelity CROT operation requires 0.1 MHz < J < 10 MHz, Eq. (6). An A 1 -like state of a donor at position r b has components in analogy to Eq. (11) , where a iterates the basis states in the µth valley. N is a normalization constant, A 1 (µ) gives the relative valley amplitudes in Eq. (4) and r a is the center of the ath basis state. We compute the A 1 -valley contribution where b sums over all donor positions, and similarly T 2 and E. The second/third row of Fig. 5e gives valley contributions of the lowest singlet/triplet. For large distances the electronic ground states of two separate donors are just A 1 like. Bringing them together breaks the A 1 symmetry. Doing so along [100], e.g. the x-axis, will energetically favour y-and z-valley orbitals that are wide along this axis, and promote E symmetry (Eq. (6) E contributions, while the directional symmetry of the x-and y-valleys is broken giving small T 2 contributions to the ground state. Again the triplet T 2 contributions oscillate for short distances. . This leads to less support in between the donors reducing the singlet-triplet difference, i.e. the exchange interaction strength as discussed above. Additionally, at d ≈ 4 nm, the triplet state becomes more T 2 -like which changes the lattice periodic structure of ρ T . The bonding/anti-bonding structure of the ρ S/T ceases to exist leading to the dent in the exchange interaction, discussed above. Along [111], the second donor interchanges between inverted tetrahedral bonds, i.e. alternates between the shifted fcc sublattices. The electron density around this donor appears to periodically flip.
Finally, we perform exhaustive 3D iterations. In Fig. 9 we compute exchange interaction strength for one donor at the center of the shown spheres with r = 5 − 15 nm and the second on the lattice site closest to each point on the surface. We exploit the symmetry of the silicon crystal to reduce the number of unique calculations that are performed by a factor of 48. We can nicely see all the features discussed above.  38 Early work by Koiller et al. 32 and Wellard et al. 33 in Fig. 10a and b used a Heitler-London approach, i.e. a single SD formed by the one-electron wave functions centered on each donor to model the singlet and triplet states. They approximated the impurity by a simple bulk-screened Coulomb potential and ignored VOC terms. They use a minimal basis with one anisotropic STO per valley per donor. The resulting exchange in Fig. 10a is roughly an order of magnitude higher than the present work. Wellard et al. 33  Still using Heitler-London, Wellard et al. 34 and Pica et al. 35 included isotropic CCCs from Ref. 15 and 14 respectively. They included VOC terms, but still used the same minimal basis set. Their EMT failed to precisely reproduce the excited valley state energies. Again, the results by Wellard et al. 34 in Fig. 10c are consistent with ours, while the exchange in Fig. 10d is over an order of magnitude smaller. J (GHz) Gonzalez-Zalba et al. 38 and Saraiva et al. 18 implemented a Full CI as the current work. However, they used a hydrogenic model, i.e. an isotropic basis and effective mass. As a result, the exchange in Fig. 10e has similar amplitudes along all directions, i.e. J [100]

VI. COMPARISON TO LITERATURE
Their basis was made from one STO-3G orbital per valley per donor and they neglected VOC. Additonally, they applied a CCC that interpolates between the bulkscreened and vacuum case to the donor ground state and assumed the excited states to be degenerate. In this model, the exchange decays quickly for all directions comparable to our [111] results. Notably, the exchange along [100] is not flat there, but exhibits a step like decay.
Wu et al. 36 used an EMT including VOC and a CCC in the style of Saraiva et al.. 18 While their basis consisted of five isotropic GTOs per valley per donor, their kinetic energy operator included the effective mass anisotropy. Again, these approximations lacked the precision to produce the excited valley state energies. To obtain the results in Fig. 10f, they used the HF method described in Sec. III B. Their [110] results are comparable to our [111] calculations. Note that the exchange seems to saturate at large distances.

VII. ELECTRIC TUNABILITY AND SWAP GATES
In the previous section we discussed a pair of exchange-coupled phosphorus donors as a function of their distance d, in the absence of electric fields. Now we investigate the electric field dependence of J, with the goal of understanding under which conditions it becomes possible tune the exchange interaction to perform two-qubit SWAP gates. [6][7][8] We add an electric field potential to our EMT Hamiltonian in Eq. (12) where E is the electric field vector. We neglect VOC terms for U E , due to their fast oscillatory behaviour. In the following, we will apply a field along the axis joining the donors and define the detuning = | E|d, where d is the donor separation. We tested the effect of polarization when p-orbitals (n i = 1) are included. We found that polarization is only a minor effect and to first order the system can be explained in a molecular orbital approximation containing covalent and ionic contributions, i.e. D 0 and D − orbitals. Fig. 11a and e show exemples of singlet/triplet energies of two donors positioned along [100] at d = 20.09 nm as a function of detuning. States with one electron on each donor are sensitive to a detuning and have a slope, while ionized states are approximately flat. With increasing detuning, the lowest singlet and triplet states (blue) increase in energy until the ionized state becomes energetically favourable. The singlet allows closed-shell configurations, i.e. i = j in Eq. (31). As a result, the ionized singlet is just A 1 -A 1 -like. The triplets in Fig. 11e do not allow such solutions and the lowest ionized state is A 1 -T 2 -like and we are missing the equivalent to the A 1 -A 1 -like singlet shown in grey. As a result, we need to invest more energy to ionize the triplet and the ionization point in Fig. 9f is shifted compared to the singlet in Fig. 11b. The valley contributions (Eq. (46)) in Fig. 9c and g tell the same story. For zero detuning,  the singlet and triplet are both A 1 like. At the ionization point, the triplet becomes A 1 -T 2 -like. Fig. 11d and h give the state contributions (Eq. (44)) as a function of detuning. At the ionization point the states occupy D − -like basis functions. We now analyze the viability of tunable exchange-coupled donors to perform SWAP gates. The SWAP oscillations of two spins prepared in |↑↓ are described by the Rabi formula 45 where J is the exchange interaction strength and ∆B z the longitudinal magnetic field gradient across the two electrons.
Tuning J relative to ∆B z allows us to tune the prefactor in Eq. (48) and switch the SWAP oscillations on or off. In a donor system, the electron and the nuclear spin are coupled by a hyperfine interaction A ≈ 120 MHz, which in a gated nanoscale device can vary from one donor to the next by an order ∆A ≈ 1 MHz. This allows the introduction of an intrinsic ∆B z , switchable by controlling the state of the two nuclear spins. 9,11 If the nuclei are parallel, a small ∆B z = ∆A ≈ 1 MHz is created by the difference in the hyperfine interaction between the two donors. A much larger ∆B z =Ā ≈ 120 MHz can be obtained by preparing the nuclei in an antiparallel state. Fig. 12a shows the exchange interaction strength (Eq. (43)) as a function of detuning for donors positioned along [100] at distances between 10 nm and 25.5 nm. Here, the smallest distance belongs to the curve with the highest J. In the unperturbed system, the exchange interaction strength is the weakest for a given distance and increases as a function of . After the singlet ionization point (see Fig. 9), the exchange interaction has a steep increase and finally saturates after the triplet ionization. Fig. 12b and d plots the prefactor in Eq. (48) F = J 2 /(J 2 + ∆B 2 z ) for ∆B z = 120 MHz and ∆B z = 1 MHz. For a high fidelity SWAP gate we demand to tune this from F = 0.01 to 0.9 indicated by the dashed lines. To perform the gate, we prepare the system anywhere below the F = 0.01 line, tune it to the F = 0.9 line, then let it evolve for T π = π/ J 2 + ∆B 2 z and finally pulse it back below the F = 0.01 line. The corresponding exchange interaction tuning ranges are indicated in Fig. 12a. To estimate the quality of the SWAP oscillations at F = 0.9, we assume quasistatic noise, where the decay time is given by 46 where we take rms = 24 × 10 −3 GHz nm −1 × d as our noise amplitude. [47][48][49] We then define the quality factor as the number of SWAP oscillations we can perform in T * Fig. 12c, e give the quality factor Q at the F = 0.9 point in Fig. 12b, d for viable distances. The quality factors are typically larger than 100 and on the same order of magnitude for parallel and antiparallel nuclear spin configurations. Picking distances larger than the minimal requirement, will result in a steeper slope of J at the driving point in Fig. 12a. This increases the sensitivity to charge noise, which leads to a shorter T * 2 and the decrease in quality factor with distance in Fig. 12c,e. Fig. 13 produces the same graphs for two tunable donors positioned along [110]. Here, the same range of exchange coupling is obtained at shorter distances due to the effective mass anisotropy (see Sec. V). Additionally, we can see the valley oscillations, which cause the irregular decrease of J with distance. We obtain quality factors of the same order of magnitude as for the [100] crystal orientation. In this configuration parallel nuclear spins yield slightly better results. However, the using such small values of exchange interaction would require slow 1 µs pulses to avoid spectral broadening beyond the value of J, whereas antiparallel spins allow much faster gate operations, on the order of 100 ns.

VIII. CONCLUSION
We presented Full Configuration Interaction simulations of exchange-coupled donors using a state-of-the-art implementation of multi-valley effective mass theory. We tested two basis sets and found the set summarized in Tab. I to be sufficient for donor separations d 5 nm. This allowed us to perform 3D iterations including 8500 donor configurations. We studied in detail two donors positioned along [100], [110] and [111]. We analysed the exchange interaction, the valley configurations and visualized the evolution of electron density as a function of distance. We found that a high-fidelity CROT gates are possible over a range of distances between 10 nm and 24 nm depending on the crystal direction. Along [100] and [110] the donor placement may straggle by 5 nm. Finally, we explored the electric tunability of exchange-coupled donors. High-fidelity SWAP gates require the system to be pulsed close to the ionization point. However, using a simple noise model we predict that high-quality SWAPs are possible. This work shows how efficient modelling can inform the design of two-qubit devices.