# Spin-lattice relaxation times of single donors and donor clusters in silicon

## Abstract

An atomistic method of calculating the spin-lattice relaxation times () is presented for donors in silicon nanostructures comprising of millions of atoms. The method takes into account the full band structure of silicon including the spin-orbit interaction. The electron-phonon Hamiltonian, and hence the deformation potential, is directly evaluated from the strain-dependent tight-binding Hamiltonian. The technique is applied to single donors and donor clusters in silicon, and explains the variation of with the number of donors and electrons, as well as donor locations. Without any adjustable parameters, the relaxation rates in a magnetic field for both systems are found to vary as in excellent quantitative agreement with experimental measurements. The results also show that by engineering electronic wavefunctions in nanostructures, times can be varied by orders of magnitude.

###### pacs:

71.55.Cn, 03.67.Lx, 85.35.Gv, 71.70.EjDue to the extremely long spin coherence times, in some cases exceeding seconds (1); (2), and the existing industrial fabrication infrastructure, silicon is well-suited to be an outstanding platform for semiconductor quantum computer technology (3); (4); (5); (6); (7); (8). Qubits hosted by donors in silicon (3) have some added advantages as they are readily available few-electron systems with a rich electronic structure and can form identical qubits (9). In the last few years, several key experimental milestones have been achieved in dopant based quantum computing, including the demonstration of electron (10) and nuclear (11) spin qubits, single spin read-out and initialization (12); (13), and the observation of spin blockade and exchange towards two qubit coupling (14); (15). Recent advances in Scanning Tunneling Microscope (STM) lithography has enabled placement of single donors with atomic scale precision (16), with the result that various functional units such as quantum wires (17), single electron transistors (SET) (13); (18), and quantum dots (19) can all be realized in-plane with densely packed donor islands. The STM approach provides the fabrication precision needed to develop test-bed quantum chips for the demonstration of quantum algorithms in a solid-state quantum computer.

One of the two most important timescales for a spin qubit is the spin-lattice relaxation time (). Recent experiments have measured times in a single donor and in a few-donor cluster indicating shorter times in the latter (12); (13). Previous theoretical works exist in the literature qualitatively describing two different spin relaxation mechanisms in a bulk donor system (20); (21). However, a comprehensive quantitative theory which combines all the different mechanisms under a unified framework and accounts for the local inhomogeneous environment of the donors in a realistic nanostructure is still lacking. Moreover, there is no theoretical work yet to explain the measured times in densely packed donor clusters. In this letter, we present a comprehensive approach to compute the times in single donors and donor clusters in silicon nanostructures based on the self-consistent atomistic tight-binding (TB) method. The computed times can explain the experimental results of Refs (12); (13) without any adjustable parameters. The times are found to depend strongly on the size and shape of the electronic wavefunctions, which suggests that quantum confinement plays an important role in the relaxation process. The calculations also provide an insight into how the times can be engineered by several orders of magnitude in these quantum devices.

Single shot measurements of the donor-bound electron spin performed in ion-implanted metal-oxide-semiconductor (MOS) devices yielded a time of 2.3 s at T and a dependency of the relaxation rate () (12). Similar measurements have since been performed in a STM patterned donor cluster in an all epitaxially engineered device with atomic scale precision, giving time of 0.4 s at T and a similar dependency of the relaxation rate (13). The donor clusters are composed of several donors with an average separation less than or equal to the single donor Bohr radii (Fig. 1b) which can host several electrons. Such a cluster provides for additional addressability within an array of qubits (13). In general, the times of donors in realistic devices and nanostructures are likely to be influenced by the inhomogeneous environment, and need to be understood from an atomic scale theory.

The full TB Hamiltonian of the silicon and the P atoms was represented with a 20-orbital basis per atom including nearest-neighbor and spin-orbit interactions (23); (24); (25). A donor was represented by a Coulomb potential of a positive charge screened by the dielectric constant of Si and subjected to an onsite cut off potential (26); (27). The model reproduces the full energy spectrum of a single donor including valley-orbit splitting (22), and also captures the single donor hyperfine (26) and spin-orbit Stark effects (28) in close agreement with experiments. Recent STM imaging of the donor wavefunction also shows excellent agreement with the tight-binding wavefunction at the atomic scale (30). The magnetic field is represented by a vector potential in a symmetric gauge and entered through a Peierls substitution. To capture multi-electron occupation of the donor clusters, a self-consistent Hartree method was employed, in which the electronic charge density was computed from the lowest energy occupied wavefunctions, . Solving the potential due to self-consistently with the tight-binding Hamiltonian until convergence enables us to obtain the binding energies and the wavefunctions of the few electrons bound to the donor cluster. This method has also successfully reproduced the experimental (i.e., the 2nd electron bound to a single P) binding energy (29); (18). Exchange and correlation terms based on local density approximation are added to the Hartree potential (31). The same methodology has been used to reproduce experimentally measured addition energies of multi-donor clusters (15). The TB Hamiltonian of about 1.4 million atoms including the Hartree potential is then solved by a parallel Block Lanczos algorithm to obtain the relevant lowest energy wavefunctions.

