  Last year, Salfi et al. made the first direct measurements of a donor wave
function and found extremely good theoretical agreement with atomistic
tight-binding [Salfi et al., Nat. Mater. 13, 605 (2014)]. Here, we show that
multi-valley effective mass theory, applied properly, does achieve close
agreement with tight-binding and hence gives reliable predictions. To
demonstrate this, we variationally solve the coupled six-valley Shindo-Nara
equations, including silicon's full Bloch functions. Surprisingly, we find that
including the full Bloch functions necessitates a tetrahedral, rather than
spherical, donor central cell correction to accurately reproduce the
experimental energy spectrum of a phosphorus impurity in silicon. We
cross-validate this method against atomistic tight-binding calculations,
showing that the two theories agree well for the calculation of donor-donor
tunnel coupling. Further, we benchmark our results by performing a statistical
uncertainty analysis, confirming that derived quantities such as the wave
function profile and tunnel couplings are robust with respect to variational
energy fluctuations. Finally, we apply this method to exhaustively enumerate
the tunnel coupling for all donor-donor configurations within a large search
volume, demonstrating conclusively that the tunnel coupling has no spatially
stable regions. Though this instability is problematic for reliably coupling
donor pairs for two-qubit operations, we identify specific target locations
where donor qubits can be placed with scanning tunneling microscopy technology
to achieve reliably large tunnel couplings.
