Ab initio calculation of energy levels for phosphorus donors in silicon

The s manifold energy levels for phosphorus donors in silicon are important input parameters for the design and modeling of electronic devices on the nanoscale. In this paper we calculate these energy levels from first principles using density functional theory. The wavefunction of the donor electron’s ground state is found to have a form that is similar to an atomic s orbital, with an effective Bohr radius of 1.8 nm. The corresponding binding energy of this state is found to be 41 meV, which is in good agreement with the currently accepted value of 45.59 meV. We also calculate the energies of the excited 1s(T2) and 1s(E) states, finding them to be 32 and 31 meV respectively.

Phosphorus donors in silicon have long been important for electronic devices but are now seen as central to the development of silicon based quantum information processing [1][2][3][4][5] . The phosphorus donor electron has been shown to have long spin coherence times in the laboratory, which make these donors excellent candidates for spin qubits 2,4 . Moreover, during the last decade a technique of phosphorus δ doping, based on scanning tunneling microscope (STM) lithography, has led to a variety of new electronic devices in silicon 6 . This δ doping technique has been used to make a quantum dot of seven donors 7 and a transistor with a gate island that consists of only one phosphorus donor 3 . Another novel electronic device is the quantized electron pump of Tettamanzi et al. 8 , which demonstrates charge pumping of single electrons through a phosphorus donor. Finally, the wavefunction of the donor electron has even recently been imaged using an STM 9 . These images have been analysed using tight binding 10 and effective mass theory 11 ; whilst the latter provides a qualitative description, the tight binding method is precise enough to pinpoint the atomic position of a single phosphorus donor in the silicon lattice. Although semi-empirical approaches have successfully been used to model the properties of these donor devices, a full ab initio treatment of the electronic structure of these donors has to-date not been possible. Here we present such a treatment.
At low doping densities it is well known that the phosphorus donor electrons occupy the lowest energy conduction band of silicon. In bulk silicon this band is sixfold degenerate but the degeneracy is lifted by a valley splitting when silicon is doped 12,13 , resulting in three nondegenerate states. These states are, in order of increasing energy, a singlet [1s(A 1 )], a triplet [1s(T 2 )], and a doublet [1s(E)] 14 . Only the ground state [1s(A 1 )] is populated 14 at liquid helium temperatures (~4K), whereas at higher temperatures (≥30 K) the populations of the excited 1s(T 2 ) and 1s(E) states become observable due to thermal broadening 13,15 .
Over a decade ago, theoretical methods for describing point defects in semiconductors were separable into two categories: "methods for deep defects and methods for shallow defects: the former defect class is treated by ab initio methods, … while for the latter class approximate one-electron theories … are used" 16 . Traditionally, shallow defects in silicon like phosphorus donors could not be treated by ab initio methods because the wavefunctions of such defects are partially delocalized. Today, however, this statement does not hold true, as in the last ten years innovations in modern computing technologies have made much larger computational resources available to scientific research. Recently it has been shown that shallow defects are now within the reach of ab initio methods such as density functional theory (DFT) 17 .
In this paper we calculate the energies of the s manifold states [ s A 1 ( ) 1 , s T 1 ( ) 2 , and s E 1 ( )] of a phosphorus donor electron in silicon from first principles using DFT. We also compute the wavefunction of the donor electron's ground state [ s A 1 ( ) 1 ]. From this we estimate the effective Bohr radius of the electron by fitting to this wavefunction. We find DFT significantly underestimates the energies of the s manifold states. This is a known problem and, as will be discussed, we correct these energies using the ground state wavefunction, via the method described in ref. 17. In this way we are able to obtain ionisation energies for the donor electron that are in good agreement with the currently accepted values. To the best of our knowledge these results are the first ab initio confirmation of the s manifold energy levels for a phosphorus donor in silicon.
The Lyman spectrum for Group V donors in silicon was first measured by Aggarwal et al. in 1965 18 . These measurements do not give the binding energy of the donor electrons but rather the energy splitting between the ground and excited states; namely, the energy splitting between the s A 1 ( ) 1 and 3p ± states. The binding energy of the phosphorus donor electron that is reported in ref. 18 was computed by "adding the theoretically calculated binding energy of 2.90 meV for the 3p ± state 12 to the energy of the transition → ± s A p 1 ( ) 3 1 ". The binding energy of the phosphorus donor electron was thereby found to be 45.31 meV 18,19 .
In 1969, Faulkner used effective mass theory (EMT) to calculate the energy levels of the ground and excited states of a donor electron for Group V donor atoms in silicon 20 . For the phosphorus donor electron the binding energies of the 3p ± and s A 1 ( ) 1 states were found to be 3.12 meV and 31.27 meV, respectively 20 . The theoretically calculated binding energy of the excited 3p ± state is in good agreement with experiment, whereas the binding energy of the s A 1 ( ) 1 state is not 20 . Later, in 1981, using the theoretical correction of Faulkner 20 and a new experimental technique that produced narrower linewidths in the excitation spectra, Jagannath et al. 19 reported a binding energy of 45.59 meV for the phosphorus donor electron. Finally, the charging energy of a single phosphorus donor in the presence of two densely-doped, phosphorus leads was recently found to be 47 ± 3 meV in experiment. This value agrees, within experimental uncertainty, to the earlier optical measurements 19 .
More recently it has been demonstrated that EMT, with effective potentials calculated from ab initio methods, is capable of reproducing the accepted values for the binding energies of the s manifold states 21 . In addition, a model for a phosphorus donor in silicon that goes "beyond effective mass theory" has been introduced 22 . In ref. 22 the binding energy was used as a fitting parameter together with non-static screening effects in a model that provided an excellent account of the s manifold of states. This study shows that the binding energy is also an important quantity for theoretical modelling. The same fact is highlighted by ref. 23, where the hyperfine Stark effect is investigated using a truncated Coulomb potential to approximate the impurity potential of an ionized phosphorus donor 23 . The truncation of the Coulomb potential was found by adjusting a free parameter "to obtain the experimental ground state energy of 45.6 meV" 23 .
The binding energies of Group V donors in silicon have also been used as input parameters to modelling of the hyperfine Stark effect with EMT 24 . EMT has been shown to be capable of reproducing the wavefunction of a phosphorus donor electron that is predicted by tight binding theory 25 . The results in ref. 25 were benchmarked against the currently accepted value for the binding energy of a phosphorus donor electron in silicon. Knowledge of the binding energy, and specifically the valley splitting, was needed to choose the exact form of the central-cell corrections, i.e. a central cell with tetrahedral, rather than spherical, symmetry 25 .
The first large-scale atomic simulations performed on a Group V donor in silicon using DFT were those presented in ref. 17. In this study the electronic properties of an arsenic donor in silicon were calculated for systems that ranged in size from 512 to 10,648 atoms. DFT has also been used to simulate phosphorus donors in silicon, with systems ranging in size from 54 to 432 atoms 26 . However, as we will show, these latter system sizes are not large enough to isolate the phosphorus donor electron from its periodic images. The confinement of the donor electron is thereby increased, which artificially raises the binding energy of the electron. The binding energy of the phosphorus donor electron was therefore unable to be reported in ref. 26. The experimental hyperfine and superhyperfine interactions for Group V donors in silicon have been reproduced accurately using DFT and a Green's functions approach 16 . However, this study was also not able to report meaningful binding energies for the donor electrons because the long-range tail of the impurity potential was ignored 16 .
In this paper we present the results of electronic structure calculations performed on a single phosphorus donor in silicon with DFT. This approach has previously been benchmarked in a number of other studies [27][28][29][30][31][32][33] . For more information on this method and its benchmarking see Appendices A and B. We have employed the SIESTA package 34,35 to carry out calculations on cubic supercells that range in size up to 10,648 atoms.
We have calculated the wavefunction, ψ, of the donor electron's ground state for cubic supercells that range in size from 512 to 10,648 atoms. A two dimensional slice of the probability density, ψ 2 , for the largest supercell studied in this work is plotted in Fig. 1. This slice is computed by evaluating the wavefunction in the silicon (001) plane that contains the phosphorus donor. The maximum of the probability density in this slice has been normalized to one. In the (001) plane the majority of the probability density can be seen to be within ~0.5 nm of the donor site, which is located at the origin in Fig. 1. The wavefunction of the donor electron has a form that is similar to an atomic s orbital. The corresponding probability density can be seen to decay to approximately 2% of its maximum value at a distance of ±1.5 nm from the donor site in the [100] and [010] crystallographic directions.
The apparent hydrogenic character of the donor electron's wavefunction is compatible with an effective Bohr model of the electron. Phosphorus is a shallow defect in silicon so it is reasonable to treat the wavefunction of the Kohn-Sham eigenvalue, calculated within DFT, as an independent single particle state that can be modelled by a simple exponential function. The wavefunction of the donor electron can therefore be described by the envelope A is a normalisation constant and ⁎ a 0 is an effective Bohr radius. It is then possible to calculate the effective Bohr radius of the donor electron by fitting its wavefunction with this envelope function. However, it is first necessary to spherically average the wavefunction of the donor electron because F(r) is radially symmetric and ψ is not. Figure 2 shows the natural logarithm of the spherically averaged probability density for the phosphorus donor electron, ln(|ψ(r)| 2 ), plotted against radial distance from the donor site, r. The domain in this figure includes the core region of the phosphorus atom, which in our model is described by a Troullier-Martins pseudopotential 36 . A pseudopotential will deviate from a Coulombic potential in the core region. The envelope function is not applicable within the core region because a hydrogenic wavefunction is not a valid solution here. We have therefore fitted the wavefunction of the donor electron on the domain . R [ , 3 0] nm, where R is termed the model radius. The model radius must be chosen such that the effects of the core region on the wavefunction do not influence the accuracy of the exponential fit. Nor can the model radius be so large that the whole of the wavefunction's exponential decay is not captured by the fit. We have set the model radius equal to the atomic nearest neighbour distance, which has a value of 0.235 nm in silicon 37 . As can be seen from Fig. 2, this value for the model radius satisfies our two requirements.
In EMT, it is possible to derive two Bohr radii for the donor electron: one corresponds to the longitudinal effective mass, m , of bulk silicon and the other to the transverse mass, ⊥ m . The geometric average of these two radii is given by ⊥ a a 2/3 1/3 . By fitting F r ln( ( ) ) 2 to ψ r ln( ( ) ) 2 , we find the effective Bohr radius to be 1.8 nm. This value is in good agreement with 2.087 nm, which is the geometric average of the two effective Bohr radii reported in ref. 38. By reconsidering Fig. 1, we can see that the effective Bohr radius can be thought of as the radial distance within which the vast majority of the probability density corresponding to the donor electron is contained.
The s manifold energy levels for the phosphorus donor electron are shown in Fig. 3. These energies are plotted relative to the conduction band minimum of bulk silicon, calculated using a supercell of 10,648 atoms, which is set to energy zero in the figure. Figure 3 illustrates how the energies of the s A 1 ( ) 1 , s T 1 ( ) 2 , and s E 1 ( ) states increase  as the size of the supercell is increased. We suggest the energy levels for the smaller supercells are artifically lowered due to the electron's interaction with its periodic images, which increases the confinement of the donor electron as the size of the supercell is decreased. The energy levels are converged to within 1 meV for a supercell of 10,648 atoms. These results justify the use of such a large supercell for the calculation of these energies.
The binding energies of the s A 1 ( ) 1 , s T 1 ( ) 2 , and s E 1 ( ) states can be calculated for each supercell by taking the difference between the energy levels and the CBM of bulk silicon. In Fig. 3 the larger supercells significantly underestimate these energies. This discrepancy is due to the fact the Kohn-Sham eigenvalues are single particle energies and do not correspond exactly to true excitations of the system. Because the conduction band of bulk silicon is unoccupied, the energy of the CBM does not correspond exactly to a true excitation of the system when calculated by DFT. This is the well-known band-gap problem of DFT 39 . It is therefore incorrect to calculate the binding energy of the donor electron by taking the difference of this energy and the CBM of bulk silicon. We now need another way to calculate the binding energies of the s manifold energy levels. This is provided by the method described in ref. 17, where the binding energy is calculated from ψ.
The method of ref. 17 allows us to calculate the binding energies of the s A 1 ( ) 1 , s T 1 ( ) 2 , and s E 1 ( ) states directly from their wavefunctions and the impurity potential of the phosphorus donor. We begin with the screened impurity potential of the phosphorus donor 17 ; where V′(q) is the Fourier transform of the unscreened impurity potential. The dielectric screening is described by a nonlinear function 22,40,41 ; with A = 1.175, α = 0.7572, β = 0.3123, γ = 2.044, and the relative permittivity of silicon ε = .
(0) 11 4. The constants A, α, β, and γ were found by fitting the above function to the q dependent dielectric screening in silicon, which was calculated from the random phase approximation 41 . The calculation of the binding energy therefore relies on the assumption that the screening of the impurity potential is well-described in the random phase approximation, within which dynamical effects on this screening are ignored. The kinetic and potential energies of the donor electron can then be computed; 3 where ψ is the wavefunction calculated from DFT. Finally, we can calculate the binding energy of the donor electron: For more information, see Appendix C.  Table 1 presents the binding energies of the s manifold states calculated using experiment, EMT, and DFT. The binding energy of the s A 1 ( ) 1 state calculated using DFT (this work) is equal to 41 meV. This energy is in good agreement with the accepted value of 45.59 meV, which has been calculated from the combination of an experimental measurement 19 and a theoretical correction 20 . In addition, we find the binding energies of the excited s T 1 ( ) 2 and s E 1 ( ) states to be 32 meV and 31 meV, respectively. These values are in excellent agreement with the other values listed in Table 1, agreeing to within 2 meV. The binding energies of the two excited states appear to be in better agreement with the accepted values for these energies than the energy of the donor electron's ground state.
In summary, we have calculated the wavefunction of a phosphorus donor electron in silicon with DFT. This wavefunction is then used to compute the effective Bohr radius of the donor electron. We employ a hydrogenic model of this electron and thereby find its Bohr radius to be 1.8 nm. In addition, we compute the binding energy of the donor electron's ground state, which is found to be in good agreement with the currently accepted value. The energies of the excited s T 1 ( ) 2 and s E 1 ( ) states are found to be in excellent agreement with the accepted values. These results constitute the first ab initio calculation of the s manifold energy levels for a single phosphorus donor in silicon.
In the future this atomistic model could be expanded beyond the bulk case, to investigate the effect of a nearby surface on the electronic and structural properties of the phosphorus donor. Recent experiments on single donor atoms are commonly performed in the presence of such an interface 9,10 . The effect of a surface on the donor's properties cannot be studied using a one-electron theory, like EMT, because these theories implicitly assume the defect is surrounded by a homogeneous bulk. In addition our ab initio model could easily be extended to investigate many donor systems, which are also of current interest 42 .