For the relaxation times computed in this work, we assume that an electron has been loaded to the ground state of a single donor or a donor cluster. In case of a bulk P donor, this represents the state at -45.6 meV below the conduction band (32). This binding energy can vary with donor depths and fields in a realistic device (22). In addition, the binding energy in a donor cluster can be sensitive to the donor numbers, electron numbers, and donor locations. Experimentally, the electron can be loaded in the ground state of the donor/cluster by bringing the Fermi level of an electron reservoir in resonance with the Zeeman split ground state, as demonstrated in Refs (12); (13).

The relaxation rate , for an electron-phonon interaction Hamiltonian can be obtained by Fermi’s Golden Rule,

(1) |

where and are the up and down spin electronic states with energy and respectively, is the angular frequency of the emitted phonon, and are the initial and final phonon states with wavevector . In the B-field range of experimental interest (12); (13), the Zeeman splitting is less than 1 meV. As a result, only acoustic phonons contribute to the spin relaxation process. then depends on the deformation potential of the crystal (, representing each of the three Cartesian directions) and the strain tensor components (33) (both of which are position dependent in the atomistic TB method), and is given by

(2) |

To evaluate , we use the relation , and compute the total change in the electron-phonon Hamiltonian due to an infinitesimal uniform lattice strain represented by (a small arbitrary constant). Since can be expressed as a change in the electronic TB Hamiltonian under a crystal deformation caused by (33); (34), is given as

(3) |

The strain dependent TB Hamiltonian expresses the TB matrix elements as functions of inter-atomic bond lengths and distortion angles depending on the relative positions of pairs of atoms in the lattice (35). This method of incorporating local strain in the TB Hamiltonian is well-established and has been shown to reproduce various experimental results (23). Although we have considered all 6 components of , we have found the off-diagonal terms to be small for single donor states located near the conduction band.

Furthermore, can also be expressed in terms of phonon creation and annhilation operators, and respectively, as

(4) |

where is the volume of the crystal, the mass density, and the phonon polarization unit vector. Using eq. 2 and 4, the matrix element of can be expressed as,

(5) |

In eq. 5, we have used 2330 kg/, while is obtained from the electron Zeeman energy . The phonon number satisfies the Boltzmann distribution. We choose a temperature range below 100mK, which guarantees us to be in the low-temperature regime where (+1) (12). Above 1K, () , and varies as (20); (21). The polarization vector takes into account the three phonon modes: one longitudinal and two transverse. Since the Zeeman splitting energy is very small ( 1 meV), linear and isotropic bulk Si phonon dispersion is assumed. The phonon wavevector is then evaluated as , where is the speed of sound in Si, and taken to be 8480 m/s for the longitudinal mode and 5860 m/s for transverse modes. The same constants were also used to interpret the experimental data (12); (13). While local vibrational modes have been observed for P atoms in Si (36), the energy corresponds to the mode frequency is at least an order of magnitude larger than the energy range in interest. The measured values do not deviate from the behavior indicates that the local vibrational modes do not contribute significantly to the spin-relaxation process.

Fig. 2 shows the spin-lattice relaxation rates, , of a single P donor and a 4P donor cluster as a function of B-field. The red squares are the measured rates for a single donor from Ref. (12), whereas the blue triangles are the measured rates for a STM patterned few donor cluster from Ref. (13). Experimentally, the exact number of donors in the cluster was unknown, but estimated to be between 2 and 5 based on STM images. From transport measurements, it was also expected that during the spin readout step at least three electrons occupied the donor cluster, while in total seven charge transitions on the donor cluster were observed (13).

