Theory of one and two donors in silicon

THEORY OF ONE AND TWO DONORS IN SILICON A. L. Saraiva,1, 2, ∗ A. Baena,2 M. J. Calderón,3 and Belita Koiller2 1 arXiv:1407.8224v1 [cond-mat.mes-hall] 30 Jul 2014 University of Wisconsin-Madison, Madison, Wisconsin 53706, USA 2 Instituto de Fisica, Universidade Federal do Rio de Janeiro, Caixa Postal 68528, 21941-972 Rio de Janeiro, Brazil 3 Instituto de Ciencia de Materiales de Madrid, ICMM-CSIC, Cantoblanco, E-28049 Madrid, Spain (Dated: August 1, 2014) We provide here a roadmap for modeling silicon nano-devices with one or two group V donors (D). 0 We discuss systems containing one or two electrons, that is, D0 , D− , D+ 2 and D2 centers. The impact of different levels of approximation is discussed. The most accurate instances – for which we provide quantitative results – are within multivalley effective mass including the central cell correction and a configuration interaction account of the electron-electron correlations. We also derive insightful, yet less accurate, analytical approximations and discuss their validity and limitations – in particular, for a donor pair, we discuss the single orbital LCAO method, the Hückel approximation and the Hubbard model. Finally we discuss the connection between these results and recent experiments on few dopant devices. I. INTRODUCTION In a letter to Peierls in 1931, Pauli stated that semiconductors “are a filthy mess” [1]. Curiously, the impurities and defects to which Pauli refers led the way to the revolution severely impacting industry and society over the last many decades, and still under way. Development and advances in semiconductor-based devices started less than 20 years after Pauli’s statement, as the transistor operated for the first time in 1947. This achievement prompted intensive research for a deeper understanding of semiconductors and the role of dopants. Currently, as “Moore’s Law” approaches its limit with devices reaching the atomic scale, quantum behavior of electrons and spins drive the advancement of the semiconductor research into the new fields of quantum electronics and spintronics [2– 4]. Doping allows controlling the sign and density of carriers in a semiconductor, a flexibility not achievable in conductors or in insulators. For the applications where current carriers perform the needed operations, donor electrons are promoted to the nearby conduction band and similarly holes from acceptors operate in the valence band. In this case the concentration of donors and acceptors and their binding energies are the only relevant quantities defining the behavior of the doped material. A macroscopic concentration of dopants leads to macroscopically observable currents or voltages, and at this point the quantum behavior becomes irrelevant. We revisit here the theoretical treatment of donors D in semiconductors from a single substitutional donor perspective, considering that its active electron(s) remains bound to the core, and similarly for donor pairs D2 , a scenario where full quantum behavior prevails. Special attention must be paid to the electronic bound state wavefunctions. This paper is organized as follows. In Sec. II we briefly review the effective mass theory (EMT) for shallow donors, highlighting its successes and limitations. In Sec. III, multivalley and anisotropy effects are analysed, details of the approach used are given, and the results for neutral donors are shown in Sec. IV. Section V reports on a donor with two bound electrons, constituting a D− negatively charged donor. The donor pair is studied in Secs. VI and VII for the one and two electron states respectively. Sec. VIII discusses simplified models for the two donor problem and makes connections with experiments. Finally, summary and conclusions are presented in Sec. IX. In order to keep the text reasonably self-contained, well established basic material and results are included with proper references provided. II. HYDROGENIC MODEL FOR DONORS IN SEMICONDUCTORS The single donor description is one of the simplest and most successful implementations of the EMT. It is based on the usual effective-medium assumptions, namely that: 1. the perturbation potential due to a substitutional donor varies slowly at the scale of the lattice constant of the host semiconductor; 2. the envelope functions of the bound states extend over long enough distances such that it is composed by a very narrow distribution of Bloch wavevectors k around the band minimum. Further assuming that the conduction band lower edge is non-degenerate and isotropic, one obtains the hydrogenic atom hamiltonian HH = − ~2 ∇2 e2 − ∗ 2m 4πǫr (1) with the electron mass substituted by the effective mass m∗ and the Coulomb potential screened by the static dielectric constant ǫ. This gives a hydrogenic atom of effective Bohr radius a∗ = a0 ǫ/m∗ and ground state binding energy E ∗ = EH m∗ /ǫ2 , with a0 = 0.053 nm and 2 GaAs:Si Binding Energies 30 ZnSe:F ZnSe:Ga ZnSe:Cl Experimental (meV) 25 ZnSe:Al 20 CdTe:Al, In 15 10 GaAs:S, Se, C, Si, Ge 5 InSb:Te 0 0 5 10 15 20 25 30 Effective Mass Theory - Hydrogenic (meV) FIG. 1: Comparison of the experimental donor binding energies for single valley materials (blue dots) with those calculated from the effective mass hydrogenic expression (black line). The agreement is very good mainly for the less bound donors. The ground state energies for donors in silicon (not shown here) are off the given scale and present a large variation (∼ 50%) as a function of the donor species, see Table I. The same spread of values is found for germanium [5]. With the N minima (valleys) at finite wavevectors, the symmetry of the wavefunction is significantly lower, and highly anisotropic effective masses are not uncommon – in silicon the longitudinal and transverse effective masses are mL = 0.916 me and mT =0.191 me at the N = 6 minima. The hydrogenic hamiltonian Eq. (1) may be adapted to account for the mass anisotropy, as discussed in Sec. III A. The mere reconciliation with the mass anisotropy still does not lead to results compatible with the experimental binding energies. The missing ingredient is the elusive valley-orbit coupling. In the absence of valley-orbit coupling, a wavefunction equivalent to Eq. (2) is obtained for each of the N valleys, leading to an N -fold degenerate ground state. The 1/r Coulomb potential singularity at r = 0 is clearly incompatible with the assumptions discussed in Sec. II – it actually mediates a finite coupling among the N valleys. This accounts for over 30% of the binding energy in silicon donors, as will be discussed in Sec. III B. We briefly present some of the most successful semiempirical approaches to describe these departures from the hydrogenic model, namely the first order perturbation intervalley coupling (Sec. III C) and the central cell correction (Sec. III D). A theoretical description of valleyorbit coupling from first principles is still lacking. A. EH = 13.6 eV the respective values for Hydrogen in vacuum. The wavefunction for the ground state electron is then the product of the hydrogenic envelope-function by the band-edge Bloch function Ψ= √ 1 πa ∗ ∗3 e−r/a eik·r uk (r), (2) The EMT provides a highly accurate description of shallow donors in single valley materials. It involves no fitting parameters, as the binding energy is dependent on the material only through ǫ and m∗ , leading to a value of E ∗ that is independent of the donor species. This is in fact clearly demonstrated experimentally for example in GaAs, where deviations from the hydrogenic expression are hard to detect [6], as illustrated in Fig. 1. In contrast, the binding energy of donors in Si and Ge is strongly dependent on the donor species [5]. III. MULTIVALLEY EFFECTS The case of multivalley semiconductors is not as simple. In a recent work [7], spatially resolved spectroscopy of isolated As donors in silicon was combined with multiband semi-empirical tight binding theory to study the valley interference. These results provide direct evidence for the rich valley structure predicted over 60 years ago by Kohn and Luttinger [8] (KL) and further outline the role of the environment on valley repopulation. Uncoupled anisotropic valleys The first successful approach to treat donors in Si was presented by KL [8] (which also set stronger formal basis for the EMT) and Kittel and Michel [9]. Taking initially the EMT assumption (ii) to be valid – namely that only small deviations from the wavevector at the band minimum kmin contribute to the wavefunction– means that any individual valley does not couple to the other five. In this way each valley is solved independently and the problem regains the hydrogenic simplicity. The band anisotropy leads to a specific effective mass hamiltonian for the different valleys. For instance, for the 2 valleys along z, it reads  2  ~2 ∂2 ∂ ~2 ∂ 2 Hz = − + + V (r). (3) − 2mT ∂x2 ∂y 2 2mL ∂z 2 The ground state is found variationally, choosing a trial (001) valley (aSi is the Si function for the kz = 0.85 a2π Si conventional lattice parameter) of the shape of a deformed 1s orbital with two radii a and b [10] φz (r) = √ 1 πa2 b e− √ (x2 +y 2 )/a2 +z 2 /b2 ik·r e uk (r). (4) For other valleys at {kµ }, µ = {x, −x, y, −y, z, −z}, the hamiltonians Hµ and the wavefunctions φµ are immediately obtained exchanging z by x or y accordingly. The deformed orbital Bohr radii a and b are taken as variational parameters to minimize the expectation value 3 of hφµ |Hµ |φµ i. For the hydrogenic potential V (r) = VH (r) = −e2 /4πǫSi r, the ground state is sixfold degenerate, and one recovers the KL result with energy EKL = −31.2 meV, and an anisotropy b/a = 0.56. B. Valley-orbit coupling The binding energy is significantly underestimated if the valley-orbit coupling is disregarded, as can be seen by comparison of EKL to the observed values for shallow donors in Si in Table I. The valley-orbit coupling breaks the sixfold degeneracy, leading to a significant reduction of the ground state energy. Even if the strength of this coupling is not known, group theory arguments give the appropriate superpo- sition of valley states induced when the spherically symmetric potential is disrupted by the tetrahedral crystal field potential. In other words, the irreducible representation associated with the ground state of a spherical potential with six-fold valley degeneracy (the 1s manifold) becomes reducible in the presence of a local tetrahedral crystal field: the sixfold degenerate level splits into states that have the symmetry of the different irreducible representations of the Td group. In this case, this leads to a singlet with A1 symmetry, a triplet with T2 symmetry and a doublet with E symmetry. Each of the six states in the 1s manifold correspond to a particular combination of the envelope-modulated Bloch functions from the six degenerate valleys φµ illustrated in Fig. 2. These states are given explicitly by [11] 1 ΨA1 (r) = √ [φx (r) + φ−x (r) + φy (r) + φ−y (r) + φz (r) + φ−z (r)] 6 1 ΨT2x (r) = √ [φx (r) − φ−x (r)] 2 1 ΨT2y (r) = √ [φy (r) − φ−y (r)] 2 1 ΨT2z (r) = √ [φz (r) − φ−z (r)] 2 1 ΨE z (r) = √ [φx (r) + φ−x (r) + φy (r) + φ−y (r) − 2φz (r) − 2φ−z (r)] 12 1 ΨE xy (r) = [φx (r) + φ−x (r) − φy (r) − φ−y (r)] . 2 The direct calculation of the valley-orbit coupling is not trivial, though. We are unaware of any successful description of the complete spectrum of the 1s manifold from first principles. This difficulty is due to the incomplete knowledge regarding the impurity potential. The picture of a 1/r point charge potential screened by a dielectric constant ǫSi is realistic only at large distances from the impurity site. Close to the nucleus, the dielectric response becomes inhomogeneous. Furthermore, the valence core shells for the impurity are different from the Si atoms, presenting a strong departure from the excess proton picture. Atomic species differ from each other at the core region, leading to the strong donor species dependence of the binding energy in silicon in contrast with other semiconductors. Early attempts of correcting the Coulomb potential adopting the ab initio dielectric function were modestly successful for the isovalent impurity P and were inaccurate for all other group V donors (for a complete discussion of early works, see Ref. [12]). Successful attempts involve empirical fitting ingredients [13– 15]. Some examples used in the literature to describe this departure from the point charge picture are shown (5) in Fig. 3. Two among the main strategies for describing the valley-orbit induced splitting are described next. C. First order perturbation theory: the intervalley coupling The EMT assumptions (i) and (ii) are clearly incompatible with the Coulomb potential, which is singular at the donor site r = 0, implying that different valley states may couple. The single-valley KL description neglects this singularity, i.e. H defined in Eq. (3) does not include the singular part of the Coulomb potential VH . We may then include the intervalley coupling as a perturbation to the KL description. The total Hamiltonian would be H = H + H ′ with H ′ accounting for the singular part of the potential around the core region r → 0, i.e. VH (r → 0) plus additional corrections as discusssed below. For simplicity in notation we are assuming an isotropic approximation (further discussed in Sec. IV) which implies that Hz = Hx = Hy = H. Now, the single-valley approximation implies 4 (a) Single Valley Approximation +kZ - kX - kY +kY With Valley-Orbit (6x) EKL (a) E (2x) T (3x) (b) TB (c) +kX - kZ 𝐴1 (1x) (b) A1|A1 y x x 2 rcc y Tx|Tx x x FIG. 2: (a) Representation in k-space of the six constant energy surfaces around the degenerate conduction band minima in Si. (b) Splitting of the sixfold degenerate conduction band minimum by a tetrahedral crystal environment combined with the singular donor potential. (c) Charge distribution for the A1 and Tx states along the x-axis (large) and at the (001) Si crystal plane (inset). Note that A1 has significant charge density at the donor site, while the Tx state has a node. The lower lines curving down give the plain screened Coulomb potential VH (solid line) and the central-cell corrected potential Vimp (dashed line) which is more attractive for distances from the donor nucleus smaller than ∼ rcc . Due to the larger charge density of the A1 state within rcc , its energy is well below the T2 and E states. hkµ |H|kν i = EKL δµ,ν , where |kµ i stands for φµ . The valley-coupling part, H ′ , is added and treated following standard degenerate perturbation theory, i.e., we solve the perturbed Hamiltonian H = H + H ′ restricted to the degenerate manifold basis {φµ }. It is easy to see by symmetry that hkµ |H ′ |kµ i is a negative constant (independent of µ) which we call (−δ), with this the matrix elements of the hamiltonian written in the {φµ } basis may be obtained hkµ |H|kµ i = hkµ |H|kµ i + hkµ |H ′ |kµ i = EKL − δ, hkµ |H|k−µ i = hkµ |H ′ |k−µ i = ∆k , hkµ |H|kν i = hkµ |H ′ |kν i = ∆⊥ , if ν 6= ±µ. (6) Direct calculation of the eigenvectors of the hamiltonian FIG. 3: Comparison between different models adopted in the literature for the central cell correction. The figure gives QP (r) such that VP (r) = −QP (r)VH (r) calculated for P in Si. The different curves correspond to (a) Ref. [15], (b) Vimp (r) as in Eq. (7), (c) Ref. [13], and the triangles give the approximation used in tight-binding [14]. The limit of a point charge QP = 1 is shown as a dotted horizontal line. Note that, for the tight-binding, QP (r) = 1 for all sites except the impurity site where Vimp is taken as a constant. In this case we indicate QP (r) = 0 in the figure, which eliminates the 1/r divergence at r → 0, but a constant U0 should be added. matrix yield a spectrum split in a singlet A1 at energy EKL − δ + ∆k + 4∆⊥ , a triplet T2 at EKL − δ − ∆k , and a doublet E at EKL − δ + ∆k − 2∆⊥ . A possible semi-empirical solution would be to treat the three matrix elements δ, ∆k and ∆⊥ as parameters chosen to reproduce the 1s manifold spectrum of each donor species. In Refs. [16,17] the values ∆k = −1.52 meV and ∆⊥ = −2.16 meV are suggested for P donors in Si. This choice reproduces the relative energy splitting among the different 1s manifold levels. An overall negative shift of δ = 4.22 meV for P is required in order to get the actual experimental energies. The small values of ∆k , ∆⊥ and δ with respect to EKL justify the perturbative treatment described here. The observed experimental energies can also be understood from central cell corrections, which take into account the incomplete screening of the positive donor potential, as detailed in Sec. III D. The adopted central cell correction involves a single, instead of three, empirical parameter characterizing each donor species. D. Effective pseudopotential: the central cell correction Another strategy to describe the donors beyond the single valley theory is to build a pseudopotential Vimp (r) that effectively mimics the complex environment of the crystal and the screening due to the valence electrons of the donor. The suffix imp stands for the impurity type, so imp=P, As, Sb or Bi. Multivalley semiconductors are more strongly subject to this inner core potential than single valley semiconductors because the valley degree of freedom allows for the electron to concentrate closer to the nucleus. The rationale for that is simple – each wavefunction φµ spreads in k-space around kµ in a form 5 determined by the Heisenberg uncertainty principle. The spread of the superposition of various φµ in k-space tends to be larger than on the non-degenerate case, so that the real space spread may be smaller without costing any kinetic energy. Thus, multivalley electronic states tend to concentrate around the nucleus, benefiting from its potential energy and therefore lowering its total energy. A consequence is that a larger portion of the ground state energy comes from the potential very close to the nucleus, which is precisely the region where the point charge model fails. Formally the coupling between valleys, hkµ |H ′ |kν i, could be obtained directly from the perturbation due to the donor breaking the effective mass assumptions, modeled above by the term H ′ . The full donor potential, Vimp (r), has been treated in the literature within several models, a few of which are given in Fig. 3. For example, a constant replacement for the screened Coulomb potential within a central cell region of radius R has been adopted in Ref. [18], namely, Vimp (r < R) = V0 , Vimp (r > R) = VH (r), with R taken as half the nearest neighbors distance in Si and V0 is a negative energy chosen to reproduce the experimental spectra. A conceptually similar scheme, widely adopted in tight-binding (TB) calculations, is presented in Ref. [19] for GaAs and in Refs. [14,20] for Si. It consists of taking the usual screened Coulomb potential ∼ 1/ri to correct the on-site energies at atomic positions ri 6= 0. For the impurity site the Coulomb correction would diverge, which is incompatible with TB, so a finite correction U0 is chosen to adjust the correct binding energies of each donor. Realistic descriptions of the donor Coulomb potential should be consistent with the expected limits: Vimp (r → ∞) = −e2 /4πǫSi r = VH (r) and Vimp (r → 0) = −e2 /4πǫ0 r. A convenient interpolation [21] is adopted here, see Fig. 3, namely     e2 1 1 1 Vimp (r) = − e−r/rcc + − 4πr ǫSi ǫ0 ǫSi = VH (r) + Vcc (r) (7) This involves a single donor-dependent parameter, rcc , giving the range around the donor site where the central cell correction is effective, as opposed to the three quantities required for the description given in Eq. (6). Values for rcc appropriate for P, As, Sb and Bi donors in Si are given in Table I. Note that while Bohr radii are of the order of a few nm, rcc is a factor of 10 smaller. Therefore the deviation of the plain Coulomb potential from the central-cell corrected one is restricted to a small region around the donor, as illustrated in Fig. 2(b). IV. NEUTRAL DONOR – D0 CENTER The coupling matrix elements calculation may be simplified by assuming spherically symmetric envelopes. It is possible to conciliate the simpler hydrogenic description 0 0 Donor Eexp (meV) [22] Ecalc (meV) rcc (nm) acc (nm) 0 -45.58 (A1 ) Eexp 0 P -33.90 (T2 ) -33.12 0.115 1.106 -32.60 (E) -33.12 As 0 Sb0 Bi 0 -53.77 (A1 ) -32.65 (T2 ) -31.40 (E) 0 Eexp -32.15 -32.15 0.126 0.815 -42.71 (A1 ) -33.00 (T2 ) -30.60 (E) 0 Eexp -32.04 -32.04 0.108 1.241 -71.00 (A1 ) -31.92 (T2 ) N/A (E) 0 Eexp N/A N/A 0.144 0.58 0 TABLE I: Experimentally measured [22] values Eexp of the 1s manifold energy states for the neutral isolated donor are com0 pared to the calculated values Ecalc obtained from the central cell corrected potential Eq. (7). The energy of the central cell corrected state A1 is fitted to the corresponding experimental value by choosing an appropriate rcc . The E and T2 levels are assumed degenerate. The calculation of the levels E and T2 involves Ẽ, the center of gravity of the sixfold 1s manifold, which is taken from experiment (see text), but is not available for Bi. Note that the experimental and calculated values for the T2 and E energies do not depend strongly on the donor species and are very similar to the single valley result EKL , in contrast to the value of the ground state A1 energy. As far as we know, the Bi E level is not available from experiment, but the T2 level is consistent with these observations. For each donor, the effective Bohr radius corrected by the central cell acc is also given. in Eq. (2) with the ground state energy EKL obtained from Eq. (4) by choosing a value of the effective mass m∗Si = 0.3 me . Note that m∗Si is very close to the geometric mean (mL × m2T )1/3 = 0.32. This approximation results, for the uncoupled valleys, into a sixfold degenerate ground state with ground state energy EKL . Now we take into account the multivalley structure of the Si conduction band and the central cell correction described in Eq. (7). The only parameter we are left with is rcc , the central cell potential characteristic length, chosen such that the ground state A1 experimental energy is reproduced for each of the donor species, as given in Table I. The T2 and E states are much less sensitive to the effect of Vcc since both have nodes at the impurity site. This is confirmed by their approximate degeneracy and energy very close to EKL (see Table I). We take the approximation of exact fivefold degeneracy of these levels, obtained assuming ∆k = ∆⊥ = ∆, so that E(A1 ) = EKL − δ + 5∆ and E(T2 ,E) = EKL − δ + ∆. The mean value of the six levels is taken from experiments and identified with 6 Ẽ = EKL − δ which leads to δ = EKL − Ẽ. It then follows that ∆ =[EKL −E(A1 )]/5 allowing to determine the excited states energies. Combining spherically symmetric envelopes with the Bloch functions obtained from the expansion coefficients given in Ref. [23] yield the electronic densities for the A1 and Tx2 bound donor states of P in Si shown in Fig. 2(c). Note that the charge distribution of the Tx2 state has a preferential alignment along the x direction and an exact node at the origin (idem to E states). Due to this node, any contact interaction (interaction potentials that are only finite at the immediate vicinity of the donor impurity) will have a reduced effect on the T2 and E states – this justifies why these states are not strongly affected by the central cell correction as well as indicates that low hyperfine coupling is to be expected in these states. For this reason, the T2 and E states energies are very similar to EKL . The results are summarized in Table I. V. TWO ELECTRONS IN ONE DONOR – D− -CENTER In analogy with the hydrogen atom, singly ionized negative (H− ) or positive (H+ ) states are also of practical and theoretical interest. The negative ion H− binding energy is defined as the energy required to remove one H− electron from the ion EB = EH 0 − EH − . For hydro− H gen EB = 0.055EH . The analogous D− center is specially important in Si nanoelectronics because its binding energy is very small and can be strongly affected by gates [24]. The presence of this bound state has been identified in transport experiments in single-atom transistors [25]. We may also regard D− centers as He-like systems [26], the main technical difficulty to obtain the spectrum consists in the evaluation of electron-electron repulsion terms. Variational wavefunctions with a large number of parameters, as used by Hylleraas [27] and others, are H− able to reproduce the experimental value of EB . In general, the ground state wavefunction is composed of an inner orbital with a Bohr radius similar to the electron bound to a neutral donor, and an outer orbital with a much larger radius, due to the screening of the Coulomb potential produced by the occupied inner orbital. Within a single valley and isotropic mass approach for negatively ionized centers in Si, we can just rescale the D− exact results for H − : EB = 0.055E ∗ ∼ 1.7 meV, which gives a good estimate for Phosphorus and Arsenic cenP− As− ters: EB = 1.7 meV, EB = 2.05 meV as measured in photoconductivity experiments [28]. The inclusion of the mass anisotropy and the valley degeneracy in Si would increase the binding energy of the D− centers even if valley-orbit interactions are not considered [29]. If no central cell corrections are included, each electron remains in a different valley. The lowest energy configuration is attained when the two electrons are located in perpendicular valleys (as mk and m⊥ are rotated in perpendicular valleys). Mass anisotropy then implies that the envelopes corresponding to the two orbitals do not overlap as much as they would for isotropic wavefunctions reducing the electron-electron repulsive interaction. Valley-orbit interactions and central-cell corrections have been treated in the literature for the D− problem under different approximations [18,30,31]. Larsen [30] assumed that only the inner orbital electron is subject to valley-orbit coupling, behaving like a neutral donor electron. The isotropic condition ∆k = ∆⊥ was also used. The effect of the valley-orbit coupling on the inner electron is to spread it in a symmetric (A1 ) combination of valleys. If both electrons have the same valley composition, the binding energy decreases with respect to them being in perpendicular valleys. Oliveira and Falicov [18] went a step further by including the full 1s multiplet in the description of the bound electrons and a constant valley-orbit coupling different from zero only in a centralcell region. More recent calculations [31] included the valley degeneracy but not the valley-orbit splitting. Here we treat the D− donors on the same level of approximation as the neutral donors, hence neglecting the mass anisotropy. We take the simplest possible form for the spin singlet ground state trial function (normalized): ψD − = 1 π p 2a31 a32  e−r1 /a1 e−r2 /a2 + e−r2 /a1 e−r1 /a2  , (8) and consider the central-cell corrected potential Vimp (r), as defined in Eq. (7). The correction term Vcc is completely determined by the value of rcc chosen from the corresponding neutral donor data (see Table I). Both a1 and a2 are calculated variationally to minimize the expectation value of the energy for the two-electrons hamiltonian: HD− = K1 + K2 + Vimp (r1 ) + Vimp (r2 ) + e2 /ǫSi r12 where Ki is the kinetic energy of electron i and Vimp (ri ) is the central cell corrected potential. The last term is the electron-electron repulsion. The energies are calculated following the prescription for He atoms given in Ref. [26]. If central cell corrections were ignored, the variational wavefunction in Eq. (8) would lead to a binding energy D− EB = 0.027E ∗ = 0.8 meV, see Ref. [24]. The inclusion of the central cell correction allows for a significant improvement of the binding energy, see Table II. The Table also shows the value of the charging energy U = EH − − 2EH 0 which can be measured in transport spectroscopy experiments in single atom transistors [24]. Note the very good agreement between calculated and experimental (when available) values in Table II, indicating the quantitative validity and transferability of the central cell correction given here and validating our predictions for Sb and Bi. 7 − − Eexp EB(calc) EB(exp) Ucalc Uexp a1 a2 Donor Ecalc − P -47.10 -47.29 1.526 1.7 43.86 43.00 1.041 4.851 As− -56.44 -55.81 2.69 2.05 51.06 51.71 0.737 3.516 Sb− -44.00 N/A 1.29 N/A 41.44 N/A 1.177 5.349 Bi− 7.55 N/A 63.45 N/A 0.48 2.16 -78.55 N/A TABLE II: Calculated and experimental [28] energies and variational parameters for the negatively charged donors. All energies are given in meV and lengths in nm. EB refers to the first ionization energy, namely, the energy required to ionize the less bound electron EB = E 0 − E − . The charging energy U is defined as U = E − − 2E 0 . The parameters a1 and a2 are the Bohr radii obtained variationally. The radius a1 corresponds to the inner orbital, and is very similar to acc in Table I. The parameter a2 is the radius of the outer orbital which is much less bound than the inner one. To the best of our knowledge, the experimental values for Sb and Bi are not known. VI. ONE ELECTRON AND TWO DONORS – D+ 2 CENTER The problem of two donors sharing a single electron is analogous to the problem of an ionized H2+ molecule in vacuum. If it was identical it would be possible to solve this problem exactly through the transformation into spheroidal coordinates [32]. However, as discussed for single donors, the effective mass anisotropy and the valley orbit coupling in silicon complicate the spectrum. Many accounts of these effects are available in the literature. Early approaches employed an altered form of the spheroidal coordinates to accomodate mass anisotropy to some extent [33]. We steer the analogy with the H2+ molecule along another direction, namely writing the eigenstates of this problem as bonding and anti-bonding molecular orbitals based on the single donor (atomic) problem, i.e., taking molecular wavefunctions as linear combination of atomic orbitals (LCAO). The two donors, referred to as A and B, located at positions rA and rB are separated by a vector R = rB − rA . An obvious difference between the substitutional donor pair and the molecular problems is the meaning and range of possible values of R. For donors, R is fixed at the sample fabrication and doping stage and coincides with a lattice vector, while in the free molecule situation, R = |R| is given by the minimization of the molecule energy. The effective mass hamiltonian of a donor pair reads HDD (r) = − ~2 ∇2 + Vimp (rA ) + Vimp (rB ). 2m∗ (9) The matrix elements of this hamiltonian are calculated in the LCAO basis, i.e., in the basis set defined by the single donor wavefunctions {A1 ,T2 , E} given in Eq. (5), centered at rA and rB . We expect this basis to give a fair account of the lowest electronic states. The hamiltonian is hence a 12×12 matrix, which for convenience we break into four 6 × 6 blocks, H= " HAA HAB HBA HBB # (10) Blocks HAA and HBB are not strictly diagonal since the problem of two donors has lower symmetry than a single donor. But a direct calculation shows that the offdiagonal terms in these 2 blocks are vanishingly small compared to other non-zero terms for interdonor distances larger than the lattice constant aSi . We therefore take the diagonal blocks to be diagonal, which shows that the donor eigenstates basis set is more convenient than the {φµ } for the pair treatment. These blocks, therefore, represent the on-site energies while the tunnel coupling terms contribute to the blocks HAB and HBA . The diagonal blocks [HAA ] = [HBB ] considering a basis ordered as in the sequence in Eq. (5):  εonsite 0 0 0 0 0 CC  0 εonsite 0 0 0 0    sv   onsite 0 εsv 0 0 0   0 [Hii ] =  . onsite  0  0 0 ε 0 0 sv   onsite  0 0 0 0 εsv 0  0 0 0 0 0 εonsite sv (11) The onsite energy for the A1 state εonsite is the ground CC state energy for the neutral donor including the central cell contribution from the A donor, corrected by ′ a long range classical term Vimp = hΨA1 (rA )|VH (rB ) + VCC (rB )|ΨA1 (rA )i to take into account the B donor potential. For T2 and E, εonsite is the single valley EKL sv corrected by VH′ = hΨT /E (rA )|VH (rB )|ΨT /E (rA )i. Following the argument in Sec. IV, the central-cell correction on T2 and E is neglected because these states have nodes at all Si sites.  The off-diagonal blocks are related by the hermiticity † condition HAB = HBA . Each term is a summation over integrals of the type 8 hφµ (rA )|HDD |φν (rB )i = Z F (rA )e−ikµ .rA u∗µ (rA )HDD F (rB )eikν .rB uν (rB )d3 r. These matrix elements can be straightforwardly calculated using the plane wave expansion of the Bloch functions (see Ref. [23]), and no further approximation is needed. Nevertheless, it is useful to simplify these integrals through some well tested approximations: (i) neglecting the matrix elements for µ 6= ν , since rapidly oscillatory integrands are involved; (ii) taking u∗µ uµ ≈ 1, as suggested in Ref. [23]. Under these assumptions we obtain B,j ikµ ·R hφ̃A,i tsv (R), µ |HDD (r)|φ̃ν i = δµ,ν e where tsv is the single valley tunnel coupling Z tsv (R) = F (rA )HDD F (rB )d3 r. (13) 𝜀𝑢⊥ 𝜀𝑢∥ 𝜏𝑢⊥1, 𝜏𝑢⊥2 𝜏𝑢∥ 𝜀𝑔∥ 𝜀𝑔⊥ E T (b) Molecular Regime 𝑹 ≈ 𝒂𝒄𝒄 E T E T E 𝐴1 𝐴1 𝐴1 T 𝜏𝑔∥ 𝜏𝑔⊥1, 𝜏𝑔⊥2 𝐴1 𝛼𝑢 𝛼𝑔 (c) (14) The effect of the hopping blocks HAB is analogous to the H2 molecule – at distances much larger than acc , this block is negligible and the states centered around sites A and B are degenerate, while for distances comparable to the one-atom wavefunction extension, gerade and ungerade combinations of the localized orbitals form, leading to split energies, see Fig. 4(a) and (b). At smaller interdonor distances, the gerade - ungerade splitting energy becomes comparable to the A1 -T valley-orbit splitting and a level inversion of the first excited state is obtained. This result is discussed theoretically in Ref. [34] and confirmed experimentally in Ref. [35]. Unlike the H2 molecule problem, it is possible for the antisymmetric state (referred to as antibonding in the context of H2 in vacuum) to be the ground state here, since the hoppings are not necessarily real negative numbers due to the oscillatory phase exp(ikµ · R). Still, we refer to the lower molecular orbital as gerade and the higher as ungerade. The same oscillatory phase may also lift the degeneracies of the T2 and E states. For instance, if the pair alignment is along the x direction, the states T2y and T2z are still equivalent, while the state T2x will have a symmetricantisymmetric splitting that oscillates as a function of R. This effect is further enhanced in the presence of the effective mass anisotropy, as seen in Ref. [34]. The effective mass anisotropy also impacts the interdonor distance at which the valley inversion of the first excited state occurs. It is known from the H2 problem that the molecular orbital approximation gives accurate results only if the variational wavefunction radius is taken to minimize the expectation value of the complete hamiltonian containing the two protons. We do the same here, obtaining a variational radius aD+ (R), which converges to the single 2 atom orbit at large distances aD+ (R → ∞) = acc . 2 W𝐞𝐚𝐤 𝐂𝐨𝐮𝐩𝐥𝐢𝐧𝐠 𝐑𝐞𝐠𝐢𝐦𝐞 𝑹 ≫ 𝒂𝒄𝒄 (a) (12) FIG. 4: Schematic behavior of the molecular orbitals spectra from the 1s manifold with different valley compositions in the (a) weak coupling regime for a single electron and (b) molecular (strong coupling) regime for a single electron. (c) Variational radii calculated to minimize the energy of the donor molecule oriented along the h100i, h110i, and h111i directions (squares, triangles and circles respectively). The variational radii for the ionized aD+ (empty symbols) and neutral aD0 2 2 (solid symbols) molecule are shown. The resulting energies for the ground [E0 (1e− )] and the first excited [E1 (1e− )] states are plotted with symbols in Fig. 5(a). The squares, triangles and circles correspond to the three different molecule orientations considered h100i, h110i and h111i. The oscillations with R are due to the intervalley interference which produce the incommensurate oscillations on the single donor wavefunctions, illustrated in Fig. 2(c). This leads to subtle oscillations in the total energy, but significant oscillations in the energy difference ∆0−1 = E1 (1e− ) − E0 (1e− ) [Fig. 5(b)]. VII. TWO ELECTRONS IN TWO DONORS – D02 CENTER The quantum mechanical solution to the problem of two impurities with two interacting electrons is significantly less studied, even though it is the most commonly found configuration since it is neutral. The electronelectron correlations are hard to describe and add significant complexity to the picture described in the previous section. The most systematic way to calculate with high 9 (a) (b) (c) (d) FIG. 5: Ground and first excited state energies for one and two electrons. (a) The ground state is the bonding αg orbital and the first excited is the antibonding αu in the weak coupling regime (R ≫ acc ) or the bonding τg in the strongly interacting regime (R ≪ acc ), see Fig. 4. The dashed lines are the LCAO results using single-valley central cell corrected 1s wavefunctions, and taking the Hückel approximation for the off-diagonal integrals. (b) The energy separation between ground and first excited states ∆0−1 , which coincides with the tunnel coupling in the weak coupling regime. Squares, triangles and circles correspond to the h100i, h110i and h111i molecule orientations respectively. (c) Two electron energies. The ground state is a spin singlet composed by two electrons in the bonding αg orbital. The first excited state is a spin triplet, with one electron in αg and the second electron in either αu or τg (in the weakly or strongly interacting regimes, respectively). The dashed lines show the Hubbard approximation. (d) The singlet-triplet separation ∆ST , which coincides with the spin-spin Heisenberg exchange coupling in the case of two electrons. accuracy the two electron energies is to build a full Configuration Interaction (CI) wavefunction from the molecular orbitals discussed for the D+ 2 problem. Now the suitability of truncating our basis set at the 1s manifold becomes clear. This manifold contains 6 atomic orbitals (one for each valley configuration), leading to 12 molecular orbitals and 144 two-electron Slater determinants. This way, the CI matrix for this problem would contain 20736 elements, each consisting of a few Coulomb integrals – six-dimensional improper numerical integrals with an integrand diverging at all points where the electron-electron distance is null r12 = 0. Even a small increment of the atomic basis set increases the numerical demand of this problem significantly. Still the problem stated above is too hard to be solved without any algebraic maneuvers and approximations. The first step is to identify elements that are null by symmetry. In the case of two electrons, one can always write down the wavefunction as a product of orbital and spinorial parts (with three or more particles these properties are necessarily entangled). Since the total spin operator and its projection both commute with the hamiltonian, it is wise to rewrite the molecular orbital basis set to be composed of eigenstates of these operators (singlets and triplets). This constitutes the so-called SACI (Spin Adapted Configuration Interaction) [36], which leads to a CI matrix composed by one 78 × 78 singlet block and one 66 × 66 triplet block, with null off-diagonal blocks. This cuts the computational effort by roughly a factor of two. These dense singlet and triplet matrices are still too challenging if all the Coulomb integrals are calculated numerically. A very robust quantum chemistry approximation to overcome this problem is to expand the Slatertype orbitals into a series of N gaussians with given radii. This approach converges quickly with N , and for the calculations performed here for interdonor distances up to R = 20 nm, the energies were well converged with N = 3. For larger distances the number of gaussians increases due to the ill described tail decay [37]. The results for the numerical calculation of the singlet and triplet energies are presented in Fig. 5(c) while the Bohr radii aD20 resulting from the variational minimization are in Fig. 4(c). The singlet triplet separation ∆ST , in Fig. 5(d) reveals the same oscillations due to valley interference [15,16] as the single electron excitation energy ∆0−1 . The CI method is considered the standard model of modern quantum chemistry. In principle, arbitrarily accurate results may be obtained if a large enough single particle basis set is adopted. In practice, though, a truncated basis set could lead to largely incorrect results – specially if the chosen basis set is inadequate. For instance, in Sec. V we adopted a simple two electron wavefunction, Eq. (8), which describes with a single orbital the energy of the D− center with great accuracy. If instead a regular CI basis set is constructed from the multiorbital LCAO method, the variational energy obtained is much higher and the result is qualitatively wrong – an unbound D− state is obtained, in contrast to the experimentally measured bound D− . VIII. A. DISCUSSION Simplified Models and effective Hamiltonians In certain studies, the description of donor states in Si based on microscopic model hamiltonians as presented in the previous sections may not be necessary or useful. Instead, simple models may be more adequate. In this 10 section we investigate possibilities to map the problem into simplified models. 1. D+ 2 : Single Orbital LCAO and the Hückel Approximation. Inclusion of higher excited orbitals improves the accuracy of the calculations, but in the simplest picture a single 1s orbital is considered, so that only two molecularorbitals are obtained. Another simplifying assumption is to disregard the effect of valley physics. This approximation disregards the oscillations in the tunnel coupling. Instead, the tunnel coupling becomes a real monotonic function of the interdonor distance, and the molecular states are the symmetric and antisymmetric combinations of the atomic orbitals
√ B |αg i = |AA 1 i + |A1 i / 2,
√ B |αu i = |AA 1 i − |A1 i / 2. (15) (16) The energies of these two states are E = Eonsite ± |t|, where Eonsite is the onsite energy of the atomic orbital and t is the tunnel coupling, or hopping parameter. The onsite energy is readily obtained. We may break down the hamiltonian as H = HA + VB , so that HA is A identical to the single donor case, and hAA 1 |HA |A1 i = E0 is simply the single electron in a single donor energy, known with great accuracy from experiments, hence we A take E0 = Eexp . The second term hAA 1 |VB |A1 i could in principle be calculated exactly, but since our wavefunction description is not particularly accurate, specially for R ≈ acc , we take a simplified point charge interacA 2 2 tion, such that hAA 1 |VB |A1 i = −e /4πǫSi R. e /4πǫSi ≈ 133.5 meV nm. The tunnel coupling could also be explicitly calculated, but a simpler result is obtained using the Hückel approximation. The argument behind this approximation is that the atomic orbital is an exact eigenstate of the single donor part HA , and the second center potential VB may be treated approximately as a percentual correction over the HA . In our case, we go further in the approximation A B A and take hAB 1 |VB |A1 i = 0, leading to t = hA1 |H|A1 i ≈ A B hA1 |HA |A1 i = E0 S(R, acc ), where the overlap function A −R/acc (1 + R/acc + R2 /3a2cc ). is S(R, acc ) = hAB 1 |A1 i = e A comparison between these approximations and the results for the complete numerical calculation, shown in Fig. 5(a) and (b), reveals that this approximation is well suited for estimating the order of magnitude of the tunnel coupling (estimated as the separation between the ground and first excited states ∆01 in our complete multiorbital LCAO method of Sec. VI), as well as the ground and first excited energies to an accuracy of approximately 10-15%, as seen in Fig. 5. The spread of the tunnel coupling in the numerical results is a product of valley interference, which is disregarded here. 2. D02 : Hubbard Model. The energy for two electrons may be calculated with similar arguments. The electron-electron repulsion supresses the probability of double occupation of the same donor. We take into account the correction to the singlet ground state due to the virtual occupation of this excited state, but discard any contributions from the triplet with both electrons at the same site since it has a much higher onsite energy. This is one of the ingredients of the so-called Hubbard approximation. This way, the triplet energy is the easiest to calculate, consisting only of the binding energies of the two electrons to each respective donor and the classical attraction/repulsion across donors/electrons. Hence, Etrip (2e− ) = ε(1,1) = 2E0 − e2 /4πǫSi R. The singlets are obtained diagonalizing the hamiltonian written in the basis {|s√A (r1 )sA (r2 )i, (|sA (r1 )sB (r2 )i + |sB (r1 )sA (r2 )i)/ 2, |sB (r1 )sB (r2 )i}, which reads   √ ε(2,0) + U 2t 0 √ √   [H]singlets =  2t ε(1,1) 2t  . (17) √ 0 2t ε(0,2) + U We take the coupling between the (1,1) and (2,0)/(0,2) states to be determined by the single particle hopping t, which is calculated within the Hückel approximation discussed previously. This approximation is not necessarily accurate, since the size of the wavefunction of (2,0)/(0,2) is most likely comparable to the size of the doubly occupied single donor ion D− . The most important parameter for the Hubbard model is the onsite charging energy U , which we take to be the experimental charging energy of the D− single dopant discussed in Section V (with exception of Sb and Bi, for which we adopt the results calculated in Sec. V). By symmetry we have ε(2,0) = ε(0,2) in the absence of detuning fields. The electrostatic arrangement of the (1,1) and the (2,0)/(0,2) charge configurations is significantly different, and therefore the onsite energy of these two states is not the same. On the other hand, the large same-site charging energy U leads to strong Coulomb blockade, and the occupation of the (2,0) and (0,2) states is only virtual. We argue that the error imposed by the approximation ε(1,1) = ε(2,0) has a very small impact in the ground state energy [composed almost exclusively by the (1,1) configuration] as well as the singlet-triplet separation ∆ST . These assumptions are confirmed a posteriori, see the comparison between the numerical and the approximated values in Fig. 5(c) and (d). B. Implications for experiments These results indicate possible directions for the design and characterization of one and two donor quantum devices. They may also bring new insights in specific 11 E(1e− ) = E0 − 133.5 R E(2e− ) = 2E0 − − t(R, acc ) 133.5 R − J(R, acc ) t(R, acc ) = E0 e−R/acc (1 + cc ) J(R, acc ) = 2 t(R,a U R acc + R2 ) 3a2 cc 2 Donor E0 (meV) acc (nm) U (meV) P -45.58 1.11 43.0 As -53.77 0.82 51.7 Sb -42.77 1.24 41.4 Bi -71.00 0.58 63.45 TABLE III: Explicit expressions (see Sec. VIII A) for the energies of the ionized donor molecule E(1e− ) and the neutral donor molecule E(2e− ) as a function of the distance R between the donors and the Bohr radius calculated for the neutral donors acc . Input parameters for the numerical estimates are taken from Tables I and II. The quantity 133.5 is given in units of meV nm (see text). aspects of larger systems such as donor clusters [38], impurity chains [39] and δ-doped systems [40]. Recent experimental results confirm the adequacy and validity of the present theory [35,41]. 1. tion of the donor interaction, since it is larger than that of a single donor D− state. This regime may be adopted for few dopant single electron transistors, leading to an enhanced charging energy and therefore to a more stable operation at room temperature. Transport through coupled donors. One of the most reliable techniques for probing single or few impurities is quantum transport spectroscopy [25, 41–45]. From these measurements the energy differences between occupation numbers (the chemical potential) are accessed directly. From the expressions in Table III, one can see that the energy of one and two electron states contain the same long distance electrostatic term. Therefore, the difference between these two energies – i.e. the first ionization energy – indicates there is quantum interaction between the two donors. If the tunnel coupling between the donors is not comparable to the binding energy, then the first ionization energy is essentially the same one would obtain for a single donor, i.e. E(2e− ) − E(1e− ) ≈ E0 . We refer to this as the weak coupling regime [45]. This is the most convenient regime for electrical control of the charge distribution and spin-spin interactions necessary for quantum computation with electron spins. The opposite regime consists of small distances – referred to as the molecular regime [41] – in which an electronic cloud is shared between the sites A and B with no classically forbidden region between the two donors. The first ionization energy departs strongly from the single donor binding energy, in such a way that the identification of donor pairs in this regime is easier. The charging energy also provides an important independent confirma- 2. Design of quantum devices. Silicon quantum devices take advantage of the quantum behavior of electrons in this semiconductor for specific tasks. Targetting parameters determined by quantum mechanics is a very challenging task. Purely quantum mechanical quantities, such as tunnel coupling and Heisenberg exchange coupling, depend explicitly on the wavefunction overlap and are therefore extremely sensitive to the asymptotic decay tails of the orbitals. Our non-perturbative approach provides potentials and wavefunctions tailored for each chemical species of the group V substitutional donors. In Fig. 6 we show (a) the tunnel coupling within single orbital LCAO theory and (b) the Heisenberg exchange coupling calculated within the Hubbard model. Both differ significantly from the results obtained using the KL wavefunction. The KL wavefunction overestimates the real tunnel and exchange coupling since it does not account for the effect of the valley-orbit coupling in shrinking the envelope function. The simplified expressions in Sec. VIII A, on the other hand, involve some approximations that overestimate and some that underestimate the tunnel and exchange couplings. For instance, taking a spherical effective mass of 0.3me leads to a smaller wavefunction overlap compared to the transverse effective mass. On the other hand, disregarding valley interference and the wavefunction size dependence on the interdonor distance R leads to an increase in the overlap. We expect the net effect of these approximations to have a lesser negative impact than the approximations in the KL theory, besides providing a path to differentiate and compare impurity chemical species. We highlight particularly important values of the tunnel and the exchange coupling. A modest impurity concentration generates impurity bands that are separated by the Mott gap due to the onsite electron-electron repulsion U . When the tunnel coupling – which determines the impurity band width – becomes comparable to U , the system becomes a metal at half-filling. This is the Mott metal-insulator transition. The detailed interdonor distance at which this transition occurs for each impurity is hard to predict – strong electronic screening is expected at the phase transition. We highlight the distance at which the unscreened repulsion U matches the tunnel coupling t. Another characteristic interdonor distance is that at which the exchange coupling becomes comparable to the nuclear hyperfine coupling. At these distances, the interaction between neighboring nuclei mediated through its electrons becomes optimally feasible [46]. This distance is very different amongst different species, as observed in 12 t≅U 70meV 10meV 60meV 50meV Tunneling t 1meV 40meV 100µeV 0.0 0.5 10µeV Sb Bi 1µeV Heisenberg Exchange J 1.0 P As J≅A 1meV 121 Sb 123 Sb 1µeV Sb 1neV Bi 4 6 8 P As 10 12 14 16 18 20 R (nm) FIG. 6: Quantum interactions for different impurity species. The tunnel rate is plotted according to the single orbital LCAO simplified expression in Table III. The data points for which t = U are marked, adopting the values of U listed in Table III. The exchange coupling is obtained from the Hubbard model, with the expression also shown in Table III. The marked values are the data points for which J = A, where the hyperfine coupling A is taken from Ref. [48]. The values of the nuclear spin I and the hyperfine coupling for P, As, Sb121 , Sb123 and Bi are, respectively, I =1/2, 3/2, 5/2, 7/2 and 9/2; and A =117.5 MHz, 198.3 MHz, 186.8 MHz, 101.5 MHz and 1.4754 GHz. Fig. 6. Moreover, this value is markedly different from the hydrogenic estimate adopted originally by Kane [47], which would set this distance at R = 23 nm independently of the chemical species. On the other hand, the possibility of shuttling a single electron across the donor pair through a time dependent electric field (in analogy with the quantum operation of electrostatic quantum dots) is most efficient at R > 20nm, for which the tunnel rate is in the range t/h = 10 kHz–10 MHz. IX. SUMMARY AND CONCLUSIONS Current progress in single or few dopants characterization [3] requires higher degree of sophistication where not only the nature (donor or acceptor) but actually identifying the atomic species of each impurity should be accessible to experiment. It is thus desirable to develop simple and reliable theories underlying the physics ∗ 1 Electronic address: also@if.ufrj.br Letter to Peierls, 29 September 1931, Wolfgang Pauli Wissenschaftlicher Briefwechsel mit Bohr, Einstein, Heisen- of such systems. We have presented a comprehensive study of donors in Si, which includes a central cell corrected potential obtained from consistency requirements with experiment. A single species-dependent parameter rcc reproduces data for neutral P, As, Sb and Bi levels within the 1s sixfold lower energy manifold in Si. The parameter rcc characterizes the range of the correction potential, and suffices to differentiate each species and respective potential. We verify by experimental comparison (when available, otherwise by plausibility arguments) that the proposed central cell corrected potential for a given species in Si (D0 ) is transferrable to other contexts of this species in Si, such as different charge states (D− ) and/or other atomic 0 arrangements, e.g. donor pairs (D+ 2 and D2 ). All approximations involved in different levels of the study are exposed and the range of validity or limitations are explicit. More realistic studies exist in the literature, at the cost of treating only specific centers. For example, while we assume isotropic envelopes, a recent study [34] preserves the Si band anisotropy effects in an analysis restricted to D+ 2 (one-electron donor pair ions), with which our results are in fair agreement. Among the calculated properties, the singlet-triplet difference ∆ST [see Fig. 5(d)] in the neutral donor pair deserves special attention for donor-based qubits applications. This parameter may be identified with the exchange coupling between the electrons, which is the basic two-qubit entanglement mechanism in the first Sibased quantum processing proposal by Kane [47]. We obtain the spatial range of the wavefunctions, quantified by acc , significantly contracted relative to previous estimates based on the screened point-charge potential. One implication is that sizeable exchange coupling between electrons bound to a donor pair requires interdonor distances significantly smaller than previously expected and attempted. This is probably one of the reasons why exchange-coupled pairs were only reported very recently for As [41] and P [35] in Si. In both cases the models detailed here were instrumental to understand the reported results. Acknowledgements. The authors thank M. Friesen and M. Eriksson for fruitful discussions. AS, AB, and BK performed this work as part of the Brazilian National Institute for Science and Technology on Quantum Information and also acknowledge partial support from the Brazilian agencies FAPERJ, CNPq, CAPES. AS also acknowledges the William F. Vilas Trust for financial support. MJC acknowledges support from MINECO-Spain through grant FIS2012-33521. AS, BK and MJC acknowledge support from a bilateral CNPq (Brazil)- CSIC (Spain) grant. 2 berg u.a. Band II: 1930-1939, Springer, 1985, p. 94 Awschalom D D, Bassett L C, Dzurak A S, Hu E L, Petta J R 2013 Science 339 1174 13 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 Zwanenburg F A, Dzurak A S, Morello A, Simmons M Y, Hollenberg L C L, Klimeck G, Rogge S, Coppersmith S N, and Eriksson M A 2013 Review of Modern Physics 85 961. Morley G W, arXiv:1407.6250. Ramdas A K and Rodriguez S 1981 Rep. Prog. Phys. 44 1297. Fetterman H R, Larsen D M, Stillman G E, Tannenwald P E and Waldmant J 1971 Phys. Rev. Lett. 26 975 Salfi J, Mol J A, Rahman R, Klimeck G, Simmons M Y, Hollenberg L C L and Rogge S 2014 Nature Materials 13, 605 Luttinger J M and Kohn W 1955 Phys. Rev. 97 869. Kittel C and Mitchel A H 1954 Phys. Rev. 96 1488. For the single donor problem it is convenient to take the donor position as the coordinate origin r = 0. Kohn W 1957 ”Shallow impurity states in silicon and germanium.” Solid State Physics 5 257-320. Pantelides S T, Sah C T 1974 Phys. Rev. B 10, 621. Nara H, Morita A 1966 J. Phys. Soc. Jpn. 21, 1852. Martins A S, Boykin T B , Klimeck G, Koiller B 2005 Phys. Rev. B 72, 193204. Wellard L C, Hollenberg L C L 2005 Phys. Rev. B 72, 085202. Koiller B, Hu X and Das Sarma S 2002 Phys. Rev. B 66, 11520. Baena A, Saraiva A L, Koiller B, Calderón M J 2012 Phys. Rev. B 86, 035317. Oliveira L E, Falicov L.M 1986 Phys. Rev. B 33, 6990. Menchero J G, Capaz R B, Koiller B, and Chacham H 1999 Phys. Rev. B 59, 2722. Martins A S, Menchero J G, Capaz R B, Koiller B 2002 Phys. Rev. B 65, 24205. Wang J S -Y and Kittel C 1973 Phys. Rev. B 7, 713. Grimmeiss H G, Janzen E, and Larsson K 1982 Phys. Rev. B 25, 2627. Saraiva A L, Calderón M J, Capaz R B, Hu XD, Das Sarma S, Koiller 2011 Phys. Rev. B 84, 155320. Calderón M J, Verduijn J, Lansbergen G P, Tettamanzi G C, Rogge S, Koiller B 2010 Phys. Rev. B 82, 075317. Sellier H, Lansbergen G P, Caro J, Rogge S, Collaert N, Ferain I, Jurczak M, Biesemans S 2006 Phys. Rev. Lett. 97, 206805. Bethe H A, Salpeter E E 1957 Quantum mechanics of one and two electron atoms, Academic Press, New York. Hylleraas E A 1929 Z. Physik 54, 347. 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 Narita S, Shinbashi T, and Kobayashi M 1982 J. Phys. Soc. Jpn. 51, 2186. Kamimura H 1979 Journal of Non-crystalline solids 32, 187. Larsen D M, 1981 Phys. Rev. B 23, 5521. Inoue J, Nakamura J, Natori A 2008 Phys. Rev. B 77, 125213. Slater J C, 1963 ”Quantum Theory of Molecules and Solids”, Vol. 1. Miller A and Abrahams E 1960 Phys. Rev. B 120, 745. Klymenko M V and Remacle F 2014 J. Phys. Condens. Matter 26 065302. Dehollain J P, Muhonen J T, Tan K Y, Saraiva A, Jamieson D N, Dzurak A S and Morello A 2014, Phys. Rev. Lett. 112, 236801. Wang J and Hines C 2009, Journal of Physics: Conference Series 185, 012053. Saraiva A L, Calderón M J, Koiller B 2007 Phys. Rev. B 76, 233302. Weber B, Tan Y H M, Mahapatra S, Watson T F, Ryu H, Rahman R, Hollenberg L C L, Klimeck G and Simmons M Y 2014 Nature Nanotechnology 9, 430. Prati E, Hori H, Guagliardo F, Ferrari G and Shinada T 2012 Nature Nanotechnology 7, 443. Shamim S, Mahapatra S, Scappucci G, Klesse W M, Simmons M Y and Ghosh A 2014 Phys. Rev. Lett. 112, 236602. Gonzalez-Zalba M F, Saraiva A, Calderón M J, Heiss D, Koiller B, and Ferguson A J arXiv:1312.4589. Gonzalez-Zalba M F, Heiss D, and Ferguson A J 2012 New Journal of Physics 14, 023050. Lansbergen G P, Rahman R, Wellard C J, Woo I, Caro I J, Collaert N, Biesemans S, Klimeck G, Hollenberg L C L, Rogge S 2008 Nature Physics 4, 656. Pierre M, Wacquez R, Jehl X, Sanquer M, Vinet M, and Cueto O 2009 Nature Nanotechnology 5, 133. Roche B, Riwar R -P, Voisin B, Dupont-Ferrier E, Wacquez R, Vinet M, Sanquer M, Splettstoesser J and Jehl X 2013 Nature Communications 4, 1581 Kalra R, Laucht A, Hill C, Morello A 2014 Phys. Rev. X 4, 021044. Kane B E 1998 Nature 393, 133. Feher G 1959 Phys. Rev. 114, 1219.