Density functional method
The electronic structure calculations were performed with density functional theory (DFT) using the SIESTA package 34,35 . Table 2 lists each of the supercell sizes that have been studied by number of atoms and the dimensions of the supercells in real space. These calculations have been performed using periodic boundary conditions. We have employed the Perdew-Burke-Ernzerhof (PBE) exchange correlation (XC) functional in the generalised gradient approximation (GGA) 43 . Application of the GGA to phosphorus-doped silicon systems in the past has produced results that are in good agreement with experiment 44 . The total energies of each of the supercells were converged to within 0.1 meV using a planewave energy cutoff of 300 Ry and a Fermi-Dirac occupation function at a temperature of 0 K. Atomic potentials were described by norm-conserving Troullier-Martins pseudopotentials 36 .
We have variationally solved the Kohn-Sham equations using a basis set of localised atomic orbitals that was optimised for phosphorus-doped silicon using the simplex method 28 . The basis set was double-ζ polarised and was comprised of 13 radial functions. In ref. 30, localised single-ζ and double-ζ polarised bases, and a delocalised planewave basis were used to calculate the valley splitting for a phosphorus δ doped monolayer in silicon. Despite  the higher precision of the planewave basis, the double-ζ polarised basis was shown to "[retain] the physics of the planewave description" 30 . We relaxed the crystallographic structure of bulk silicon using this basis set and found the lattice constant to be 5.4575 Å. This value is in good agreement with the experimental value of 5.431 Å 45 . The overestimation of the lattice constant by approximately 0.5% is lower than the usual systematic deviation of the lattice constant that is expected from the PBE XC functional, which is a 0.7% deviation 46 .