The red solid line in Fig. 2 represents the calculations performed in this work for a bulk P donor. The calculated rates show a dependence of , and also yield similar magnitudes of (2.5 s at 2 ) as the experiment. The B field in this calculation is applied along the [110] direction, consistent with the experiment (12). To understand the effect of donor number and electron number on , we have simulated donor clusters comprising of 2 to 4 donors with various electron occupation. In Fig. 2, we show the results of the 4P cluster with 1, 3 and 5 bound electrons (the green, cyan, and the blue solid lines, respectively). While the dependency holds in all cases, the rates vary considerably in magnitude, and increase with the number of electrons. Our calculations show that higher measured relaxation rates of the cluster come from a 5e occupation in a 4P cluster, which is also consistent with the experimental finding of the electron number being 3 (13).

We have intentionally chosen an odd number of electrons because the relaxation between a net 1/2 and -1/2 spin is assumed, which requires an unpaired number of electrons. This is also consistent with the experimental measurements, where no spin read-out signal was observed for alternate electron occupation in the cluster (13). The calculations also reveal the startling fact that if we have only one electron in a 4P cluster, the relaxation rate is actually smaller than the bulk P, and the times can be increased from few seconds to hundreds of seconds. The physical reason that determines the magnitude of the times is also investigated in this work.

To understand the impact of electron number on the times observed we determined the electron probability densities of the outermost electrons for varying cluster sizes and electron number in Fig. 3. The size and shape of the wavefunctions depends on the number of donors and electrons and their locations within the cluster. For the same number of electrons, more donors result in a more tightly bound wavefunction because of the stronger potential of the larger number of positively charged donor cores. For the same number of donors, as more electrons are added, the wavefunction spreads out more as the donor core becomes more strongly screened. Electron-electron repulsion also causes the wavefunctions of the outermost electron to spread out more. We have extracted the Bohr radii of these wavefunctions by fitting an exponential decay function to the tail of probability densities () along the x-axis through the center of the clusters.

Fig. 4a shows the relaxation rates (at 2 T) as a function of the Bohr radii for the same clusters as in Fig. 3. It is observed that donor clusters with larger Bohr radii (i.e., those clusters with more electrons and fewer donors) result in higher relaxation rates. Since the acoustic phonon wavelength corresponding to a Zeeman energy of 0.2 meV (at 2 T) is about 100 nm, the phonon wavelength is much larger than the electronic wavelengths in this system. A larger electronic wavefunction therefore interacts with the phonons more (38). Perhaps more importantly, a larger wavefunction also implies less quantum confinement in the system, and hence a smaller energy gap between excited and ground state (i.e. the valley-orbit gap) (39), which relates well to the Hasegawa theory for a bulk donor (20), in which increases with this energy gap. Since our calculations show a valley-orbit gap ranging from 30 meV to 5 meV as cluster wavefunction increases in radii, we expect a strong dependence of on radii as well. All the calculations presented here also include exchange and correlation effects. We observed that the inclusion of the exchange and correlation effects results in slightly larger wavefunctions, as the electrons experience greater net repulsion. This causes the relaxation rates to increase slightly and move closer to the experimental data in Fig. 2.

Within a lithographic template for a donor cluster, there can be some uncertainties in the exact locations of the donor (16). However, if the cluster only comprises of a few donors, all the positional configurations can be enumerated, and the most compact and the most dispersed donor clusters can be identified. Since it is computationally time intensive to simulate all possible donor configurations within a cluster, we have simulated the two extreme cases for donor clusters of 1 to 4 donors (Fig. 4c). In Fig. 4b, we have plotted the relaxation rates as a function of electron number for different donor clusters with these two positional configurations (marked by crosses within an ellipse). As expected, the times have some associated variations with donor locations within a cluster. However, the dependency on the total donor and total electron numbers is stronger. This suggests that measurements can be used to infer information about donor and electron numbers in STM patterned donor devices, as a non-invasive metrology technique. Such information can be useful for engineering pulses to control single- and two-qubit operations in experiments.

Previous effective mass calculations of Hasegawa (20) and Roth (21) predicted two different spin relaxation mechanisms due to an effective g-factor shift. Hasegawa’s mechanism predicts this effective g-factor shift due to the strain induced redistribution of the donor wavefunction among the six conduction band valleys (20), while Roth’s mechanism predicts a single valley g-factor shift due to a strain induced mixing of higher conduction bands (21). Both mechanisms inherently depend on the spin-orbit interaction in silicon, which reduces the bulk donor g-factor slightly below 2. Both Hasegawa and Roth’s theory show a dependency of for a bulk donor, but are rather qualitative in nature. The above approaches are also limited to a single bulk donor for which Kohn-Luttinger wavefunctions can be used (32). Our TB method captures both mechanisms under the same framework, as a full bandstructure description is used from the atomic orbital basis including spin-orbit interaction. The method is also general and applies to any nanostructures and semiconductors for which accurate TB models can be developed.

