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

Benjamin Joecker  Centre for Quantum Computation and Communication Technology,
School of Electrical Engineering & Telecommunications, UNSW, Sydney, NSW,
2052, Australia  Andrew D. Baczewski  Center for Computing Research, Sandia
National Laboratories, Albuquerque, NM, 87123, USA  John K. Gamble  Center for
Computing Research, Sandia National Laboratories, Albuquerque, NM, 87123, USA
Microsoft Quantum, 1 Microsoft Way, Redmond, WA, 98052, USA  Jarryd J. Pla
Centre for Quantum Computation and Communication Technology, School of
Electrical Engineering & Telecommunications, UNSW, Sydney, NSW, 2052,
Australia  André Saraiva  Centre for Quantum Computation and Communication
Technology, School of Electrical Engineering & Telecommunications, UNSW,
Sydney, NSW, 2052, Australia  Andrea Morello  Centre for Quantum Computation
and Communication Technology, School of Electrical Engineering &
Telecommunications, UNSW, Sydney, NSW, 2052, Australia

###### Abstract

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 28Si, it has demonstrated exceptional coherence
timesMuhonen  _et al._ (2014). High fidelity initialization, readout Morello
_et al._ (2010), and single-qubit operationsPla  _et al._ (2012); Muhonen  _et
al._ (2015); Dehollain  _et al._ (2016) are well established. Several
proposals for multi-qubit gates rely on the exchange interaction, using either
strong tunable exchange to perform SWAP oscillations Kane (1998); Hollenberg
_et al._ (2006); Hill (2007), or weak exchange and microwave pulses to perform
a CROT gate Hill (2007); Kalra  _et al._ (2014). Experimental work by He et
al.He  _et al._ (2019) demonstrated fast SWAP oscillations in the strong
exchange regime, using two donor clusters fabricated with a scanning tunnel
microscope. Madzik et al.Mądzik  _et al._ (2020) 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 31P 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.