Benchmarking of density functional method
To reduce the computational expense of performing these electronic structure calculations, we have used a k point grid that contains only a single k point: the Γ point, i.e. k = (0, 0, 0). For the supercell of 10,648 atoms, an increase in the size of the k point grid would result in these calculations being computationally impractical. When the number of k points is increased up to 8 × 8 × 8 for the supercell of 512 atoms, we find the eigenvalue of the s A 1 ( ) 1 state at the Γ point converges to a value that is approximately 5 meV greater than that of the Γ point calculation. The eigenvalues of the s T 1 ( ) 2 and s E 1 ( ) states converge to values that are approximately 1 meV greater than the result of their respective Γ point calculations. We expect these changes in the eigenvalues of the system to decrease as the size of the supercell is increased because the size of the corresponding Brillouin zone will decrease. Previous calculations of an arsenic donor in silicon with DFT have also been restricted to the Γ point 17 .
We geometrically optimised the ionic positions of supercells that ranged in size from 64 to 4096 atoms. We found the maximum displacement of a silicon atom was largest for the supercell of 64 atoms. When the size of the supercells is increased up to 4096 atoms the maximum displacement decreased to less than 0.02 Å. This displacement is equivalent to less than 0.5% of the lattice constant of bulk silicon. We therefore conclude it is unnecessary to relax the ionic positions of silicon atoms beyond their bulk values for supercells larger than 4096 atoms.
The conduction band minimum (CBM) of silicon is located at π ≈ . k a 0 85(2 / ), along each of the cardinal axes of reciprocal space, inside the face centred cubic Brillouin zone. Because silicon is an indirect bandgap semiconductor, the energy of the lowest conduction valley at the Γ point is not equal to the energy of the CBM. This is a result of the dispersion of the energy bands. We calculate the eigenvalues of the phosphorus donor electron at |k| = 0, not π ≈ . k a 0 85(2 / ), and therefore it is necessary to offset the computed energies of the s A 1 ( ) 1 , s T 1 ( ) 2 , and s E 1 ( ) states to find their value at π ≈ . k a 0 85(2 / ). The size of the Brillouin zone is decreased when the size of the supercell is increased. Decreasing the size of the Brillouin zone causes the bands, and therefore CBM, to be folded towards the centre of the zone, i.e. the Γ point, in a process known as band folding 30 . Consequently, the amount by which the energies of the s A 1 ( ) 1 , s T 1 ( ) 2 , and s E 1 ( ) states must be offset, to account for the parabolic dispersion of the band, is different for each supercell. The folding of the lowest conduction valley is plotted in Fig. 4a for supercells that range in size from 8 to 4096 atoms. The value of the offset for each supercell can be computed by taking the difference between the energy of the valley at the Γ point ( Γ E ) and the conduction-band minimum (CBM). These energies have been plotted for all supercells in Fig. 4b.
The value of − Γ E CBM decreases as the size of the supercell is increased. As shown in Fig. 4b, this relationship is not monotonic: the CBM is not always folded closer to the Γ point as the size of the Brillouin zone is decreased. Figure 4a shows the lowest conduction valley for bulk silicon only. If the dispersion of this band does not change significantly upon doping with phosphorus, then the difference − Γ E CBM can be used to correct the computed energies of the s A 1 ( ) 1 , s T 1 ( ) 2 , and s E 1 ( ) states. The positions of the conduction valleys on the k x axis, in Fig. 4, have been computed by folding the band structure of bulk silicon. The unfolded band structure was calculated using an eight atom simple cubic unit cell and a k point grid of 6 × 6 × 6. For the sake of clarity, we do not show the part of the bands that are reflected back into the Brillouin zone at the zone boundary. Neither do we show the conduction valleys of supercells with more than 4096 atoms in Fig. 4. The reflection of the bands at the zone boundary is a consequence of the fact that a solution in one Brillouin zone must be a solution in all Brillouin zones 30 .
The energy of the lowest conduction valley of bulk silicon at the Γ point ( Γ E ) is shown in Fig. 5 for every supercell studied in this work. As expected, the value of Γ E is different for each supercell. The conduction band minima plotted in Fig. 5 are for bulk silicon and have been calculated by substracting − The CBM for bulk silicon is not expected to change as the size of the supercell is increased. We therefore use the CBM of the supercell containing 10,648 atoms as a point of reference by setting it to energy zero in the figure.
We find the CBM for each of the supercells do not agree when the energies are corrected for band folding only. We also need to account for the differences in the valence band maximum (VBM) of each supercell. The conduction band minima are shifted by the difference between the VBM of each supercell and the VBM of the supercell containing 10,648 atoms. Once this is done, the CBM in Fig. 5 agree to within 4 meV. The remaining discrepancies in the conduction band minima could be caused by the differing k point grids that were used to calculate the quantity − Γ E CBM and the conduction band minima plotted in Fig. 5, or similar errors in the VBM itself. The VBM is not affected by band folding because it appears at the Γ point in the Brillouin zone.