In conclusion, we have presented an atomistic approach to calculate the phonon induced spin lattice relaxation rates in donors in silicon. The times agree very well with recently measured values on single donors and donor clusters, and help to explain the variation of with the numbers of donors and electrons, and the donor locations. The values of were found to have a strong dependency on the size of the electronic wavefunctions. This also provides a way to engineer larger times by using donor clusters with large number of P donors and single electron. An atomistic description of the times and their variations in inhomogeneous environment provides crucial information in the design of silicon qubits.

###### Acknowledgements.

This research was conducted by the Australian Research Council Centre of Excellence for Quantum Computation and Communication Technology (project No. CE110001027), the US National Security Agency and the US Army Research Office under contract No. W911NF-08-1-0527. Computational resources on nanoHUB.org, funded by the NSF grant EEC-0228390, were used. M.Y.S. acknowledges a Laureate Fellowship.Electronic address: hsuehy@purdue.edu, rrahman@purdue.edu

### References

- A. M. Tyryshkin et al., Nature Materials 11, 143 (2011).
- M. Steger et al., Science 336, 1280 (2012).
- B. E. Kane et al., Nature 393, 133 (1998).
- L. C. L. Hollenberg et al., Phys. Rev. B 74, 045311 (2006).
- R. deSousa et al., Phys. Rev. A 70, 052304 (2004).
- M. J. Calderon et al., Phys. Rev. Lett. 96, 096802 (2006).
- M. Friesen et al., Phys. Rev. B 67, 121301 (2003).
- C. D. Hill et al., Phys. Rev. B 72, 045350 (2005).
- F. A. Zwanenburg et al., Rev. Mod. Phys. 85, 961 (2013).
- J. J. Pla et al., Nature 489, 541 (2012).
- J. J. Pla et al., Nature 496, 334 (2013).
- A. Morello et al., Nature 467, 687 (2010).
- H. Büch et al., Nature Communications 4, 2017 (2013).
- B. Weber et al., Nano letters 12, 4001 (2012).
- B. Weber et al., Nature Nanotechnology 9, 430 (2014).
- S. R. Schofield et al., Phys. Rev. Lett. 91, 136104 (2003).
- B. Weber et al., Science 335, 64 (2012).
- M. Fuechsle et al., Nature Nanotechnology 7, 242 (2012).
- M. Fuechsle et al., Nature Nanotechnology 5, 502 (2010).
- H. Hasegawa, Phys. Rev. 118, 1523 (1960).
- L. Roth, Massachusetts Institute of Technology Lincoln Laboratory Reports, April 1960 (unpublished).
- R. Rahman et al., Phys. Rev. B 80, 165314 (2009).
- G. Klimeck et al., IEEE Trans. Electron Dev. 54, 2079-2089 (2007).
- G. Klimeck et al., Computer Modeling in Engineering and Science Volume 3, No. 5, pp 601-642 (2002).
- J. C. Slater and G. F. Koster, Phys. Rev. 94, 1498 (1954).
- R. Rahman et al., Phys. Rev. Lett. 99, 036403 (2007).
- S. Ahmed et al., Springer Encyclopedia of Complexity and Systems Science, Springer-Verlag GmbH, Heidelberg, pp. 5745, ISBN: 978-0-387-75888-6 (2009).
- R. Rahman et al., Phys. Rev. B 80, 155301 (2009).
- R. Rahman et al., Physical Review B 84, 115428 (2011).
- J. Salfi et al., Nature Materials, DOI: 10.1038/NMAT3941 (2014).
- S. Lee et al., Phys. Rev. B 84, 205309 (2011).
- W. Kohn and J. M. Luttinger, Phys. Rev. 98, 915 (1955).
- B. K. Ridley, Quantum processes in semiconductors (Oxford University Press, 1982).
- J. M. Ziman, Electrons and Phonons (Oxford University Press, 1960).
- T. Boykin et. al., Phys. Rev. B 66, 125207 (2002).
- A. S. Barker, Jr. and A. J. Sievers, Rev. Mod. Phys. 47, S1 (1975).
- F. Mohiyaddin et al., Nanoletters 13, 1903 (2013).
- R. Hanson et al., Rev. Mod. Phys. 79, 1217 (2007).
- R. Rahman et al., Phys. Rev. B 83, 195323 (2011).