Effective mass theory (EMT) was originally developed to find approximate
solutions for the complex wave-function of a donor-bound electron in
siliconKohn and Luttinger (1955); Luttinger and Kohn (1955) and has been
continuously developedNing and Sah (1971); Pantelides and Sah (1974); Shindo
and Nara (1976); Hui (2013); Saraiva  _et al._ (2015). It has been applied to
donors near an interfaceCalderón  _et al._ (2009); Baena  _et al._ (2012) and
single-electron donor moleculesHu  _et al._ (2005); Klymenko and Remacle
(2014); Gamble  _et al._ (2015); Klymenko  _et al._ (2017); Saraiva  _et al._
(2015). It is widely used to model physical properties, such as valley
splittingSaraiva  _et al._ (2011); Gamble  _et al._ (2016), tunnel couplingHu
_et al._ (2005); Gamble  _et al._ (2015), and both StarkFriesen (2005);
Debernardi  _et al._ (2006) and hyperfinePica  _et al._ (2014a) shifts.

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 Szabo and Ostlund (2012); Jensen (2017) to
the multi-valley effective mass theory frame work. Widely used is the Heitler-
London approachKoiller  _et al._ (2001); Wellard  _et al._ (2003); Wellard and
Hollenberg (2005); Pica  _et al._ (2014b). 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.Koiller  _et al._
(2001) and Wellard et al.Wellard  _et al._ (2003) approximated the donor
potential by a simple Coulomb potential and ignored the valley-orbit coupling
(VOC) terms. Extending on this work, Wellard et al.Wellard and Hollenberg
(2005) and Pica et al.Pica  _et al._ (2014b) 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.Wu and
Fisher (2020) 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 basisKettle  _et al._ (2006). Including all possible configurations
is called Full CI and it is the optimal solution in a given valley-orbit
basisSzabo and Ostlund (2012); Jensen (2017); Saraiva  _et al._ (2015).
Gonzalez-Zalba et al.González-Zalba  _et al._ (2014) and Saraiva et al.Saraiva
_et al._ (2015) also implemented Full CI to model exchange-coupled 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.Gamble  _et al._ (2015, 2016). They solved the coupled six-valley Shindo-
Nara equations, including silicon’s full Bloch functionsShindo and Nara
(1976). 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 $\sim$1.3\text{\times}{10}^{6}\text{\,}$$ 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\text{\,}\mathrm{nm}$ and $24\text{\,}\mathrm{nm}$
depending on the crystal direction. Along [100] and [110] the donor placement
may straggle by $5\text{\,}\mathrm{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.Gamble  _et al._ (2015, 2016).
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.

##  II Effective Mass Theory

Figure 1: (a) Illustration of an exchange-coupled pair of phosphorus donors
(red) in a lattice of silicon atoms (silver) (b) Experimentally measured
energy levels of a negativeNarita  _et al._ (1982) (red) and neutral
donorGrimmeiss  _et al._ (1982) (blue). The degeneracy of the six conduction
band minima (‘valleys’) illustrated in panel (c) is lifted by the valley-orbit
coupling. (d) Constant energy surfaces $E=$-1\text{\,}\mathrm{eV}$$ of the
impurity potential with the tetrahedral central cell correction in Eq. (7)
with parameters given in Tab. 1.

This section gives a brief review of the underlying effective mass theory
established by Gamble et al.Gamble  _et al._ (2015). The cornerstone of EMT is
that low-energy conduction electron states $\psi(\vec{r})$ in silicon only
have support around the six valley minima (Fig. 1c):

$\displaystyle\psi(\vec{r})$
$\displaystyle=\sum_{\mu}F^{\mu}(\vec{r})\varphi_{\mu}(\vec{r}).$ (1)

Here, we sum over the six valley minima at $\vec{k}_{\mu}$ with the
corresponding Bloch functions
$\varphi_{\mu}(\vec{r})=u_{\vec{k}_{\mu}}(\vec{r})e^{i\vec{k}_{\mu}\vec{r}}$.
The envelope functions $F^{\mu}(\vec{r})$ are subject to the multivalley EMT
equations, a system of coupled Schrödinger-like equations:

$\displaystyle[\mathbf{T}_{\mu}+U(\vec{r})]F^{\mu}(\vec{r})+\sum_{\mu}V^{VO}_{\mu\nu}(\vec{r})F^{\nu}(\vec{r})=EF^{\mu}(\vec{r}).$
(2)

The kinetic energy operator of the $\mu$th valley $\mathbf{T}_{\mu}$ gives a
different effective mass to motions along ($m_{\parallel}$) and perpendicular
to ($m_{\perp}$) $\vec{k}_{\mu}$, e.g.
$\mathbf{T}_{+x}=-\frac{\hbar^{2}}{2m_{\parallel}}\frac{\partial^{2}}{\partial
x^{2}}-\frac{\hbar^{2}}{2m_{\perp}}(\frac{\partial^{2}}{\partial
y^{2}}+\frac{\partial^{2}}{\partial z^{2}})$. $U(\vec{r})$ is the potential of
the impurity. The EMT equations are coupled via the valley-orbit potential

$\displaystyle
V^{VO}_{\mu\nu}(\vec{r})=\varphi^{\ast}_{\mu}(\vec{r})\varphi_{\nu}(\vec{r})U(\vec{r}),$
(3)

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^{\mu}(\vec{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$, …),

$\displaystyle A_{1}\quad\frac{1}{\sqrt{6}}(1,1,1,1,1,1),$ (4) $\displaystyle
T_{2}\quad\begin{cases}\frac{1}{\sqrt{2}}(1,-1,0,0,0,0)\\\
\frac{1}{\sqrt{2}}(0,0,1,-1,0,0)\\\
\frac{1}{\sqrt{2}}(0,0,0,0,1,-1),\leavevmode\nobreak\ \text{and}\end{cases}$
(5) $\displaystyle E\;\quad\begin{cases}\frac{1}{2}(1,1,-1,-1,0,0)\\\
\frac{1}{2}(1,1,0,0,-1,-1).\end{cases}$ (6)

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 $\varphi^{\ast}_{\mu}(\vec{r})\varphi_{\nu}(\vec{r})$ and to find
the right binding potential $U(\vec{r})$ for our donor. The latter consists of
the bulk-screened Coulomb potential
$U_{c}(\vec{r})=-e^{2}/(4\pi\epsilon_{0}\epsilon_{r}r)$, where $e$ is the
electron charge, $\epsilon_{0}$ is the permittivity of free space,
$\epsilon_{r}$ the dielectric constant of the host material, and $r=|\vec{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
$\epsilon_{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.Gamble  _et al._ (2015)

$\displaystyle
U_{cc}(\vec{r})=A_{0}e^{-r^{2}/(2a^{2})}+A_{1}\sum_{i}e^{-|\vec{r}-b\vec{t}_{i}|^{2}/(2c^{2})},$
(7)

where we sum over all tetrahedral bond directions
$\vec{t}_{i}\in\\{(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. 1.

The Bloch functions can be expanded in a Fourier series,

$\displaystyle\varphi_{\mu}(\vec{r})=u_{\vec{k}_{\mu}}(\vec{r})e^{i\vec{k}_{\mu}\vec{r}}=e^{i\vec{k}_{\mu}\vec{r}}\sum_{\vec{G}}A_{\vec{G}}^{\mu}e^{i\vec{G}\vec{r}},$
(8)

where the sum goes over reciprocal lattice vectors $\vec{G}$. The Fourier
coefficients $\\{A^{\mu}_{\vec{G}}\\}$ were calculated by Gamble et al.Gamble
_et al._ (2015) using density functional theory. Additionally, they
established that the series of products of Bloch functions

$\displaystyle{\varphi_{\mu}}^{\ast}(\vec{r})\varphi_{\nu}(\vec{r})=\sum_{\vec{G},\vec{G}^{\prime}}\left(A^{\mu}_{\vec{G}^{\prime}}\right)^{\ast}A^{\nu}_{\vec{G}}e^{i(\vec{k}_{\nu}-\vec{k}_{\mu}+\vec{G}-\vec{G}^{\prime})\vec{r}}$
(9)

is well-converged at $|\vec{G}-\vec{G}^{\prime}|\leq 4.4\times
2\pi/$0.543\text{\,}\mathrm{nm}$$ and can be truncated at this point.

To solve Eq. (2), we expand the envelope functions $F^{\mu}(\vec{r})$ over a
finite set of orbitals $F^{\mu}_{a}(\vec{r})$,

$\displaystyle\psi(\vec{r})$
$\displaystyle=\sum_{\mu}\sum_{a}C^{\mu}_{a}F^{\mu}_{a}(\vec{r})\varphi_{\mu}(\vec{r})$
(10) $\displaystyle=\sum_{\mu}\sum_{a}C^{\mu}_{a}\chi^{\mu}_{a}(\vec{r}).$
(11)

We also define $\chi^{\mu}_{a}$ as a valley orbit basis state. In this basis,
the EMT Hamiltonian has the following form:

$\displaystyle\begin{split}h^{\mu\nu}_{a,b}=&\delta_{\mu,\nu}\braket{F^{\mu}_{a}(\vec{r})}{\mathbf{T}_{\mu}+U(\vec{r})}{F^{\nu}_{b}(\vec{r})}\\\
&+\braket{F^{\mu}_{a}(\vec{r})}{V^{VO}_{\mu\nu}(\vec{r})}{F^{\nu}_{b}(\vec{r})}\end{split}$
(12)

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,

$\displaystyle\mathbf{h}\vec{C}_{i}=\epsilon_{i}\mathbf{S}\vec{C}_{i},$ (13)

where
$S^{\mu\nu}_{a,b}=\braket{\chi^{\mu}_{a}}{\chi^{\nu}_{b}}=\delta_{\mu,\nu}\braket{F^{\mu}_{a}}{F^{\nu}_{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.

###  III.1 The donor molecule

The Hamiltonian of the donor molecule is given by

$\displaystyle\mathbf{H}=\mathbf{h}(1)+\mathbf{h}(2)+\mathbf{r}_{12}(1,2)$
(14)

Here, $\mathbf{h}(i)$ the one-electron effective mass Hamiltonian of the $i$th
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

$\displaystyle\mathbf{r}_{12}(1,2)=\frac{e^{2}}{4\pi\epsilon_{r}\epsilon_{0}|\vec{r}_{1}-\vec{r}_{2}|}$
(15)

where $\vec{r}_{i}$ is the position operator of the $i$th 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.

###  III.2 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.

####  III.2.1 Slater Determinants

We describe each electron by a spin orbital $\phi_{i}$, a one-electron wave
function that is the product of a spatial orbital $\psi_{i}$ and the spin,
which can either be up $v_{\uparrow}$ or down $v_{\downarrow}$. In the
simplest case we can have only one spin-up
$\phi_{\uparrow}(1)=\psi_{\uparrow}(1)v_{\uparrow}(1)$ and one spin-down
$\phi_{\downarrow}(1)=\psi_{\downarrow}(1)v_{\downarrow}(1)$ orbital. The many
electron wave-function $\Phi(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:

$\displaystyle\Phi(1,..,i,...,j,...,N)=-\Phi(1,..,j,...,i,...,N)$ (16)

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:

$\displaystyle\Phi(1,2)$
$\displaystyle=\frac{1}{\sqrt{2}}\begin{vmatrix}\phi_{\uparrow}(1)&\phi_{\downarrow}(1)\\\
\phi_{\uparrow}(2)&\phi_{\downarrow}(2)\end{vmatrix}$ (17)
$\displaystyle=\frac{1}{\sqrt{2}}\left(\psi_{\uparrow}(1)v_{\uparrow}(1)\psi_{\downarrow}(2)v_{\downarrow}(2)-\psi_{\downarrow}(1)v_{\downarrow}(1)\psi_{\uparrow}(2)v_{\uparrow}(2)\right)$

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.
$\psi_{\uparrow}=\psi_{\downarrow}$, is called Restricted Hartree-Fock (RHF).

####  III.2.2 Hartree-Fock Equations

To find the molecular orbitals (MOs) $\phi_{i}$ that minimize the energy of a
single Slater Determinant, we need to solve a set of effective one-electron
equations, the HF Equations:

$\displaystyle\mathbf{F}_{i}\phi_{i}=\epsilon_{i}\phi_{i}.$ (18)

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

$\displaystyle\mathbf{F}_{i}=\mathbf{h}+\mathbf{J}_{\uparrow}+\mathbf{J}_{\downarrow}-\mathbf{K}_{i}.$
(19)

Here, we introduced the Coulomb $\mathbf{J}_{i}$ and exchange $\mathbf{K}_{i}$
operators, which we define via the following integrals

$\displaystyle\braket{\phi_{j}}{\mathbf{J}_{i}}{\phi_{j}}$
$\displaystyle=\braket{\phi_{i}(1)\phi_{j}(2)}{\mathbf{r}_{12}}{\phi_{i}(1)\phi_{j}(2)}$
(20) $\displaystyle\braket{\phi_{j}}{\mathbf{K}_{i}}{\phi_{j}}$
$\displaystyle=\braket{\phi_{i}(1)\phi_{j}(2)}{\mathbf{r}_{12}}{\phi_{j}(1)\phi_{i}(2)}.$
(21)

The integral over the Coulomb operator is just the classical repulsion between
the charge distributions $\phi_{j}^{2}$ and $\phi_{i}^{2}$. Note how the
$\mathbf{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

$\displaystyle\mathbf{F}=\mathbf{h}+2\mathbf{J}-\mathbf{K}.$ (22)

We will focus on this case in the following.

####  III.2.3 Self-Consistent Field Iterations

The Fock operator itself depends on the orbitals $\phi_{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

$\displaystyle\mathbf{F}\vec{C}_{i}=\epsilon_{i}\mathbf{S}\vec{C}_{i},$ (23)

the generalized eigenvalue problem, where
$S^{\mu,\nu}_{a,b}=\braket{\chi^{\mu}_{a}}{\chi^{\nu}_{b}}$ and
$F^{\mu,\nu}_{a,b}=\braket{\chi^{\mu}_{a}}{\mathbf{F}}{\chi^{\nu}_{b}}$ are
$\mu$-$\nu$-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 localKoiller  _et al._
(2001); Wellard  _et al._ (2003); Saraiva  _et al._ (2011),

$\displaystyle\begin{split}\braket{\chi^{\mu}_{a}\chi^{\kappa}_{c}}{\mathbf{r}_{12}}{\chi^{\nu}_{b}\chi^{\lambda}_{d}}&=\delta_{\mu\nu}\delta_{\kappa\lambda}\braket{F^{\mu}_{a}F^{\kappa}_{c}}{\mathbf{r}_{12}}{F^{\mu}_{b}F^{\kappa}_{d}}\\\
&=\delta_{\mu\nu}\delta_{\kappa\lambda}(ab|cd)_{\mu\kappa},\end{split}$ (24)

as the fast oscillatory part of the Bloch function,

$\displaystyle\varphi^{\ast}_{\mu}(1)\varphi_{\nu}(1)=u_{-\vec{k}_{\mu}}(\vec{r}_{1})u_{\vec{k}_{\nu}}(\vec{r}_{1})e^{i(\vec{k}_{\nu}-\vec{k}_{\mu})\vec{r}_{1}},$
(25)

suppresses the contribution of mixed valley charge distributions, i.e. for
$\vec{k}_{\mu}\neq\vec{k}_{\nu}$. We also introduced the four-index tensor
$(ab|cd)_{\mu\kappa}$, a convenient chemistry notation. Given an initial
orbital $\psi_{i}(1)=\sum_{\mu}\sum_{a}C^{\mu}_{a,i}\chi^{\mu}_{a}(1)$, we
find

$\displaystyle\braket{\chi^{\mu}_{a}}{\mathbf{J}}{\chi^{\nu}_{b}}=\delta_{\mu\nu}\sum_{\kappa}\sum_{c,d}{C^{\kappa}_{d,i}}^{\ast}C^{\kappa}_{d,i}(ab|cd)_{\mu\nu},$
(26)
$\displaystyle\braket{\chi^{\mu}_{a}}{\mathbf{K}}{\chi^{\nu}_{b}}=\sum_{c,d}{C^{\mu}_{c,i}}^{\ast}C^{\nu}_{d,i}(ad|cb)_{\mu\nu}.$
(27)

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.

The self-consistent field iterations now go as follows: (i) Obtain an inital
guess for the state $\psi_{0}$. (ii) Form the Fock matrix by application of
Eq. (26) and Eq. (27). (iii) Diagonalize the Fock matrix, i.e. solve Eq. (23).
(iv) Find the new ground state coefficients $C^{\mu}_{a,0}$ and iterate from
step (ii) until the cycle converges.

The upper panel of Fig. 2 shows the fully converged $|\psi_{i}(x,y,z=0)|^{2}$
for two phosphorus donors spaced by $d=$7.059\text{\,}\mathrm{nm}$$ along
[100]. The basis used is defined in Tab. 1 and contains 3 orbitals per valley
per donor resulting in a 36-dimensional problem.

###  III.3 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 $\Psi$ over a set of SDs

$\displaystyle\Psi=\sum_{i}a_{i}\Phi_{i}.$ (28)

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.

Figure 2: Illustration of the Configuration Interaction (CI) method for the
exemplary case of the singlet ground state of two phosphorus donors spaced by
$d=$7.059\text{\,}\mathrm{nm}$$ along [100]. The upper panel ((b), (d) and (f)
offline) shows the converged molecular orbitals (MO) $|\psi_{i}(x,y,z=0)|^{2}$
for the basis set in Tab. 1, where the slider iterates through $i$. In the
lower panel, the we perform the CI calculation including configurations in Eq.
(31) up to and including the MO in the top panel. We show the resulting
electron density
$\rho(\vec{r}_{1})=\int|\Psi^{S}(\vec{r}_{1},\vec{r}_{2})|^{2}d^{3}\vec{r}_{2}$.
Only including the first MO ((a) offline) corresponds to the Restricted
Hartree Fock (RHF) reference state, while including all MOs is the Full CI
result ((e) offline). Adding more MOs allows the CI to account better for the
electron correlation and lowers the CI energy. The radio button optionally
shows the shift of the CI electron density compared to the RHF reference
state. (Offline: Panel (c) gives an $z=y=0$ linecut as a function of included
number of MOs on the y-axis)

####  III.3.1 Spin-Adapted Configurations

Our Hamiltonian in Eq. (14) does not contain any spin operators and
$\mathbf{S}^{2}$ and $\mathbf{S}_{z}$ commute with $\mathbf{H}$. As a result,
the exact eigenfunctions of $\mathbf{H}$ are also eigenfunctions of the spin
operators. In the two-electron case, these are the well-known singlet $S$
($\mathbf{S}^{2}=0$) and triplet $T$ ($\mathbf{S}^{2}=2$) spin functions:

$\displaystyle S(1,2)$
$\displaystyle=\frac{1}{\sqrt{2}}(v_{\uparrow}(1)v_{\downarrow}(2)-v_{\downarrow}(1)v_{\uparrow}(2))$
(29) $\displaystyle T(1,2)$
$\displaystyle=\begin{cases}v_{\uparrow}(1)v_{\uparrow}(2)\\\
\frac{1}{\sqrt{2}}(v_{\uparrow}(1)v_{\downarrow}(2)+v_{\downarrow}(1)v_{\uparrow}(2))\\\
v_{\downarrow}(1)v_{\downarrow}(2)\end{cases}$ (30)

In the previous section, we solved the RHF Roothaan’s equations and found a
set of MOs $\psi_{i}$. The RHF ground state configuration was simply given by
$\psi_{\uparrow}=\psi_{\downarrow}=\psi_{0}$ in Eq. (17), which corresponds to
a spin singlet as expected. In principle, we can build a set of excited SDs
$\Phi_{ij}$ by promoting the up and down spin orbitals to excited MOs,
$\psi_{i}$ and $\psi_{j}$. However, configurations $i\neq j$ are not
eigenfunctions of the total spin operator $\mathbf{S}^{2}$. We use spin-
adapted configurations (SACs) instead, which conveniently split $\mathbf{H}$
into singlet and triplet blocks. The spatial part of the singlet SACs are

$\displaystyle\Phi^{S}_{ij}(1,2)=\begin{cases}\psi_{i}(1)\psi_{j}(2)&i=j\\\
\frac{1}{\sqrt{2}}(\psi_{i}(1)\psi_{j}(2)+\psi_{j}(1)\psi_{i}(2))&i\neq
j\end{cases}$ (31)

and for the triplet

$\displaystyle\Phi^{T}_{ij}(1,2)=\frac{1}{\sqrt{2}}(\psi_{i}(1)\psi_{j}(2)-\psi_{j}(1)\psi_{i}(2)),$
(32)

where the latter requires $i\neq 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:

$\displaystyle\Psi^{S}(1,2)=\sum_{i\leq
j}a^{S}_{ij}\Phi^{S}_{ij}(1,2)\,S(1,2)$ (33)
$\displaystyle\Psi^{T}(1,2)=\sum_{i<j}a^{T}_{ij}\Phi^{T}_{ij}(1,2)\,T(1,2)\ $
(34)

Note that using Full CI, it does not matter if we employ solutions of Eq. (13)
or Eq. (23).

####  III.3.2 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

$\displaystyle h_{ij}$
$\displaystyle=\braket{\psi_{i}}{\mathbf{h}}{\psi_{j}}=\sum_{\mu,\nu}\sum_{a,b}{C^{\mu}_{a,i}}^{\ast}C^{\nu}_{b,j}h^{\mu\nu}_{ab}$
(35) $\displaystyle(ij|kl)$
$\displaystyle=\sum_{\mu,\nu}\sum_{a,b,c,d}{C^{\mu}_{a,i}}^{\ast}C^{\mu}_{b,j}{C^{\nu}_{c,k}}^{\ast}C^{\nu}_{d,l}(ab|cd)_{\mu\nu}.$
(36)

Performing the transformation of the four-index tensor one index at a time
reduces the complexity to $\mathcal{O}(N_{\rm basis}^{5})$. With this, we can
easily compute the the SAC matrix elements, which in the triplet case they are
given by

$\displaystyle\begin{split}\mathbf{H}^{T}=\braket{\Phi^{T}_{ij}}{\mathbf{H}}{\Phi^{T}_{kl}}=&\delta_{jl}h_{ik}-\delta_{jk}h_{il}-\delta_{il}h_{jk}+\delta_{ik}h_{jl}\\\
&+(ik|jl)-(il|jk)\end{split}$ (37)

The singlet elements $\mathbf{H}^{S}$ are a little more complex due to the
cases $i=j$ in Eq. (31). By solving the eigenvalue problem

$\displaystyle\mathbf{H}^{S/T}\vec{a}_{n}^{S/T}=E_{n}^{S/T}\vec{a}_{n}^{S/T},$
(38)

we find the singlet/triplet $E_{n}^{S/T}$ energies as well as the
corresponding coefficients $a^{S/T}_{ij}$ (Eq. (33) /Eq. (34)) of the $n$th
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

$\displaystyle\rho(\vec{r}_{1})=\int|\Psi^{S}(\vec{r}_{1},\vec{r}_{2})|^{2}d^{3}\vec{r}_{2}$
(39)

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_{\rm RHF}=$-98.50\text{\,}\mathrm{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_{\rm
CI}=$-107.03\text{\,}\mathrm{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 statesGamble  _et
al._ (2015). The resulting energy is lower at $E_{\rm
CI}=$-109.64\text{\,}\mathrm{meV}$$. The energy difference between RHF and
Full CI is called correlation energy and notably accounts for
$10\text{\,}\mathrm{\char 37\relax}$ 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^{\mu}_{a}(\vec{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 anisotropicVoisin  _et al._ (2020); a $+z$-valley
orbital centered at the origin has the form:

$\displaystyle\text{STO}(\vec{r};\mathcal{N},n_{x},n_{y},n_{z},\alpha_{\perp},\alpha_{\parallel})=\mathcal{N}x^{n_{x}}y^{n_{y}}z^{n_{z}}e^{-\sqrt{\alpha_{\perp}(x^{2}+y^{2})+\alpha_{\parallel}z^{2}}},$
(40)

where $\mathcal{N}$ is a normalization factor, the $n_{i}$ are the exponents
of polynomial prefactors and the $\alpha$ 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 analyticallyGamble  _et al._ (2021). The GTO
analogue to the STO above has the following shape

$\displaystyle\text{GTO}(\vec{r};\mathcal{N},n_{x},n_{y},n_{z},\alpha_{\perp},\alpha_{\parallel})=\mathcal{N}x^{n_{x}}y^{n_{y}}z^{n_{z}}e^{-\alpha_{\perp}(x^{2}+y^{2})-\alpha_{\parallel}z^{2}}.$
(41)

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 accuracyJensen (2017). This
makes the transformation in Eq. (36), which scales as $\mathcal{O}(N_{\rm
basis}^{5})$, 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-$n$G basis function then has the following
form

$\displaystyle\text{STO-$n$G}(\vec{r};\mathcal{N},n_{x},n_{y},n_{z},\alpha_{\perp},\alpha_{\parallel})=\mathcal{N}\sum_{i=1}^{n}\text{GTO}(\vec{r};\mathcal{N}_{i},n_{x},n_{y},n_{z},\alpha_{\perp}\cdot\beta_{i},\alpha_{\parallel}\cdot\beta_{i}),$
(42)

where the $\mathcal{N}_{i}$ and $\beta_{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.

Table 1: Small basis set parameters resulting in the energies on the right.
The last orbital is added to fit the $D^{-}$ energy.

Basis Parameters

# | $(n_{x},n_{y},n_{z})$ | $\alpha_{\perp}({\rm nm}^{-2})$ | $\alpha_{\parallel}({\rm nm}^{-2})$
---|---|---|---
1 | $(0,0,0)$ | 0.246 | 0.950
2 | $(0,0,0)$ | 2.408 | 5.458
3 | $(0,0,0)$ | 0.0944 | 0.478

Central Cell Parameters

$A_{0}$ |  $-1.395\text{\,}\mathrm{meV}$
---|---
$A_{1}$ |  $-2717.0\text{\,}\mathrm{meV}$
$a$ |  $0.127\text{\,}\mathrm{nm}$
$b$ |  $0.194\text{\,}\mathrm{nm}$
$c$ |  $0.0972\text{\,}\mathrm{nm}$

Energy Levels

Having decided on the style of basis, we need to find a set of STO-$n$Gs that
models all features of the two-donor system with sufficient accuracy. We
follow Gamble et al.Gamble  _et al._ (2015) 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. 1 together with the orbitals #1 to #2 produce the D0 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 T2 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. 2
achieve good accuracy for all D0 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.

Table 2: Large basis set parameters resulting in the energies on the right.
The last orbital is added to fit the $D^{-}$ energy.

Basis Parameters

# | $(n_{x},n_{y},n_{z})$ | $\alpha_{\perp}({\rm nm}^{-2})$ | $\alpha_{\parallel}({\rm nm}^{-2})$
---|---|---|---
1 | $(0,0,0)$ | 0.227 | 0.864
2 | $(0,0,0)$ | 1.121 | 3.406
3 | $(2,0,0)$ | 17.04 | 103.9
4 | $(0,2,0)$ | 17.04 | 103.9
5 | $(0,0,2)$ | 17.04 | 103.9
6 | $(0,0,0)$ | 0.0695 | 0.302

Central Cell Parameters

$A_{0}$ |  $-1.626\text{\,}\mathrm{meV}$
---|---
$A_{1}$ |  $-2330.8\text{\,}\mathrm{meV}$
$a$ |  $0.140\text{\,}\mathrm{nm}$
$b$ |  $0.210\text{\,}\mathrm{nm}$
$c$ |  $0.102\text{\,}\mathrm{nm}$

Energy Levels

Figure 3: Basis Set analysis. (a) Direct comparison of the Large Basis (LB in
Tab. 2) and Small Basis (SB in Tab. 1) results as a function of donor
separation along three crystal directions. The first row shows the lowest
singlet (S) and triplet (T) energies for the LB and SB results as a grey
shadow. The second shows the energy difference $\Delta E$ between the LB and
SB for the S and T respectively. Both converge to the same value, i.e .
$\Delta E=0$ for large $d$. The third is the resulting difference in exchange
interaction strength on a log scale. (b) The first row shows a $y=z=0$ line
cut of the LB (full) and SB (dashed) anisotropic basis functions along (left)
and perpendicular (right) to the valley orientation. For the SB the blue, red
and yellow curves belong to the orbitals #1, #2 and #3 in Tab. 1, while for
the LB they belong to the orbitals #1, #2 and #6 in Tab. 2. Note that for the
LB the $n_{z}=n_{y}=2$ states vanish in the chosen line cut and only orbital
#3 is shown in purple. The contributions of these basis states (Eq. (44)) to
the lowest singlet and triplet states as a function of donor separation along
[100] are shown in the bottom panels. Here we use the same colour code as in
the top panels, e.g. the contribution of the state function visualized by a
full blue line in the top panels, belong to the full blue line curves in the
bottom panels. However, the purple curve accounts for all $2s$-orbitals.

In the following, we will test our two basis sets on two exchange-coupled
donors and compare the results to literature to get a better understanding of
their limitations. The first row of Fig. 3a shows the LB results for lowest
singlet and triplet energy states as a function of donor separation along the
[100], [110] and [111] crystal directions, in comparison to the SB results as
a gray shade. The differences between both basis sets are very subtle. For
both sets, the energies approach twice the energy of an isolated donor for
large distances, while the singlet-triplet splitting decays. The second row
shows the difference $\Delta E$ between the singlet/triplet energies
calculated using the SB and the LB, while the third row compares the exchange
interaction strength

$\displaystyle J=E^{T}_{0}-E^{S}_{0},$ (43)

on a logarithmic scale. The sensitivity to the choice of basis set decreases
with distance and the singlet and triplet show a similar trend in sensitivity,
i.e. amplitude of $\Delta E$. As a function of distance, the singlet and
triplet alternate in being more sensitive, but converge to the same value, i.e
. $\Delta E=0$ for large $d$. The shape of the basis functions and especially
the overlap of the flanks determine these fine features.

The first row of Fig. 3b shows a $y=z=0$ line cut of the $+x$\- and
$+z$-valley LB (full lines) and SB (dashed lines) basis functions in the top
panel, i.e. the line cuts parallel and perpendicular to the valley direction.
The bottom panel gives their contributions to the lowest singlet and triplet
states for two donors separated by $d$ along the [100] crystal axis in the
same colorcode. The state contributions (SC) of the $a$th basis state are
computed as

$\displaystyle
SC_{a}(\Psi^{S/T})=\sum_{\mu}\int|\braket{F^{\mu}_{a}(\vec{r}_{1})\varphi_{\mu}(\vec{r}_{1})}{\Psi^{S/T}(\vec{r}_{1},\vec{r}_{2})}|^{2}d^{3}\vec{r}_{2},$
(44)

where $F^{\mu}_{a}$ is the $a$th basis state in the $\mu$th valley, which
differs from the $a$th 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\lessapprox$5\text{\,}\mathrm{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_{\rm
B}=$0.53\text{\,}\mathrm{\SIUnitSymbolAngstrom}$$. 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\text{\,}\mathrm{\SIUnitSymbolAngstrom}$$. 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 $\Delta 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
$\Delta E$ nicely converging to zero.

Figure 4: (a) Comparison between the small basis (SB) and large basis (LB)
Full CI exchange interaction strength results in frequency units as a function
of donor separation along the [100], [110] and [111] crystal directions. The
SB results are shown as a grey shadow. (b) Double logarithmic plot of the
difference in exchange interaction strength between SB and LB sets as a
function of the exchange strength. The diagonal dashed lines mark the $\Delta
J/J=100\%$ and $10\%$ regions.

Fig. 4a gives our Full CI multi-valley EMT results for the LB (colored) and SB
(grey shade). The curves show exponential decay for
$d>$5\text{\,}\mathrm{nm}$$ and the typical valley oscillations for crystal
orientations not in $\braket{100}$. Additionally,
$J_{[100]}(d)>J_{[110]}(d)>J_{[111]}(d)$ due to the effective mass anisotropy
(see Sec. V). 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\text{\,}\mathrm{nm}$$ both basis sets are equally viable. Fig. 4b shows
the difference $\Delta J$ as a function of $J$. The diagonal dashed lines mark
the $\Delta 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 supercomputerGad . A LB iteration typically costs
$120\text{\,}\mathrm{S}\mathrm{U}$ (‘Service Units’ or computer core hours),
while a SB step costs $4.7\text{\,}\mathrm{S}\mathrm{U}$. We here showed that
the SB is sufficient for $d\gtrapprox$5\text{\,}\mathrm{nm}$$ and we can use
it to perform exhaustive 3D iterations (see Fig. 6) including $N\approx 8500$
iterations run with approximately $40\text{\,}\mathrm{k}\mathrm{S}\mathrm{U}$.
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.

Figure 5: (a) Exchange interaction strength $J=E^{T}_{0}-E^{S}_{0}$ as a
function of donor separation along three crystal orientations marked by the
correspondingly colored arrow in the silicon unit cell in (b). The dashed
lines indicate the exchange interaction strength range required for high-
fidelity CROT gates. (c) The $+x$-valley contribution
$|F^{+x}(\vec{r})\varphi_{+x}(\vec{r})|^{2}$ to the full $A_{1}$ ground state
$|\psi(\vec{r})|^{2}$ in (d). They are plotted on a log scale on the surface
of a sphere with $r=$10\text{\,}\mathrm{nm}$$ around the donor. (e) Analysis
panel as function of distance for the crystal orientations in (b). The first
row gives the lowest three energy states, a singlet and two triplets. The
second/third row give the valley contributions computed via Eq. (46) for the
lowest singlet and triplet. At roughly $d\approx$4\text{\,}\mathrm{nm}$$, the
triplet has an avoided crossing and gains $T_{2}$ contributions. Figure 6:
Electron density Eq. (39) calculations for one donor placed at the origin and
one moved along [100] (x-axis) with the slider. Shown is a $z=0$ cutplane. The
radio-button selects between the electron density of the lowest singlet ((a)
and (b) offline), triplet ((c) and (d) offline) or their difference ((e) and
(f) offline. (e) averages over the $y$-axis and iterates over all donor
distances $d$ instead). Figure 7: Electron density Eq. (39) calculations for
one donor placed at the origin and the second moved along [110] with the
slider. Shown is a $z=0$ cutplane with $r_{\parallel}$ being the
$x$-$y$-bisector. The radio-button selects between the electron density of the
lowest singlet ((a) and (b) offline), triplet ((c) and (d) offline) or their
difference ((e) and (f) offline. (e) averages over the $r_{\perp}$-axis and
iterates over all donor distances $d$ instead). Figure 8: Electron density Eq.
(39) calculations for one donor placed at the origin and the second moved
along [111] with the slider. Shown is a cutplane defined by the axis joining
the two donors ($r_{\parallel}$) and their cross-product with the z-axis
$r_{\perp}$. The radio-button selects between the electron density of the
lowest singlet ((a) and (b) offline), triplet ((c) and (d) offline) or their
difference ((e) and (f) offline. (e) averages over the $r_{\perp}$-axis and
iterates over all donor distances $d$ instead).

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^{\mu}(\vec{r})$ (Eq. (1))Voisin  _et al._ (2020) as
outlined in the following. Each $F^{\mu}(\vec{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}(\vec{r})\varphi_{+x}(\vec{r})|^{2}$ on the surface of a sphere with
$r=$10\text{\,}\mathrm{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\approx$4\text{\,}\mathrm{nm}$$, especially for [100]. At
this point the two lowest spin triplet states have an avoided crossingKlymenko
and Remacle (2014). 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 Dehollain  _et al._ (2014).

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_{1}$-$T_{2}$\- triplet
(red). Along [110] and [111], we can additionally observe energy oscillations.
The coupling between valley states depends on the symmetry of the system. The
inter-donor valley-orbit coupling strength (Eq. (3)) crucially depends on the
overlap of the Bloch functions on the individual donors. The plane wave
component of the Bloch functions is not lattice periodic and results in
interference patterns when the donors are moved relative to each other. This
causes the energy oscillations for lower symmetry cubic crystal orientations,
i.e. for all directions but the the $\braket{100}$.

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\text{\,}\mathrm{MHz}$<J<$10\text{\,}\mathrm{MHz}$$ Hill (2007); Kalra
_et al._ (2014), marked by the dashed lines in Fig. 5a. Depending on the
crystal direction this condition can be met starting from $d_{\rm
min}=$10.34\text{\,}\mathrm{nm}$$ along [111] up to $d_{\rm
max}=$23.892\text{\,}\mathrm{nm}$$ along [100]. Along [100] and [110], high-
fidelity CROT gates are possible for donors placed over a range of
$5\text{\,}\mathrm{nm}$. However, along [111] we find that the amplitude of
the valley oscillations is of the order of the allowed exchange interaction
strength range, and requires for precise donor placement. There is no distance
range in which high-fidelity CROT gates are possible independent of the
crystal orientation.

To analyse the valley contributions we define state projectors that have the
valley structure as defined in Eq. (4) to Eq. (6). An $A_{1}$-like state of a
donor at position $\vec{r}_{b}$ has components

$\displaystyle\psi_{A_{1}}^{b}(\vec{r}_{1})=\sum_{\mu}\sum_{a}\mathcal{N}A_{1}(\mu)\delta(\vec{r}_{a}-\vec{r}_{b})e^{i\vec{k}_{\mu}\vec{r}_{b}}F^{\mu}_{a}(\vec{r}_{1})\varphi_{\mu}(\vec{r}_{1})$
(45)

in analogy to Eq. (11) , where $a$ iterates the basis states in the $\mu$th
valley. $\mathcal{N}$ is a normalization constant, $A_{1}(\mu)$ gives the
relative valley amplitudes in Eq. (4) and $\vec{r}_{a}$ is the center of the
$a$th basis state. We compute the $A_{1}$-valley contribution

$\displaystyle
A_{1}(\Psi^{S/T})=\sum_{b}\int|\braket{\psi_{A_{1}}^{b}(\vec{r}_{1})}{\Psi^{S/T}(\vec{r}_{1},\vec{r}_{2})}|^{2}d^{3}\vec{r}_{2},$
(46)

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)). At the point where the lowest triplet has an
avoided crossing with the first excited triplet branch, it becomes roughly
half $T_{2}$ like. Going along [111], favours no axis in particular, but lifts
the directional symmetry along $x$, $y$ and $z$, introducing $T_{2}$
contributions (Eq. (5)), while leaving $E$ states untouched. For short
distances the triplet oscillates between large and small $T_{2}$ contributions
when crossing with the first excited state. Along [110] the $z$-valleys have a
strong overlap resulting in $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.

Figure 9: Exhaustive exchange interaction strength evaluation of two exchange-
coupled donors in various spatial configurations. One donor is position in the
center of the sphere, the second donor at the lattice site closest to each
point on the surface. We show the resulting exchange interaction strength on a
log scale. It is highly oscillatory and reproduces the anisotropy of the
valley-orbit basis.

Fig. 6, 7 and 8 give the electron densities Eq. (39) for a donor placed at the
origin and a second donor moved along the [100], [110] and [111] crystal axis.
The cut-planes are defined by the axis joining the two donors and their cross-
product with the z-axis. The slider moves the second donor, i.e. varies the
inter-donor distance. The radio button changes between the electron density of
the singlet, triplet or their difference. For large distances, they are
predominately $A_{1}$-like and resemble the one-electron ground stateGamble
_et al._ (2015). Additionally, the bonding singlet has more support in the
region between the donors, than the anti-bonding triplet. Along [110] and
[111], the amplitude of this discrepancy oscillates as a function of
distances, due to the changing overlap of the Bloch functions at each donor
position. We can observe the anisotropic structure of the wave functions on
each donor that decay very fast along [111] and [110], while being wide along
[100]. 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\approx$4\text{\,}\mathrm{nm}$$, the triplet state becomes
more $T_{2}$-like which changes the lattice periodic structure of $\rho_{T}$.
The bonding/anti-bonding structure of the $\rho_{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\text{\,}\mathrm{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.

##  VI Comparison to Literature

Fig. 10 compares our exchange results to EMT implementations in the
literature. The results of the cited works are shown in colour, while our
results are shown as a grey shadow for comparison. Additionally, we include an
experimental data point (purple) from Ref. González-Zalba  _et al._ (2014).

Early work by Koiller et al.Koiller  _et al._ (2001) and Wellard et al.Wellard
_et al._ (2003) 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.Wellard  _et al._ (2003) are fairly consistent with our
results. In both cases, the $J_{[110]}(d)$ and $J_{[111]}(d)$ levels have the
same magnitude, but lie below the $J_{[100]}(d)$. In comparison, we found that
$J_{[100]}(d)>J_{[110]}(d)>J_{[111]}(d)$.

Still using Heitler-London, Wellard et al.Wellard and Hollenberg (2005) and
Pica et al.Pica  _et al._ (2014b) included isotropic CCCs from Ref. Pantelides
and Sah (1974) and Ning and Sah (1971) 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.Wellard and Hollenberg (2005) in Fig. 10c are consistent with ours, while
the exchange in Fig. 10d is over an order of magnitude smaller.

Gonzalez-Zalba et al.González-Zalba  _et al._ (2014) and Saraiva et al.Saraiva
_et al._ (2015) 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]}(d)\approx J_{[110]}(d)\approx J_{[111]}(d)$. 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 bulk-screened
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.Wu and Fisher (2020) used an EMT including VOC and a CCC in the style
of Saraiva et al.Saraiva  _et al._ (2015). 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.2. Their [110]
results are comparable to our [111] calculations. Note that the exchange seems
to saturate at large distances.

Figure 10: EMT exchange interaction strength calculations in literature as a
function of distance along the [100], [110] and [111] crystal directions. The
results presented in Sec. V are shown as a grey shadow for easy comparison.
Additionally, we include an experimental data point for
$d=2.3\pm$0.5\text{\,}\mathrm{nm}$$ from Ref. González-Zalba  _et al._ (2014)
in purple (a), (b) Ref. Koiller  _et al._ (2001) and Wellard  _et al._ (2003)
used a Kohn-Luttinger EMT with a simple bulk-screened Coulomb potential. The
exchange is calculated with the Heitler-London method. (c), (d) Ref. Wellard
and Hollenberg (2005) and Pica  _et al._ (2014b) also used Heitler-London, but
included VOC and a CCC. (e) Gonzalez-Zalba et al.González-Zalba  _et al._
(2014) and Saraiva et al.Saraiva  _et al._ (2015) implemented a Full CI code.
Their EMT used an isotropic effective mass and orbitals, neglected VOC and
used a CCC specific to the $A_{1}$ ground state. (f) Ref. Wu and Fisher (2020)
included the effective mass anisotropy, VOC and a CCC. Their GTO basis however
is isotropic, the exchange is calculated via HF.

##  VII Electric Tunability and SWAP Gates

Figure 11: Analysis of the electrical tunability of two phosphorus donors
positioned at distance $d=$20.09\text{\,}\mathrm{nm}$$ along the [100] crystal
axis. An electric field $E_{x}$ is applied along the axis joining the donors
and we vary the detuning $\varepsilon=E_{x}d$. The first/second row looks at
the singlet/triplet states. (a), (e) The energies of the singlet (a) or
triplet (e) subspace as a function of detuning. In each panel, the energies of
the spin other subspace are shown in grey. The lowest energy states analyzed
in the following panels are marked in blue. (b), (f) The absolute difference
in probability of finding an electron on the left or right donor. (c), (g) The
valley contributions computed via Eq. (46). Upon ionization the triplet
becomes $A_{1}$-$T_{2}$-like. (d), (h) Contributions of the different basis
states computed via Eq. (44). Figure 12: Viability analysis of a SWAP gate
performed on two tunable phosphorus donors positioned along [100]. (a)
Exchange interaction strength as a function of detuning $\varepsilon$ for
distances between $10\text{\,}\mathrm{nm}$ and $25.5\text{\,}\mathrm{nm}$. $J$
increases with field until the system is ionized and $J$ saturates. The
amplitude of the SWAP operations are given by $F=J^{2}/(J^{2}+\Delta
B_{z}^{2})$ for anitparallel nuclear spins $\Uparrow\Downarrow$ with $\Delta
B_{z}=$120\text{\,}\mathrm{MHz}$$ in (b) and parallel nuclear spins
$\Uparrow\Uparrow$ with $\Delta B_{z}=$1\text{\,}\mathrm{MHz}$$ in (d). For
high-fidelity gates, we need to tune $F$ from $0.01$ (lower dashed line) to
$0.9$ (upper dashed line). The corresponding exchange tuning ranges are
indicated in (a). (c)/(e) For viable distances, we evaluate the quality factor
in Eq. (50) (number of coherent oscillations) at the $F=0.9$ points. Figure
13: Viability analysis of a SWAP gate performed on two tunable phosphorus
donors positioned along [110]. (a) Exchange interaction strength as a function
of detuning $\varepsilon$ for distances between $10\text{\,}\mathrm{nm}$ and
$25\text{\,}\mathrm{nm}$. $J$ increases with field until the system is ionized
and $J$ saturates. The amplitude of the SWAP operations are given by
$F=J^{2}/(J^{2}+\Delta B_{z}^{2})$ for $\Delta
B_{z}=$120\text{\,}\mathrm{MHz}$$ in (b) and $\Delta
B_{z}=$1\text{\,}\mathrm{MHz}$$ in (d). For high-fidelity gates, we need to
tune $F$ from $0.01$ (lower dashed line) to $0.9$ (upper dashed line). For
viable distances, we evaluate the corresponding quality factor Eq. (50)
(number of coherent oscillations) at the $F=0.9$ points in (c) and (e).

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 Kane (1998); Hollenberg  _et al._
(2006); Hill (2007). We add an electric field potential to our EMT Hamiltonian
in Eq. (12)

$\displaystyle U_{E}(\vec{r})=\vec{E}\cdot\vec{r},$ (47)

where $\vec{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
$\epsilon=|\vec{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\text{\,}\mathrm{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
$\ket{\uparrow\downarrow}$ are described by the Rabi formulaRabi (1937)

$\displaystyle
P_{\ket{\uparrow\downarrow}\rightarrow\ket{\downarrow\uparrow}}(t)=\frac{J^{2}}{J^{2}+\Delta
B_{z}^{2}}\sin^{2}\left(\sqrt{J^{2}+\Delta B_{z}^{2}}\frac{t}{2}\right),$ (48)

where $J$ is the exchange interaction strength and $\Delta B_{z}$ the
longitudinal magnetic field gradient across the two electrons. Tuning $J$
relative to $\Delta 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\approx 120$ MHz, which in a gated nanoscale device
can vary from one donor to the next by an order $\Delta A\approx 1$ MHz. This
allows the introduction of an intrinsic $\Delta B_{z}$, switchable by
controlling the state of the two nuclear spins Kalra  _et al._ (2014); Mądzik
_et al._ (2020).

If the nuclei are parallel, a small $\Delta B_{z}=\Delta
A\approx$1\text{\,}\mathrm{MHz}$$ is created by the difference in the
hyperfine interaction between the two donors. A much larger $\Delta
B_{z}=\bar{A}\approx$120\text{\,}\mathrm{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\text{\,}\mathrm{nm}$ and $25.5\text{\,}\mathrm{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 $\epsilon$. 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}+\Delta B_{z}^{2})$ for $\Delta
B_{z}=$120\text{\,}\mathrm{MHz}$$ and $\Delta B_{z}=$1\text{\,}\mathrm{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_{\pi}=\pi/\sqrt{J^{2}+\Delta B_{z}^{2}}$ 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 byDial  _et al._ (2013)

$\displaystyle T_{2}^{*}=\frac{1}{\sqrt{2}\pi(dJ/d\epsilon)\epsilon_{rms}},$
(49)

where we take
$\epsilon_{rms}=$24\text{\times}{10}^{-3}\text{\,}\mathrm{GHz}\text{\,}{\mathrm{nm}}^{-1}$\times
d$ as our noise amplitudeFreeman  _et al._ (2016); Thorgrimsson  _et al._
(2017); Harvey-Collard  _et al._ (2017). We then define the quality factor as
the number of SWAP oscillations we can perform in $T_{2}^{*}$

$\displaystyle Q=\frac{T_{2}^{*}}{T_{\pi}}=\sqrt{\frac{J^{2}+\Delta
B_{z}^{2}}{2(dJ/d\epsilon)^{2}\epsilon_{rms}^{2}}}.$ (50)

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
$\gg$1\text{\,}\mathrm{\SIUnitSymbolMicro s}$$ pulses to avoid spectral
broadening beyond the value of $J$, whereas antiparallel spins allow much
faster gate operations, on the order of $100\text{\,}\mathrm{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. 1 to be
sufficient for donor separations $d\gtrapprox$5\text{\,}\mathrm{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\text{\,}\mathrm{nm}$ and $24\text{\,}\mathrm{nm}$ depending on the crystal
direction. Along [100] and [110] the donor placement may straggle by
$5\text{\,}\mathrm{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.

##  Appendix A STO-nG fits

Figure 14: Visulatization of the STO-$n$G fits (see Eq. (42)).

##  Appendix B Carbon footprint

Table 3: Carbon emissions involved in all calculations run in this project. We included all calculations, even if the results are not shown in this publication. Estimations follow the standardised reporting table by Scientific CO2nductco (2).  Numerical Simulations |
---|---
Total Kernel Hours (h) | 181130
Thermal Design Power Per Kernel (W) | 5.75
Total Energy Consumption Simulations (kWh) | 1041.5
Average Emission Of CO2 In New South WalesDEE (2020) (kg/kWh) | 0.81
Total CO2-Emission For Numerical Simulations (kg) | 843.6
Were The Emissions Offset? | No
Transport |
Total CO2-Emission for Transport (kg) | 0
Were The Emissions Offset? | n/a
Total CO2-Emission (kg) | 843.6

###### Acknowledgements.

The research at UNSW was funded by the Australian Research Council Centre of
Excellence for Quantum Computation and Communication Technology (Grant No.
CE170100012) and the US Army Research Office (Contract No. W911NF-17-1-0200).
This research/project was undertaken with the assistance of resources and
services from the National Computational Infrastructure (NCI), which is
supported by the Australian Government. The research at Sandia National
Laboratories was supported by the Laboratory Directed Research and Development
program, Project 213048. Sandia National Laboratories is a multi-missions
laboratory managed and operated by National Technology and Engineering
Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International
Inc., for the National Nuclear Security Administration of the US Department of
Energy under contract DE-NA0003525. The views and conclusions contained in
this document are those of the authors and should not be interpreted as
representing the official policies, either expressed or implied, of the ARO or
the US Government. The US Government is authorized to reproduce and distribute
reprints for government purposes notwithstanding any copyright notation
herein.

## References

  * Muhonen  _et al._ (2014) J. T. Muhonen, J. P. Dehollain, A. Laucht, F. E. Hudson, R. Kalra, T. Sekiguchi, K. M. Itoh, D. N. Jamieson, J. C. McCallum, A. S. Dzurak, _et al._ , Nature nanotechnology 9, 986 (2014).
  * Morello  _et al._ (2010) A. Morello, J. J. Pla, F. A. Zwanenburg, K. W. Chan, K. Y. Tan, H. Huebl, M. Möttönen, C. D. Nugroho, C. Yang, J. A. van Donkelaar, _et al._ , Nature 467, 687 (2010).
  * Pla  _et al._ (2012) J. J. Pla, K. Y. Tan, J. P. Dehollain, W. H. Lim, J. J. Morton, D. N. Jamieson, A. S. Dzurak, and A. Morello, Nature 489, 541 (2012).
  * Muhonen  _et al._ (2015) J. Muhonen, A. Laucht, S. Simmons, J. Dehollain, R. Kalra, F. Hudson, S. Freer, K. M. Itoh, D. Jamieson, J. McCallum, _et al._ , Journal of Physics: Condensed Matter 27, 154205 (2015).
  * Dehollain  _et al._ (2016) J. P. Dehollain, J. T. Muhonen, R. Blume-Kohout, K. M. Rudinger, J. K. Gamble, E. Nielsen, A. Laucht, S. Simmons, R. Kalra, A. S. Dzurak, _et al._ , New Journal of Physics 18, 103018 (2016).
  * Kane (1998) B. E. Kane, nature 393, 133 (1998).
  * Hollenberg  _et al._ (2006) L. Hollenberg, A. Greentree, A. Fowler, and C. Wellard, Physical Review B 74, 045311 (2006).
  * Hill (2007) C. D. Hill, Physical review letters 98, 180501 (2007).
  * Kalra  _et al._ (2014) R. Kalra, A. Laucht, C. D. Hill, and A. Morello, Physical Review X 4, 021044 (2014).
  * He  _et al._ (2019) Y. He, S. Gorman, D. Keith, L. Kranz, J. Keizer, and M. Simmons, Nature 571, 371 (2019).
  * Mądzik  _et al._ (2020) M. T. Mądzik, A. Laucht, F. E. Hudson, A. M. Jakob, B. C. Johnson, D. N. Jamieson, K. M. Itoh, A. S. Dzurak, and A. Morello, arXiv preprint arXiv:2006.04483 (2020).
  * Kohn and Luttinger (1955) W. Kohn and J. Luttinger, Physical Review 98, 915 (1955).
  * Luttinger and Kohn (1955) J. M. Luttinger and W. Kohn, Physical Review 97, 869 (1955).
  * Ning and Sah (1971) T. H. Ning and C.-T. Sah, Physical Review B 4, 3468 (1971).
  * Pantelides and Sah (1974) S. T. Pantelides and C.-T. Sah, Physical Review B 10, 621 (1974).
  * Shindo and Nara (1976) K. Shindo and H. Nara, Journal of the Physical Society of Japan 40, 1640 (1976).
  * Hui (2013) H. Hui, Solid state communications 154, 19 (2013).
  * Saraiva  _et al._ (2015) A. Saraiva, A. Baena, M. Calderón, and B. Koiller, Journal of Physics: Condensed Matter 27, 154208 (2015).
  * Calderón  _et al._ (2009) M. J. Calderón, A. Saraiva, B. Koiller, and S. Das Sarma, Journal of applied physics 105, 122410 (2009).
  * Baena  _et al._ (2012) A. Baena, A. Saraiva, B. Koiller, and M. Calderón, Physical Review B 86, 035317 (2012).
  * Hu  _et al._ (2005) X. Hu, B. Koiller, and S. D. Sarma, Physical Review B 71, 235332 (2005).
  * Klymenko and Remacle (2014) M. Klymenko and F. Remacle, Journal of Physics: Condensed Matter 26, 065302 (2014).
  * Gamble  _et al._ (2015) J. K. Gamble, N. T. Jacobson, E. Nielsen, A. D. Baczewski, J. E. Moussa, I. Montaño, and R. P. Muller, Physical Review B 91, 235318 (2015).
  * Klymenko  _et al._ (2017) M. Klymenko, S. Rogge, and F. Remacle, Physical Review B 95, 205301 (2017).
  * Saraiva  _et al._ (2011) A. Saraiva, M. Calderón, R. B. Capaz, X. Hu, S. D. Sarma, and B. Koiller, Physical Review B 84, 155320 (2011).
  * Gamble  _et al._ (2016) J. K. Gamble, P. Harvey-Collard, N. T. Jacobson, A. D. Baczewski, E. Nielsen, L. Maurer, I. Montaño, M. Rudolph, M. Carroll, C. Yang, _et al._ , Applied Physics Letters 109, 253101 (2016).
  * Friesen (2005) M. Friesen, Physical review letters 94, 186403 (2005).
  * Debernardi  _et al._ (2006) A. Debernardi, A. Baldereschi, and M. Fanciulli, Physical Review B 74, 035202 (2006).
  * Pica  _et al._ (2014a) G. Pica, G. Wolfowicz, M. Urdampilleta, M. L. Thewalt, H. Riemann, N. V. Abrosimov, P. Becker, H.-J. Pohl, J. J. Morton, R. N. Bhatt, _et al._ , Physical Review B 90, 195204 (2014a).
  * Szabo and Ostlund (2012) A. Szabo and N. S. Ostlund,  _Modern quantum chemistry: introduction to advanced electronic structure theory_ (Courier Corporation, 2012).
  * Jensen (2017) F. Jensen,  _Introduction to computational chemistry_ (John wiley & sons, 2017).
  * Koiller  _et al._ (2001) B. Koiller, X. Hu, and S. D. Sarma, Physical review letters 88, 027903 (2001).
  * Wellard  _et al._ (2003) C. Wellard, L. Hollenberg, F. Parisoli, L. Kettle, H.-S. Goan, J. McIntosh, and D. Jamieson, Physical Review B 68, 195209 (2003).
  * Wellard and Hollenberg (2005) C. Wellard and L. Hollenberg, Physical Review B 72, 085202 (2005).
  * Pica  _et al._ (2014b) G. Pica, B. Lovett, R. N. Bhatt, and S. A. Lyon, Physical Review B 89, 235306 (2014b).
  * Wu and Fisher (2020) W. Wu and A. Fisher, arXiv preprint arXiv:2006.16127 (2020).
  * Kettle  _et al._ (2006) L. Kettle, H.-S. Goan, and S. C. Smith, Physical Review B 73, 115205 (2006).
  * González-Zalba  _et al._ (2014) M. F. González-Zalba, A. Saraiva, M. J. Calderón, D. Heiss, B. Koiller, and A. J. Ferguson, Nano letters 14, 5672 (2014).
  * Narita  _et al._ (1982) S.-i. Narita, T. Shinbashi, and M. Kobayashi, Journal of the Physical Society of Japan 51, 2186 (1982).
  * Grimmeiss  _et al._ (1982) H. Grimmeiss, E. Janzén, and K. Larsson, Physical Review B 25, 2627 (1982).
  * Voisin  _et al._ (2020) B. Voisin, J. Bocquel, A. Tankasala, M. Usman, J. Salfi, R. Rahman, M. Simmons, L. Hollenberg, and S. Rogge, Nature Communications 11, 1 (2020).
  * Gamble  _et al._ (2021) J. Gamble _et al._ ,  _Computational Materials, Chemistry, and Biochemistry: From Bold Initiatives to the Last Mile_ , edited by S. Shankar _et al._ (Springer, London, UK, 2021).
  * (43) “National computational infrastructure,” http://nci.org.au/.
  * Dehollain  _et al._ (2014) J. P. Dehollain, J. T. Muhonen, K. Y. Tan, A. Saraiva, D. N. Jamieson, A. S. Dzurak, and A. Morello, Physical review letters 112, 236801 (2014).
  * Rabi (1937) I. I. Rabi, Physical Review 51, 652 (1937).
  * Dial  _et al._ (2013) O. Dial, M. D. Shulman, S. P. Harvey, H. Bluhm, V. Umansky, and A. Yacoby, Physical review letters 110, 146804 (2013).
  * Freeman  _et al._ (2016) B. M. Freeman, J. S. Schoenfield, and H. Jiang, Applied Physics Letters 108, 253108 (2016).
  * Thorgrimsson  _et al._ (2017) B. Thorgrimsson, D. Kim, Y.-C. Yang, L. Smith, C. Simmons, D. R. Ward, R. H. Foote, J. Corrigan, D. Savage, M. Lagally, _et al._ , npj Quantum Information 3, 1 (2017).
  * Harvey-Collard  _et al._ (2017) P. Harvey-Collard, N. T. Jacobson, M. Rudolph, J. Dominguez, G. A. Ten Eyck, J. R. Wendt, T. Pluym, J. K. Gamble, M. P. Lilly, M. Pioro-Ladrière, _et al._ , Nature communications 8, 1 (2017).
  * co (2) “Scientific co2nduct,” https://scientific-conduct.github.io.
  * DEE (2020) DEE, “National greenhouse accounts factors,” (2020).