Calculation of binding energies
In this section, we give the mathematical details of the calculation of the binding energy for the donor electron in full. This method was first proposed in ref. 17 for an arsenic donor in silicon.
The binding energy of the donor electron is given by  = + E T U (6) where T is the kinetic energy and U is the potential energy of the donor electron. In (6) the potential energy is defined as 3 where ψ is the wavefunction of the donor electron and V is the impurity potential for the phosphorus donor. The impurity potential can be written as P:Si 1e:Si where V P:Si is the electric potential for a phosphorus-doped silicon system and V 1e:Si is the electric potential for an electron-doped silicon system. By an electron-doped silicon system, we mean a bulk silicon system with one electron added. In contrast, for the phosphorus-doped system, one electron is added to the system by substituting a silicon atom with a phosphorus atom. These two electric potentials can be defined as  where V ee is the electron-electron contribution to the electric potential, V XC is the exchange-correlation contribution to the electric potential, and V eN is the electron-nuclear contribution to the electric potential. Substituting (9) and (10) into (8) and rearranging we have In the equations above, the impurity potential is screened by the electron-electron and exchange-correlation terms. Next, we set V ee and V XC to zero and thereby introduce a new quantity, the unscreened impurity potential V′. The unscreened impurity potential is given by the last term in (11):  Si,1e:Si where R 0 is the ionic position of the phosphorus donor atom, R i is the ionic position of silicon atom i, and V pp P and V pp Si are the pseudopotentials of phosphorus and silicon, respectively. Substituting (13) and (14) into (12), we obtain Si,P:Si S i,1e:Si is approximate because the norm-consering Troullier-Martins pseudopotentials are nonlocal. Let R 0 = (0, 0, 0), then pp pp P S i That is, the unscreened impurity potential is given by the difference in the pseudopotentials for phosphorus and silicon. We have used only the l = 0 component of the norm-conserving pseudopotentials for phosphorus and silicon when evaluating Eq. 15. This approximation is justified given the structure of the eigenfunction for the s A 1 ( ) 1 state (cf. Fig. 1 in the main text). Electron screening can now be reintroduced using the following description. We rewrite the screened impurity potential as ref. 17 where V′(q) is the Fourier transform of the unscreened impurity potential. The dielectric screening is described by a nonlinear function 40,41 ε α β with A = 1.175, α = 0.7572, β = 0.3123, γ = 2.044, and ε = .
(0) 11 4. The constants A, α, β, and γ were found by fitting the above function to the q dependent dielectric screening in silicon, which was calculated from the random phase approximation 41 . We can then use (16) to calculate the potential energy of the donor electron using (7). Finally, to calculate the kinetic energy of the donor electron, we use the virial theorem: 3 The binding energy of the donor electron can then be calculated from the kinetic and potential energies using (6). The binding energies of the donor electron's ground state, calculated using supercells of 216 to 10,648 atoms, are shown in Fig. 6. For the supercell of 10,648 atoms, the value of the binding energy is converged to within 1 meV. The accepted value for the binding energy of the donor electron's ground state, which is equal to 45.59 meV, is shown as a dashed line in Fig. 6. The supercell of 216 atoms overestimates the binding energy of the s A 1 ( ) 1 state in the figure but the binding energy decreases as the size of the supercell is increased. This energy is within 5 meV of the accepted value for a supercell of 10,648 atoms. Unfortunately there is no systematic way of calculating the uncertainty in this energy, but it is unlikely that the uncertainty is less than 5 meV.