Multivalley effective mass theory simulation of donors in silicon

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 theory results [Salfi et al., Nat. Mater. 13, 605 (2014)]. Here, we show that multivalley effective mass theory, applied properly, does achieve close agreement with tight-binding results 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. Although 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.

Figure 1

Multivalley effective mass calculations for a single phosphorus donor in silicon. (a) Sketch of the band structure of silicon and the resulting donor physics. The conduction band valleys are initially sixfold degenerate; valley-orbit coupling causes level splitting due to the sharp confinement of an impurity potential. The resulting energy levels for phosphorus are shown. (b) Our converged donor potential, including the central cell correction, which exhibits tetrahedral symmetry. The constant-energy surfaces shown are $-0.5$ (outer contour), $-1.0$ (middle contour), and $-4.0$ eV (central contour), respectively. (c), (d) Multivalley effective mass ground state for a single phosphorus donor in silicon. (c) shows a side view, while (d) shows a top-down view of the $x\text{−}y$ plane. The silicon lattice is superposed toward the center of the plots for scale; the white curtain indicates when the envelope ${\left|F\right|}^2$ is $1\%$ of its maximum value. (e), (f) Atomistic tight-binding simulations corresponding to (c), (d), performed in nemo-3d and visualized using the atomic orbitals of Ref. [33]. The envelope curtain is copied from (c) and (d) for comparison. (g), (h) Cuts along the parallel and perpendicular directions of the envelope function in one conduction band valley. The dashed lines are the effective mass theory from the present work; the shaded bands are $\pm 4\sigma$ statistical uncertainty limits, determined by the UQ techniques described in Appendix pp2. The lower bold curves show the corresponding Kohn-Luttinger envelope functions, for comparison. (i)–(k) Cuts along the $x\mathrm{axis}$ of the entire effective mass electron density for effective mass (solid curves) and nemo-3d [dotted curve in (i)]. (i) shows the $A_1$ ground state, (j) shows one of the three degenerate $T_2$ first excited states, and (k) shows one of the two degenerate $E$ first excited states.

Figure 2

The total Bloch function density, with the silicon lattice superimposed on the plots for scale. (a) shows well-converged Bloch functions, including high-frequency terms due to their periodic parts. (b) truncates to form factors, where each pair $u_{{\mathbf{k}}_0^l}^{*}\left(\mathbf{r}\right)u_{{\mathbf{k}}_0^j}\left(\mathbf{r}\right)$ is set equal to its constant-frequency component. (c) simplifies the situation further, using trivial [$u_{{\mathbf{k}}_0^j}\left(\mathbf{r}\right)=1$] Bloch functions. As shown in (b) and (c), these represent drastic approximations, so it is not surprising that calculations using them yield results different from those using the Bloch function in (a).

Figure 3

Tunnel couplings computed for two phosphorus donors in silicon. (a)–(d) Comparison of tunnel couplings computed within multivalley effective mass (points with error bars) and nemo-3d atomistic tight-binding theory (connected points with no error bars). Here, the tunnel coupling is defined as the energy difference between the first excited state and the ground state of the one-electron, two-donor problem. (a) shows tunnel coupling along the [100] direction, (b) along [110], and (c) along [111]. (d) depicts typical random instances, not along any particular direction. In all cases, the atomistic and effective mass theory exhibits very similar trends and magnitude of oscillations. Along [111], there appears to be a phase discrepancy, likely due to differing placements of the conduction band minima (see the main text for details). The error bars on the effective mass predictions are $\pm 4\sigma$ statistical uncertainty limits, determined by the UQ techniques described in Appendix pp2. (e) Exhaustive tunnel coupling enumeration for two phosphorus donors. Here, we placed one donor at the origin and the second at every possible point within a 30 nm cube surrounding it ($\sim 1.3\times {10}^6$ instances). The spherical shells show cuts (with nearest-neighbor interpolation) of the tunnel coupling at fixed donor separation distances. The tunnel coupling is highly oscillatory, and there is no large region of stability in the tunnel coupling. The full results of the enumeration are tabulated in the Supplemental Material [42].

Figure 4

Probability of achieving large tunnel coupling with uncertain donor placement. In each panel, one phosphorus donor is placed at the origin and a second is placed at lattice sites within a surrounding 30 nm cube. The placement of the second donor is uncertain. The plotted probability is that of obtaining $t>0.1$ meV for a Gaussian distribution of donor placements as a function of the distribution's center. The lower bound of $0.1$ meV is chosen to be about an order of magnitude larger than typical dilution refrigerator electron temperatures. We performed 20 000-shot Monte Carlo simulations, sampling from a three-dimensional (3D) isotropic Gaussian distribution with varying widths: (a) corresponds to 1 nm, (b) to 5 nm, and (c) to 10 nm straggle. (a) depicts an experimentally realistic straggle for donor placement based on scanning tunneling microscopy (STM), while (b) and (c) depict the results of increasing donor straggle. The white curtain shown in each plot indicates the contour of constant probability as labeled. These results show that STM placement can ensure large tunnel coupling with high yield, while ion implantation technology can only ever achieve low yield, rendering ion implantation ineffective for deterministic use.