Spin-orbit coupling induced two-electron relaxation in silicon donor pairs

TL;DR

This paper theoretically identifies a dominant spin-orbit coupling-based relaxation mechanism in silicon donor pairs, revealing its dependence on donor configuration, temperature, and magnetic field, with implications for spintronics and quantum computing.

Contribution

It uncovers a previously overlooked relaxation process in silicon donor pairs driven by spin-orbit coupling, electron-phonon, and inter-donor interactions, with detailed symmetry analysis.

Findings

Relaxation rate scales from J^5 to J^4 with exchange interaction J.

Dependence on donor alignment and triplet spin orientation.

Relaxation occurs without magnetic field, with crossover behavior at high Zeeman energies.

Abstract

We unravel theoretically a key intrinsic relaxation mechanism among the low-lying singlet and triplet donor-pair states in silicon, an important element in the fast-developing field of spintronics and quantum computation. Despite the perceived weak spin-orbit coupling (SOC) in Si, we find that our discovered relaxation mechanism, combined with the electron-phonon and inter-donor interactions, dominantly drives the transitions in the two-electron states over a large range of donor coupling regime. The scaling of the relaxation rate with inter-donor exchange interaction $J$ goes from $J^5$ to $J^4$ at the low to high temperature limits. Our analytical study draws on the symmetry analysis over combined band, donor envelope and valley configurations. It uncovers naturally the dependence on the donor-alignment direction and triplet spin orientation, and especially on the dominant SOC source from donor impurities. While a magnetic field is not necessary for this relaxation, unlike in the single-donor spin relaxation, we discuss the crossover behavior with increasing Zeeman energy in order to facilitate comparison with experiments.

Tables

Table 1. Table 1: For each 𝐝 ^ ^ 𝐝 \hat{\mathbf{d}} , [001], [111] or [110], we find its point group and the IRs of the states and phonon modes. For 𝐝 ^ ∥ \hat{\mathbf{d}}\| [111], two subcases, with inversion symmetry ( i 𝑖 i ) and without it ( n 𝑛 n ), depend on the alternating donor positions (even though individual donors possess no inversion symmetry). To be definite, we set z ′ ∥ [ 111 ] conditional superscript 𝑧 ′ delimited-[] 111 z^{\prime}\|[111] , x ′ ∥ [ 10 − 1 ] conditional superscript 𝑥 ′ delimited-[] 10 1 x^{\prime}\|[10-1] , y ′ ∥ [ − 12 − 1 ] conditional superscript 𝑦 ′ delimited-[] 12 1 y^{\prime}\|[-12-1] , z ′′ ∥ [ 110 ] conditional superscript 𝑧 ′′ delimited-[] 110 z^{\prime\prime}\|[110] , x ′′ ∥ [ 001 ] conditional superscript 𝑥 ′′ delimited-[] 001 x^{\prime\prime}\|[001] and y ′′ ∥ [ 1 − 10 ] conditional superscript 𝑦 ′′ delimited-[] 1 10 y^{\prime\prime}\|[1-10] .

	[001] $\left(D_{2d}\right)$	$\begin{array}{c}\left[111\right]n\left(C_{3v}\right)\\ \left[111\right]i\left(D_{3d}\right)\end{array}$	[110] $\left(C_{2v}\right)$
$S$	$A_1$	$A_1$	$A_1$
$T_0$	$B_1$	$A_2\left(n\right),A_1^{\prime }\left(i\right)$	$A_2$
$T_{\pm }$	$E$	$E\left(n\right),E^{\prime }\left(i\right)$	$\begin{array}{c}{A}_1\left(T_{+}+T_{-}\right)\\ B_2\left(T_{+}-T_{-}\right)\end{array}$
	$A_1(\begin{array}{c}{ϵ}_{zz}\\ {ϵ}_{xx}+{ϵ}_{yy}\end{array})$	$A_1(\begin{array}{c}{ϵ}_{z^{\prime }{z}^{\prime }}\\ {ϵ}_{x^{\prime }{x}^{\prime }}+{ϵ}_{y^{\prime }{y}^{\prime }}\end{array})$	$A_1(\begin{array}{c}{ϵ}_{x^{\prime \prime }{x}^{\prime \prime }}\\ {ϵ}_{y^{\prime \prime }{y}^{\prime \prime }}\\ {ϵ}_{y^{\prime \prime }{y}^{\prime \prime }}\end{array})$
$\begin{array}{c}{ϵ}_{ij}\\ \end{array}$	$\begin{array}{c}{B}_1\left({ϵ}_{xx}-{ϵ}_{yy}\right)\\ B_2\left({ϵ}_{xy}\right)\\ E\left(\{{ϵ}_{xz},{ϵ}_{yz}\}\right)\end{array}$	$E\left(\begin{array}{c}\{{ϵ}_{x^{\prime }{x}^{\prime }}-{ϵ}_{y^{\prime }{y}^{\prime }},\\ 2{ϵ}_{x^{\prime }{y}^{\prime }}\}\\ \{{ϵ}_{x^{\prime }{z}^{\prime }},{ϵ}_{y^{\prime }{z}^{\prime }}\}\end{array}\right)$	$\begin{array}{c}{B}_1\left({ϵ}_{x^{\prime \prime }{z}^{\prime \prime }}\right)\\ B_2\left({ϵ}_{x^{\prime \prime }{y}^{\prime \prime }}\right)\\ A_2\left({ϵ}_{y^{\prime \prime }{z}^{\prime \prime }}\right)\end{array}$

Table 2. Table 2: Transformation matrices of spatial and TR symmetry operations for a two-donor system aligned along the [001] direction, i.e., D 2 ​ d subscript 𝐷 2 𝑑 D_{2d} group, upon T d subscript 𝑇 𝑑 T_{d} group IRs and spin basis states. For concreteness, we specify for the degenerate IRs that, E z I ∼ x 2 − y 2 similar-to subscript superscript 𝐸 𝐼 𝑧 superscript 𝑥 2 superscript 𝑦 2 E^{I}_{z}\sim x^{2}-y^{2} , E z I ​ I ∼ ( x 2 + y 2 − 2 ​ z 2 ) / 3 similar-to subscript superscript 𝐸 𝐼 𝐼 𝑧 superscript 𝑥 2 superscript 𝑦 2 2 superscript 𝑧 2 3 E^{I\!I}_{z}\sim(x^{2}+y^{2}-2z^{2})/\sqrt{3} ; T 1 ​ x ∼ R x , T 1 ​ y ∼ R y , T 1 ​ z ∼ R z formulae-sequence similar-to subscript 𝑇 1 𝑥 subscript 𝑅 𝑥 formulae-sequence similar-to subscript 𝑇 1 𝑦 subscript 𝑅 𝑦 similar-to subscript 𝑇 1 𝑧 subscript 𝑅 𝑧 T_{1x}\sim R_{x},T_{1y}\sim R_{y},T_{1z}\sim R_{z} ; and T 2 ​ x ∼ x , T 2 ​ y ∼ y , T 2 ​ z ∼ z formulae-sequence similar-to subscript 𝑇 2 𝑥 𝑥 formulae-sequence similar-to subscript 𝑇 2 𝑦 𝑦 similar-to subscript 𝑇 2 𝑧 𝑧 T_{2x}\sim x,T_{2y}\sim y,T_{2z}\sim z in terms of their transformation properties.

	$A_1$	$A_2$	$E^I,E^{II}$	$T_1^x,T_1^y,T_1^z$	$T_2^x,T_2^y,T_2^z$	${\uparrow }_z,{\downarrow }_z$
$C_{2z}$	1	1
${\rho }_{\left[1-10\right]}$	1	$-1$
${\rho }_{\left[110\right]}$	1	$-1$
$C_{2x}$	1	1
$C_{2y}$	1	1
$S_{4z}^{+}$	1	$-1$
$S_{4z}^{-}$	1	$-1$
TR	1	1

Table 3. Table 3: Transformation matrices of spatial and TR symmetry operations for a two-donor system aligned along the [111] direction, i.e., C 3 ​ v subscript 𝐶 3 𝑣 C_{3v} group without inversion or D 3 ​ d subscript 𝐷 3 𝑑 D_{3d} group with inversion (which connects α 𝛼 \alpha and β 𝛽 \beta donors and keeps the spin direction), upon T d subscript 𝑇 𝑑 T_{d} group IRs and spin basis states.

	$A_1$	$A_2$	$E^I,E^{II}$	$T_1^x,T_1^y,T_1^z$	$T_2^x,T_2^y,T_2^z$	${\uparrow }_{\left[111\right]},{\downarrow }_{\left[111\right]}$
$C_{3\left[111\right]}^{+}$	1	1
$C_{3\left[111\right]}^{-}$	1	1
${\rho }_{\left[1-10\right]}$	1	$-1$
${\rho }_{\left[10-1\right]}$	1	$-1$
${\rho }_{\left[01-1\right]}$	1	$-1$
$i$	1	$-1$
TR	1	1

Table 4. Table 4: Transformation matrices of spatial and TR symmetry operations for a two-donor system aligned along the [110] direction, i.e., C 2 ​ v subscript 𝐶 2 𝑣 C_{2v} group, upon T d subscript 𝑇 𝑑 T_{d} group IRs and spin basis states.

	$A_1$	$A_2$	$E^I,E^{II}$	$T_1^x,T_1^y,T_1^z$	$T_2^x,T_2^y,T_2^z$	${\uparrow }_{\left[110\right]},{\downarrow }_{\left[110\right]}$
$C_{2z}$	1	1
${\rho }_{\left[110\right]}$	1	-1
${\rho }_{\left[1-10\right]}$	1	-1
TR	1	1

Equations

$$H_{ee}=\sum _{i=1,2}[K({\mathbf{r}}_i)+V({\mathbf{r}}_i)]+V_{ee}({\mathbf{r}}_1,{\mathbf{r}}_2),$$

$$K=-\frac{{\hslash }^2{\nabla }^2}{2m},V=V_{\mathrm{cryst}}+\sum _{j=\alpha ,\beta }\delta V_{\mathrm{ion},j},V_{ee}=\frac{e^2}{ϵr_{12}},$$

$$H_{\mathrm{ep}}=-\sum _k\delta {\mathbf{R}}_k\cdot \nabla V(\mathbf{r}-{\mathbf{R}}_k),$$

$$S$$

$$T_0$$

$$T_{\pm }$$

$$⟨T_0|{ϵ}_{xx}-{ϵ}_{yy}|S{⟩}_{[001]},⟨T_{\pm }|{ϵ}_{xz}\mp (\pm )i{ϵ}_{yz}|S(T_0){⟩}_{[001]},$$

$$⟨T_{\pm }|{ϵ}_{x^{\prime }{z}^{\prime }}\pm i{ϵ}_{y^{\prime }{z}^{\prime }},{ϵ}_{x^{\prime }{x}^{\prime }}-{ϵ}_{y^{\prime }{y}^{\prime }}\pm 2i{ϵ}_{x^{\prime }{y}^{\prime }}|S,T_0{⟩}_{[111]n},$$

$$⟨T_{\pm }|{ϵ}_{x^{\prime }{z}^{\prime }}\pm i{ϵ}_{y^{\prime }{z}^{\prime }},{ϵ}_{x^{\prime }{x}^{\prime }}-{ϵ}_{y^{\prime }{y}^{\prime }}\pm 2i{ϵ}_{x^{\prime }{y}^{\prime }}|T_0{⟩}_{[111]i},$$

$$⟨T_0|{ϵ}_{y^{\prime \prime }{z}^{\prime \prime }}|S,T_1{⟩}_{[110]},⟨T_2|{ϵ}_{x^{\prime \prime }{z}^{\prime \prime }}|T_0{⟩}_{[110]},$$

$$⟨T_1|{ϵ}_{x^{\prime \prime }{x}^{\prime \prime }},{ϵ}_{z^{\prime \prime }{z}^{\prime \prime }},{ϵ}_{y^{\prime \prime }{y}^{\prime \prime }}|S{⟩}_{[110]},⟨T_2|{ϵ}_{x^{\prime \prime }{y}^{\prime \prime }}|S,T_1{⟩}_{[110]}.$$

$$⟨T_0|{ϵ}_{xx}-{ϵ}_{yy}|S{⟩}_{[001]}$$

$$⟨T_0|{ϵ}_{y^{\prime \prime }{z}^{\prime \prime }}|S{⟩}_{[110]}\propto \sum _{\begin{array}{c}\updownarrow =\uparrow ,\downarrow \\ \gamma =\alpha ,\beta \end{array}}(-1{)}^{\gamma }⟨{\psi }_{\gamma {\updownarrow }^{\prime \prime }}^{T_0}|{ϵ}_{y^{\prime \prime }{z}^{\prime \prime }}|{\psi }_{\alpha {\uparrow }^{\prime \prime }}^S⟩⟨{\psi }_{\overline{\gamma }{\overline{\updownarrow }}^{\prime \prime }}^{T_0}|{\psi }_{\beta {\downarrow }^{\prime \prime }}^S⟩\text{ (by σy′′, C2x′′ and σz′′)},$$

$$⟨T_{+}|{ϵ}_{xz}\mp i{ϵ}_{yz}|S/T_0{⟩}_{[001]}\propto \frac{\mp 1}{\sqrt{2}}\sum _{\gamma =\alpha ,\beta }(-1{)}^{\gamma }⟨{\psi }_{\gamma \uparrow }^{T_{+}}|{ϵ}_{xz}\mp i{ϵ}_{yz}|{\psi }_{\alpha \downarrow }^{S/T_0}⟩⟨{\psi }_{\overline{\gamma }\uparrow }^{T_{+}}|{\psi }_{\beta \uparrow }^{S/T_0}⟩\text{ (by S4z and C2z)},$$

$$⟨T_{+}|{ϵ}_{x^{\prime }{z}^{\prime }}+i{ϵ}_{y^{\prime }{z}^{\prime }}|S/T_0{⟩}_{[111]n}\propto \frac{\mp 1}{2\sqrt{2}}\sum _{\gamma ,{\gamma }^{\prime }=\alpha ,\beta }(-1{)}^{\gamma +\frac{{\gamma }^{\prime }}{2}\mp \frac{{\gamma }^{\prime }}{2}}⟨{\psi }_{\gamma {\uparrow }^{\prime }}^{T_{+}}|{ϵ}_{x^{\prime }{z}^{\prime }}+i{ϵ}_{y^{\prime }{z}^{\prime }}|{\psi }_{{\gamma }^{\prime }{\downarrow }^{\prime }}^{S/T_0}⟩⟨{\psi }_{\overline{\gamma }{\uparrow }^{\prime }}^{T_{+}}|{\psi }_{\overline{{\gamma }^{\prime }}{\uparrow }^{\prime }}^{S/T_0}⟩\text{(by C3z′)},$$

$$⟨T_{+}|{ϵ}_{x^{\prime }{z}^{\prime }}+i{ϵ}_{y^{\prime }{z}^{\prime }}|T_0{⟩}_{[111]i}\propto \frac{1}{\sqrt{2}}\sum _{\gamma =\alpha ,\beta }(-1{)}^{\gamma }⟨{\psi }_{\gamma {\uparrow }^{\prime }}^{T_{+}}|{ϵ}_{x^{\prime }{z}^{\prime }}+i{ϵ}_{y^{\prime }{z}^{\prime }}|{\psi }_{\alpha {\downarrow }^{\prime }}^{T_0}⟩⟨{\psi }_{\overline{\gamma }{\uparrow }^{\prime }}^{T_{+}}|{\psi }_{\beta {\uparrow }^{\prime }}^{T_0}⟩\text{ (by C3z′ and i)},$$

$$⟨T_{+}\pm T_{-}|{ϵ}_{A_1}/{ϵ}_{x^{\prime \prime }{y}^{\prime \prime }}|S{⟩}_{[110]}\propto \sum _{\begin{array}{c}\updownarrow =\uparrow ,\downarrow \\ \gamma =\alpha ,\beta \end{array}}(-1{)}^{\updownarrow +\gamma }⟨{\psi }_{\gamma {\uparrow }^{\prime \prime }}^{T_{1/2}}|{ϵ}_{A_1/x^{\prime \prime }{y}^{\prime \prime }}|{\psi }_{\alpha {\updownarrow }^{\prime \prime }}^S⟩⟨{\psi }_{\overline{\gamma }{\uparrow }^{\prime \prime }}^{T_{1/2}}|{\psi }_{\beta {\overline{\updownarrow }}^{\prime \prime }}^S⟩\text{ (by σz′′ and σy′′)},$$

$$⟨T_{+}\pm T_{-}|{ϵ}_{y^{\prime \prime }{z}^{\prime \prime }}/{ϵ}_{x^{\prime \prime }{z}^{\prime \prime }}|T_0{⟩}_{[110]}\propto \sum _{\begin{array}{c}\updownarrow =\uparrow ,\downarrow \\ \gamma =\alpha ,\beta \end{array}}(-1{)}^{\gamma }⟨{\psi }_{\gamma {\uparrow }^{\prime \prime }}^{T_{1/2}}|{ϵ}_{y^{\prime \prime }{z}^{\prime \prime }/x^{\prime \prime }{z}^{\prime \prime }}|{\psi }_{\alpha {\updownarrow }^{\prime \prime }}^{T_0}⟩⟨{\psi }_{\overline{\gamma }{\uparrow }^{\prime \prime }}^{T_{1/2}}|{\psi }_{\beta {\overline{\updownarrow }}^{\prime \prime }}^{T_0}⟩\text{ (by σz′′ and σy′′)},$$

$$⟨T_{+}-T_{-}|{ϵ}_{x^{\prime \prime }{y}^{\prime \prime }}|T_{+}+T_{-}{⟩}_{[110]}\propto \sum _{\begin{array}{c}\updownarrow =\uparrow ,\downarrow \\ \gamma =\alpha ,\beta \end{array}}(-1{)}^{\gamma +\updownarrow }⟨{\psi }_{\gamma {\updownarrow }^{\prime \prime }}^{T_2}|{ϵ}_{x^{\prime \prime }{y}^{\prime \prime }}|{\psi }_{\alpha {\uparrow }^{\prime \prime }}^{T_1}⟩⟨{\psi }_{\overline{\gamma }{\updownarrow }^{\prime \prime }}^{T_2}|{\psi }_{\beta {\uparrow }^{\prime \prime }}^{T_1}⟩\text{ (by σz′′ and σy′′)},$$

$$⟨{\psi }_{\alpha (\beta )\uparrow }^{T_0}|{ϵ}_{xx}-{ϵ}_{yy}|{\psi }_{\alpha \uparrow }^S⟩=-⟨{\psi }_{\alpha (\beta )\downarrow }^{T_0}|{ϵ}_{xx}-{ϵ}_{yy}|{\psi }_{\alpha \downarrow }^S⟩$$

$$\{{\psi }_k\}=\sum _{i=1}^6{v}_{k,i}{F}_{k,i}^{nlm}(\mathbf{r}){\psi }_{k,i}^{{\Delta }_j}(\mathbf{r}),$$

$${\alpha }_{\uparrow }$$

$$⟨T_0|{ϵ}_{xx}-{ϵ}_{yy}|S{⟩}_{[001]}^{(1)}\propto \sum _{\nu =1s,3d_{\pm 1}}\frac{⟨{\alpha }_{{\Delta }_1,1s,A_1}^{T_0}|{\overline{H}}_{\mathrm{int}-\mathrm{d}}|{\alpha }_{{\Delta }_1,\nu ,T_{2z}}^{T_0}⟩⟨{\alpha }_{{\Delta }_1,\nu ,T_{2z}}|L_z|{\alpha }_{{\Delta }_1,1s,E_z^I}⟩⟨{\alpha }_{{\Delta }_1,1s,E_z^I}|{ϵ}_{xx}-{ϵ}_{yy}|{\alpha }_{{\Delta }_1,1s,A_1}⟩}{({\mathcal{E}}_{{\Delta }_1,1s,A_1}-{\mathcal{E}}_{{\Delta }_1,\nu ,T_2})({\mathcal{E}}_{{\Delta }_1,1s,A_1}-{\mathcal{E}}_{{\Delta }_1,1s,E})},$$

$$M_{\kappa ,\lambda }$$

$$[001]:\alpha \uparrow$$

$$\beta \downarrow$$

$$\alpha \downarrow$$

$$\beta \uparrow$$

$$[111]:\alpha \uparrow$$

$$\alpha \downarrow$$

$$\alpha$$

$$[110]:\alpha \uparrow$$

$$\beta \downarrow$$

$$\alpha \downarrow$$

$$\beta \uparrow$$

$$S(T_0{)}^0:\frac{(1-P_{12})}{2\sqrt{1\pm {\chi }_{A_1}^2}}[{\alpha }_{A_1\uparrow }({\mathbf{r}}_1){\beta }_{A_1\downarrow }({\mathbf{r}}_2)\mp {\alpha }_{A_1\downarrow }({\mathbf{r}}_1){\beta }_{A_1\uparrow }({\mathbf{r}}_2)],$$

$$S(T_0{)}^{\prime }:\frac{(1-P_{12})}{2\sqrt{2(1\pm {\chi }_{A_{1}X}^2\pm {\chi }_{A_1}{\chi }_X)}}[e^{ia}{\alpha }_{A_1\uparrow }({\mathbf{r}}_1){\beta }_{X\downarrow }({\mathbf{r}}_2)+e^{ib}{\alpha }_{X\uparrow }({\mathbf{r}}_1){\beta }_{A_1\downarrow }({\mathbf{r}}_2)\mp e^{ic}{\alpha }_{A_1\downarrow }({\mathbf{r}}_1){\beta }_{X\uparrow }({\mathbf{r}}_2)\mp e^{id}{\alpha }_{X\downarrow }({\mathbf{r}}_1){\beta }_{A_1\uparrow }({\mathbf{r}}_2)],$$

$${\mathcal{E}}_{S(T_0{)}^0}$$

$${\mathcal{E}}_{S(T_0{)}^{\prime }}$$

$$⟨S(T_0{)}^{\prime }|S(T_0{)}^0⟩$$

$$⟨S(T_0{)}^{\prime }|H_{ee}|S(T_0{)}^0⟩$$

$$T_1$$

Full text

I Introduction

A pair of coupled donors in an enriched 28Si host crystal forms the most cleanly defined two-qubit system or singlet-triplet ($S$-$T$) qubit in solid state materials Zwanenburg_RMP13 ; Shankar_PRB15 . The nearly noise-free crystal environment Tyryshkin_PRB03 ; Muhonen_NatNano14 ; Itoh_MRS14 together with the reproducible donor properties fixed by nature Kohn_SSP57 places an essential role for the intrinsic energy relaxation among the lowest two-electron states by electron-phonon (e-ph) interaction in studying the coherence limit of quantum computation involving such qubits Borhani_PRB10 ; Hsueh_PRL14 . A proper treatment of this problem with a transparent and physical basis is imperative, enabling its various modifications such as proximity with interface Calderon_PRB07 , gate potential Dehollain_PRL14 ; Gonzales-Zalba_NanoLett14 , or quantum dots Foote_APL15 ; Harvey-Collard_arxiv15 which are necessary for qubit control and operation.

Electron spin relaxation of single donor states in Si has been studied since the 1950s. Much interest then was associated with the ensemble spin resonance experiments which could remarkably map out the detailed donor wavefunction by its hyperfine interaction with the donor or 29Si nucleus Fletcher_PR54 ; Feher_PR59 . Theoretical attention was thus originally paid to the simultaneous flip of electron and donor nucleus spins ($T_{X}$) Pines_PR57 through hyperfine interaction, until later the importance of spin-orbit coupling (SOC) driven spin relaxation was established Hasegawa_PR60_Roth_PR60_Roth_LincolnLab60 ; Wilson_PR61 ; Castner_PR63 ; Castner_PR67 . However, the two-electron relaxation on a coupled donor pair has been much less studied over the decades, with only a few works in the context of concentration-dependent spin relaxation Sugihara_JPCS68 or of quantum computation Borhani_PRB10 . Both of these previous studies focused on the hyperfine interaction to mix $S$ and $T$ states, while neglecting SOC as too weak in Si.

Here we unravel a distinct general relaxation mechanism among two-electron states to shed important light on this long-standing physical problem (i.e., electron relaxation in high density coupled donors), motivated by the development of Si donor-based spin qubits. This robust mechanism has a parametrically dominating $J^{5}$ dependence on the 2-donor exchange coupling $J$ for the $S\leftrightarrow T$ transition, and becomes always dominant at sufficiently strong coupling. This mechanism relies only on the intrinsic SOC in the system to couple spins, with the donor SOC dominating the Si host SOC, but not on any structure-induced Rashba or Dresselhaus field. This particular dominance of the donor SOC is absent in the classical single-donor spin relaxation by $g$-factor modulation Hasegawa_PR60_Roth_PR60_Roth_LincolnLab60 . The effectiveness of the two SOCs, i.e., of the host and the impurity, is considered by taking full account of the three physical ‘layers’ in the problem: the Bloch bands, the hydrogenic-like donor envelopes, and the multi-valley configurations. Within each layer rigourous selection rules are enforced for e-ph, SOC and inter-donor perturbations, while at the same time, the decoupling between the underlying fast-oscillatory physics and the slowly varying localized envelopes realizes an efficient leading-order quantitative estimation. The unique donor and spin alignment dependence naturally falls out of our analytical treatment. Finally, we elucidate the crossover from low to high Zeeman energy as a reduction of the two-electron problem into a single-electron one, although the presence of spin splitting is not a priori essential for our relaxation mechanism.

The rest of the paper is organized as follows. In Sec. II, we lay out the theoretical framework, and focus on developing a suitable strategy in order to have a clear insight and efficient calculation of the two-electron state relaxation matrix elements. The general idea follows the tradition of utilizing symmetry properties to reduce the complexity as much as possible, and to proceed with the physically correct grasp (not to overlook dominant contributions or cancellations) using perturbative approaches. To keep the flow of theoretical development, we defer several particularly long derivation details to the Appendices. In Sec. III, we make approximate quantitative evaluations for various relaxation transitions and lifetimes based on well-known physical parameters of Si band and donor energy structures, phonon dispersion, and deformation potentials. We also analyze the relaxation trend in two-electron states with external magnetic field going from zero to large values. A summary is given in Sec. IV.

II Theoretical derivation

To properly treat the relaxation problem driven by SOC in a combined system of Si host and coupled double donors, the potential energy part of the Hamiltonian ($V$) for the donor outer electrons includes both the Si crystal potential ($V_{\rm cryst}$) and the difference between Si atoms and the substitutional donor ions ($\delta V_{{\rm ion},j}$, where $j=\alpha,\beta$ labels the donors). Combined with the kinetic energy ($K$) and electron-electron interaction ($V_{ee}$), the Hamiltonian for the electron states in the absence of e-ph interaction is,

[TABLE]

where $r_{12}\equiv|\mathbf{r}_{1}-\mathbf{r}_{2}|$, $\epsilon$ is the dielectric constant in Si, $V_{\rm cryst}=V_{{\rm cryst},0}+V_{\rm cryst,so}$ contains spin-independent and SOC parts, and $\delta V_{{\rm ion},j}(\mathbf{r})=-e^{2}/\epsilon|\mathbf{r}-\mathbf{R}_{j}|+\delta V_{\rm ion,cc}(\mathbf{r}-\mathbf{R}_{j})$ contains the Coulomb potential as well as the short-range central-cell correction at the core of the ion. $\delta V_{\rm ion,cc}=\delta V_{\rm ion,cc,0}+\delta V_{\rm ion,cc,so}$ and $\delta V_{\rm ion,cc,so}$ easily dominates the donor SOC effect due to the fast changing core potential Weisz_PR66 , compared with that from the smooth Coulomb potential away from the central cell region. The e-ph interaction that induces the electronic transition has the general form within the harmonic approximation,

[TABLE]

where the lattice displacement at atom $k$ is a superposition of phonon states, $\delta\mathbf{R}_{k}=\sum_{\mathbf{q}}\sqrt{\hbar/2\omega^{\lambda}_{\mathbf{q}}M}(a^{\lambda}_{\mathbf{q}}\bm{\xi}^{\lambda}_{k,\mathbf{q}}e^{i\mathbf{q}\cdot\mathbf{R}_{k}})+$H.c, with $M$ being the mass of Si crystal, $\omega^{\lambda}_{\mathbf{q}}$, $a^{\lambda}_{\mathbf{q}}$ and $\bm{\xi}^{\lambda}_{k,\mathbf{q}}$ the eigen-frequency of the $\lambda$th phonon band at wavevector $\mathbf{q}$, the associated annihilation operator and the normalized phonon polarization vector, respectively Maradudin_SSP63 . Here we use the long-wavelength acoustic phonons in bulk Si which can compensate the small energy difference ($\lesssim$ meV) in the two-electron states. As we will see, these two-electron relaxations require the third-order perturbation theory even using the bound Bloch states as basis. As a result, the relevant intermediate states encompass both large and small energy scales due to the various fast (SOC and e-ph) and slow (exchange) perturbative terms. We are forced to go beyond the usual conductance-band-bottom effective mass approximation Kohn_SSP57 to include other Bloch bands, in addition to the central cell correction we have explicitly taken into account. This study thus represents a complicated single-particle problem which combines crystal and atomic physics along with a relatively simple many-body (two-body) problem. It is important to choose a systematic treatment for efficient calculations and avoiding pitfalls.

We start with a set of two-electron representations satisfying the symmetry of the specific two-donor system. While these states mainly comprise the Heitler-London states out of unperturbed single-donor ground states, both nonpolar and polar (‘ionic’) mixtures are included. These excited components are needed for incorporating nonzero e-ph and SOC coupling, and are mixed in by the inter-donor interaction. The dominant effect comes from the leading-order mixture where one of the two electrons occupies the single-donor ground state. As a result, the general states follow the concise expressions,

[TABLE]

where $\mathcal{P}_{12}$ exchanges electrons 1 and 2, the one-electron wavefunction $\psi_{\alpha(\beta)}$ is mainly located around donor $\alpha(\beta)$, $\chi=\langle\psi_{\alpha}|\psi_{\beta}\rangle$, superscripts ($S,T_{0,\pm}$) distinguish different mixture components which are described in details later [Eq. (22)], and finally spins $\uparrow$ and $\downarrow$ (which generally are quasi-spins that contain small opposite spin component due to SOC-induced mixing) are along the $\alpha-\beta$ donor-alignment direction ($\hat{\mathbf{d}}$), dictated by symmetry.

We next classify the two-electron states and the e-ph interaction by the total system symmetries for three crystallographic directions $\hat{\mathbf{d}}$, listed in Table 1. For brevity, we use strain element $\epsilon_{ij}$ in place of the e-ph interaction $\sim[r_{i}dV/dr_{j}+(i\leftrightarrow j)]/2$ which transform the same way symmetry-wise. For $\hat{\mathbf{d}}\|$[110], $S$ or $T_{-}+T_{+}$ denotes the eigenstate with the majority $S$ or $T_{-}+T_{+}$ component respectively since they belong to the same irreducible representation (IR) .

This top-down symmetry approach immediately identifies the allowed relaxation processes and the associated specific acoustic phonon modes, independent of the quantitative treatments one adopts. The nonvanishing phonon-induced relaxation matrix elements are,

[TABLE]

where $T_{1/2}\equiv T_{+}\pm T_{-}$ in Eqs. (10) and (11). We also have taken into account the fact that $S$ and $T$ states (e-ph operator) are two(one)-body objects, removing all $\langle T_{-}|\epsilon_{ij}|T_{+}\rangle$ where $\langle\psi^{T_{-}}_{\alpha(\beta)\downarrow}|\psi^{T_{+}}_{\alpha\uparrow}\rangle$ are always strictly zero, solely by the $C_{2(3)}$ operation in the $\hat{\mathbf{d}}\|$001 case.

To evaluate the allowed relaxation matrix elements, the first simplifications come from substituting the forms of $S$ and $T$ states with Eqs. (4)-(6). Even though the specific $\psi$’s are altered away from the single-donor wavefunction by inter-donor interaction, they still obey precise symmetry relations among themselves due to the point group of a given two-donor system. Utilizing these relations, the matrix elements are reduced into products of single-electron ones. A representative set of steps is shown for $\langle T_{0}|\epsilon_{xx}-\epsilon_{yy}|S\rangle_{[001]}$,

[TABLE]

where for brevity we use ‘$\propto$’ and omit the normalization factor $1/\sqrt{1\pm\chi^{2}}$, $\overline{\gamma}$ or $\overline{\updownarrow}$ denotes the opposite donor or spin respectively, and at exponent $\alpha,\uparrow\equiv 0$ and $\beta,\downarrow\equiv 1$. The key symmetry operations used for the reduction are explicitly marked in the parenthesis. Similarly, we can substitute the expressions from Eqs. (4)-(6) in the rest of relaxation matrix elements in Eqs. (7)-(11), and utilize available symmetry operations to make the simplifications.

[TABLE]

where $\uparrow^{\prime}(\uparrow^{\prime\prime})$ indicates that the spin is along the $z^{\prime}(z^{\prime\prime})$ rather than $z$ direction, and again $T_{1/2}\equiv T_{+}\pm T_{-}$ in Eqs. (17)-(19). The transition matrix element magnitudes from $S$ or $T_{0}$ to $T_{-}$ state are the same as those to $T_{+}$ as obtained by switching between $\epsilon_{xz}\!\pm\!i\epsilon_{yz}$ ($\hat{\mathbf{d}}\|[001]$) or $\epsilon_{x^{\prime}z^{\prime}}\!\pm\!i\epsilon_{y^{\prime}z^{\prime}}$ ($\hat{\mathbf{d}}\|[111]$). In $\hat{\mathbf{d}}\|[111]$, one can just substitute $\epsilon_{x^{\prime}z^{\prime}}\!\pm\!i\epsilon_{y^{\prime}z^{\prime}}$ with $\epsilon_{x^{\prime}x^{\prime}}-\epsilon_{y^{\prime}y^{\prime}}\pm 2i\epsilon_{x^{\prime}y^{\prime}}$. We note that the additional inversion in the $[111]i$ case equates each of the two pairs in Eq. (15) of $[111]n$, canceling $\langle T_{\pm}|\epsilon_{x^{\prime}z^{\prime}}\!\!\pm\!i\epsilon_{y^{\prime}z^{\prime}}\!|\!S\rangle_{\![111]i}$ while leading to Eq. (16). This set of equations arising from very general symmetry considerations constitutes one of the key results in this work.

We apply the perturbation theory to quantify the single-electron matrix elements in terms of SOC and exchange coupling constants, in addition to the deformation potential coupling (i.e., e-ph interaction with acoustic phonons). The first two perturbations are necessary, as without them the transitions reduce to those between pure opposite spins or single-donor spin relaxation which must vanish. In particular, the same-spin e-ph matrix element in Eq. (12) from $\langle T_{0}|\epsilon_{xx}-\epsilon_{yy}|S\rangle_{[001]}$ requires a $z$-component SOC operator, since

[TABLE]

by the $\sigma_{[110]}$ reflection symmetry. It is similar for the same-spin transitions of $\hat{\mathbf{d}}\|[110]$ in Eqs. (13), (17)-(19) except with the $\epsilon_{A_{1}}$ or $\epsilon_{x^{\prime\prime}z^{\prime\prime}}$ modes, while all the rest manifestly require SOC to flip the spin.

To proceed within the perturbation theory, we choose our basis states to be the spinless donor states, i.e., the eigenstates of $H_{e,0}=K+V_{\rm cryst,0}-e^{2}/\epsilon r$,

[TABLE]

which are identified by three indices: the Si bulk band ($\Delta_{1},\Delta_{2^{\prime}},\Delta_{2},\Delta_{1^{\prime}}$ and $\Delta_{5}$) Lax_PR61 , the donor envelope with an orbital number $nlm$ (with ellipsoidal effective mass) Kohn_PR55 , and the $T_{d}$ (tetrahedral) group IR ($A_{1},A_{2},E,T_{1}$ and $T_{2}$) Bradley_Cracknell72 which determines $v_{k,i}$ ($\sum_{i}|v_{k,i}|^{2}=1$) considering the participating $\Delta_{j}$ and $nlm$. These three indices fix the energy level in a roughly descending order. The three perturbations in addition to $H_{e,0}(\mathbf{r}_{1}-\mathbf{r}_{\alpha})+H_{e,0}(\mathbf{r}_{2}-\mathbf{r}_{\beta})$ can now be explicitly seen in $H_{ee}+H_{\rm ep}$ [from Eqs. (1) and (3)]: the SOC part, $V_{\rm crysta,so}$; the obvious e-ph part, $H_{\rm ep}$; and the inter-donor part, $H_{\rm int-d}=V_{ee}+\delta V_{\rm ion,\alpha}(\mathbf{r}_{2})+\delta V_{\rm ion,\beta}(\mathbf{r}_{1})$, which is a two-body interaction and its matrix element can be integrated over one variable (e.g., $\mathbf{r}_{2}$) to obtain the mixture components it contributes to a single electron wavefunction [e.g., $\psi(\mathbf{r}_{1})$ in Eqs. (4)-(6)]. The effective single-electron interaction from the inter-donor coupling, defined in this way, is denoted $\overline{H}_{\rm int-d}$, as used in Eq. (23). Before exhaustively working out all possible selection rules among this multitude of states, we examine the essential physics of coupling strengths for different perturbation interactions and select the stronger couplings efficiently.

First, we focus on the inter-donor interaction, including direct Coulomb and exchange terms. It couples donor states made of different bulk $\Delta$ bands very weakly. As we know, in single donor ground states, bands other than the conduction $\Delta_{1}$ band are routinely neglected due to their fast-oscillating difference and the slowly varying nature of the Coulomb interaction Luttinger_PR55 ; Kohn_SSP57 . Here the coupling by inter-donor interaction is even weaker as the inter-donor distance is several times the Bohr radius. Within the same bulk band, it can couple different donor envelopes as well as valley configurations effectively, as their differences are (partly) slowly varying. Second, the e-ph coupling between the same or different $\Delta_{j}$’s are efficient when allowed by symmetry, as the interaction involves periodic ion potentials. This gives rise to various intraband and interband deformation potentials. Once the interaction matches the symmetry difference of the two bulk bands, it can only couple the same envelopes as no extra symmetry from the phonon mode compensates for the different envelope symmetries. However, it may couple different valley configurations as the intravalley e-ph coupling may change from valley to valley. Last, we discuss the SOC of two different types arising from the host and the donor footnote_2SOC . For Si host SOC, it couples symmetry-allowed $\Delta$ bands strongly but different envelopes negligibly just like the e-ph coupling. The donor impurity SOC may couple donor envelopes in the same band effectively Sugihara_JPSJ63 ; Shimizu_JPSJ70 ; Castner_PR67 . Additionally, both SOCs can couple different valley configurations allowed by symmetry. However, the host SOC together with the e-ph coupling connecting two same $\Delta$ bands (imposed by the inter-donor coupling) largely recovers the “Elliott-Yafet” (E-Y) cancellation that occurs in the bulk Si spin-phonon interaction Yafet_SSP63 ; Song_PRB12 , and suppresses their effect by a large factor of the relevant phonon wavelength divided by lattice constant. Without this suppression, the spin mixing caused by the host SOC would be about $\Delta^{\rm hst}_{\rm SOC}/\Delta\mathcal{E}^{\rm hst}\sim$ 40 meV/ 4eV (interband coupling) and comparable to that by the P donor SOC, $\Delta^{\rm dnr}_{\rm SOC}/\Delta\mathcal{E}^{\rm dnr}\sim$ 0.03 meV/12 meV (inter-bound-state coupling). For more details see Eq. (23). The spin splitting in the donor states, on the other hand, depends strongly on the donor types (P, As, Sb) and comes largely from the donor SOC Aggarwal_PR65 ; Castner_PR67 . It measures the SOC contribution from the $T_{d}$ potential deviating from a spherical one, i.e, the SOC contribution that breaks the inversion symmetry and hence annuls the first-order-in-wavevector E-Y cancellation Yafet_SSP63 ; Song_PRB12 .

Understanding the potentially dominant couplings, one can identify the symmetry-allowed ones among them. The e-ph and SOC couplings follow conventional single-particle selection rules, as discussed below. The two-body inter-donor interaction, however, requires a different treatment. We identify the allowed $T_{d}$ IRs for the mixed single-donor states, such that under every two-donor symmetry operation the resulting $\psi_{\gamma\updownarrow}$ transforms the same way as that for the ground state. Symmetry-wise there are totally 10 components under the $T_{d}$ group: $A_{1},A_{2},E,T_{1}$ and $T_{2}$ IRs (here we do not include SOC in this mixture, as we intend to account for inter-donor and SOC effects separately). By checking the character tables of $D_{2d}$, $C_{3v}(D_{3d})$ or $C_{2v}$ groups for the symmetry operations with $\hat{\mathbf{d}}\|$ [001], $[111]n(i)$ or [110] respectively, one can obtain the allowed mixture components by comparison. To not get distracted from the development of the central physical idea, we list all the results systematically in Appendix A, providing the explicit operation matrices for clarity. For instance, $\psi_{\alpha\uparrow}$ in Eqs. (4)-(6) for $\hat{\mathbf{d}}\|[001]$ follows,

[TABLE]

where ‘$\psi$’ is omitted for shortness, and both the nonpolar ($\gamma=\alpha$) and polar ($\beta$) mixtures are included with time-reversal (TR) compatible phases and small real coefficients $\delta^{\gamma}_{i}$’s. Each $\delta$ is distinct for different states except in the same IR (e.g., $T_{\pm}$ of [001] in Table 1). Without crystal anisotropy, its difference between $S$ and $T$ is due to the exchange part (as opposed to the direct Coulomb), and $\delta^{\beta}_{0}=0$ for $T$ states due to Pauli exclusion. $\delta^{\gamma}_{i}$’s can be obtained perturbatively due to the $H_{\rm int-d}$ interaction (see Appendix B). The magnitude of a general $\gamma_{X}$ mixture into the unperturbed $\alpha_{A_{1}}$ state is on the order of exchange interaction between two-electron states, $\alpha_{A_{1}}(\mathbf{r}_{1})\beta_{A_{1}}(\mathbf{r}_{2})$ and $\alpha_{A_{1}}(\mathbf{r}_{2})\gamma_{X}(\mathbf{r}_{1})$, divided by the energy difference of these two states.

The perturbation theory for relaxation matrix elements throughout Eqs. (12)-(19) then proceeds in a straightforward manner following the above prescriptions. We find the symmetries of the relevant basis states and interaction operators, and determine various perturbation integrals, as shown in the technical details in Appendix C. Here we illustrate the key common aspects by analyzing the representative Eq. (12) in more detail. Both terms require inter-donor interaction amounting to two overlap factors, so neither of them may be neglected a priori. As shown in Eq. (20), SOC is necessary. We find that it is dominated by donor SOC ($\sim\lambda_{\rm soc}\mathbf{L}\cdot\mathbf{s}$, where $\mathbf{L}$ denotes the operator conjugating to the spin in SOC, i.e., $\mathbf{L}\propto\bm{\nabla}V\times\mathbf{p}$). The donor $L_{z}$ and $\epsilon_{xx}-\epsilon_{yy}$ bring $\Delta_{1}$ band back to itself allowing the remaining inter-donor coupling (we note again that the second electron of $\mathbf{r}_{2}$ in the $H_{\rm int-d}$ integral is always in the donor ground states). $\epsilon_{xx}-\epsilon_{yy}$, moreover, connects two available valley configurations comprising $\Delta_{1}$-$1s$ states. Together, the $\gamma=\alpha$ term in Eq. (12) contains a perturbation expansion,

[TABLE]

plus another one with reversed ordering of interactions, where $\overline{H}_{\rm int-d}$ couples $A_{1}$ and $T_{2z}$ envelopes in the same donor [as expected from Eq. (22)] of the $S$ instead of $T_{0}$ state. That leaves only the exchange part of $H_{\rm int-d}$ effective ($\sim J_{A_{1}T_{2}}$). $\mathcal{E}_{\Delta_{1},1s,A_{1}}-\mathcal{E}_{\Delta_{1},1s,E/T_{2}}(\mathcal{E}_{\Delta_{1},3d_{\pm 1},T_{2}})\approx-12(40)$ meV Aggarwal_PR65 , and $\langle E^{I}_{z}|\epsilon_{xx}\!-\!\epsilon_{yy}|A_{1}\rangle=\sqrt{2/3}\Xi_{u}$ with $\Xi_{u}\approx 8.77$ eV Yu_Cardona_Book . Two major donor SOC couplings emerge: one is between two $1s$ configurations and relates to the donor spin splitting $\Delta^{\rm dnr}_{\rm SOC}$($\sim$ 0.03, 0.1, 0.3 meV for P, As, Sb donors respectively; this impurity core effect is deduced from experiments and goes beyond effective mass approximation) Castner_PR67 ; Song_PRL14 . The other is between $1s$ and $3d_{\pm 1}$ where $\langle 3d_{\pm 1}|L_{z}|1s\rangle\neq 0$ within a single ($x$ or $y$) valley due to anisotropy of the $1s$ envelope Kohn_PR55 . We find the $\nu=1s$ component can be safely used for an order of magnitude estimate footnote_3d . For the $\gamma=\beta$ term in Eq. (12), the $J_{A_{1}T_{2}}$ factor is replaced by an exchange term with the polar state $\beta_{A_{1}}\beta_{T_{2}}$ multiplied by another overlap factor which, therefore, has the same order of magnitude. By re-ordering the interactions, we see that two other perturbation expansions of similar or smaller magnitudes are allowed [shown in Eqs. (36) and (37)].

More perturbation terms of similar magnitudes exist, representing the combined interaction operators. The typical example is the Yafet term in the E-Y spin flip mechanism (purely opposite spin states coupled by the SOC part of e-ph interaction) Yafet_SSP63 . Others include phonon-modulated exchange interaction. It is not practical to enumerate each term, however. Even in the much simpler pure bulk Si, numerous comparable leading-order terms contribute Song_PRB12 . In the single-donor relaxation, for example, the Yafet term is not considered Hasegawa_PR60_Roth_PR60_Roth_LincolnLab60 . More importantly, it is also not essential to do so since these terms possess the same symmetry dependence on operator components (e.g., $\epsilon_{ij},s_{i}$) according to the method of invariants Luttinger_PR56 ; Bir_Pikus_Book , justifying this treatment. The net numerical prefactors are difficult to evaluate exactly, and should be left to be extracted experimentally aided by our transparent expressions.

III Quantitative analysis and discussion

Following the derivation shown above, we obtain the leading-order magnitude $|M|$ for relaxation channels $\langle T_{0}|...|S\rangle_{[110]}$, $\langle T_{\pm}|...|S\rangle_{[111]n}$, and $\langle T_{+}+T_{-}|...|S\rangle_{[110]}$ (indexed by $\kappa=2,3,4$ whereas $\langle T_{0}|...|S\rangle_{[001]}$ by $\kappa=1$) via donor SOC,

[TABLE]

where $\lambda$ is the phonon mode, $\Delta\mathcal{E}_{1(2)}=-12$ meV, channel-dependent factor $\mathcal{F}_{1,2,3,4}=1,\frac{1}{\sqrt{2}},\frac{3\sqrt{2}i}{8},\frac{2}{3\sqrt{2}}$, respectively, and $\epsilon^{ph}_{\kappa,\lambda}=\tilde{\epsilon}_{\kappa,\lambda}\sqrt{\hbar(n_{q_{\lambda}}+1)/2\rho\omega_{\lambda}(q)}$. $\tilde{\epsilon}_{\kappa,\lambda}$ for $\kappa=1,2$ is $\frac{1}{2}(\tilde{\epsilon}_{xx}-\tilde{\epsilon}_{yy})$, for $\kappa=3$ is $\frac{1}{3\sqrt{2}}(e^{i\frac{\pi}{6}}\tilde{\epsilon}_{xx}-i\tilde{\epsilon}_{yy}-e^{-i\frac{\pi}{6}}\tilde{\epsilon}_{zz}+e^{-i\frac{\pi}{6}}\tilde{\epsilon}_{xy}+i\tilde{\epsilon}_{xz}-e^{i\frac{\pi}{6}}\tilde{\epsilon}_{yz})$ and for $\kappa=4$ is $\tilde{\epsilon}_{zz}+\frac{1}{2}(\tilde{\epsilon}_{xx}+\tilde{\epsilon}_{yy})$ where $\tilde{\epsilon}_{ij}\equiv i(q_{i}\xi_{j}+q_{j}\xi_{i})/2$, and then phonon polarization $\bm{\xi}(\mathbf{q})$ is projected into $\lambda=$ LA, TA1, TA2 by elastic continuum approximation Ehrenreich_PR56 . $\rho=2.33$ g/cm3, $\hbar\omega_{\lambda}=\hbar v_{\lambda}q$ equals the relaxation energy, $v_{\rm LA(TA)}=8.7(5)\times 10^{5}$ cm/s, and $n_{q}$ is the phonon distribution. Equation (24) is not expected to be exact as discussed above. Our goal is to provide the leading-order estimate for $M$ and its functional dependence on $J$ and temperature. We find for $\langle T_{\pm}|...|S\rangle_{[001]}$ and $\langle T_{+}-T_{-}|...|S\rangle_{[110]}$ the leading-order term vanishes by symmetry and effective-mass analysis, as shown in Appendix C, and our mechanism is overshadowed by hyperfine-induced relaxation Borhani_PRB10 . $M$ between $T$ states are also obtained with a critical difference regarding the inter-donor coupling. Unlike between $S$ and $T$ states, its effect stems from the anisotropic SOC part of $H_{\rm int-d}$ and scales down roughly by a large factor $|J/I|$ where $I=\mathcal{E}_{T_{0}}-\mathcal{E}_{T_{\pm}}$.

Our result differs from the existing hyperfine-coupling theory on two-donor relaxation in Refs. Sugihara_JPCS68 and Borhani_PRB10 in a critical way that the spin mixing scales with $|$SOC coupling/spinless coupling$|$ ($\sim$ 0.03 meV/12 meV) instead of $|$hyperfine($A$)/exchange($J$)$|$ ($A\sim 0.2\mu$eV and $J$ is on the order of 1 meV for $d\sim 8$ nm). Thus our mechanism provides a comparable or stronger relaxation channel at the higher end of the two-qubit coupling regime ($\leq 10$ nm) Kane_Nat98 ; Koiller_PRL02 , and is therefore extremely relevant for the quantum computing platforms involving $S$-$T$ donor-pair qubits Klymenko_JPCM14 ; Dehollain_PRL14 ; Gonzales-Zalba_NanoLett14 . Also our stronger scaling on $J$ (by $J^{2}$) offers a clear experimental distinction between the two mechanisms.

Figure 1(a) quantifies the relaxation time $\tau$ for exchange coupling between 0.1 and 2 meV, by $\tau^{-1}=\frac{2\pi}{\hbar}\int\frac{d^{3}q}{(2\pi)^{3}}\sum_{\lambda}|M_{\lambda}|^{2}\delta(J-\hbar v_{\lambda}q_{\lambda})$ via Eq. (24) and $\int d\Omega_{\mathbf{q}}|\tilde{\epsilon}_{\kappa,\lambda}|^{2}$ footnote_tau . $\tau\propto J^{-5}$ ($k_{B}T\ll J$) or $TJ^{4}$ ($k_{B}T\gg J$). The dependence on the phonon mode, which is rigorous from Eqs. (12)-(19), is made explicit here. The numerical factor due to inter-donor interaction, however, is crudely averaged as exchange splitting $J$. Various exchange terms from different perturbation expansions [as discussed between Eqs. (23) and (24)] lead to superposed oscillations over donor vector $\mathbf{d}$, unrealistic for a perturbation calculation like ours to specify. To have a rough idea about the relative trends between our mechanism and the hyperfine-induced one, we also add the hyperfine result from Ref. Borhani_PRB10 in Fig. 1(a). Although the spin mixing is generally stronger in our mechanism, the phonon-mode related factors $\tilde{\epsilon}_{\kappa,\lambda}$ reduce our relaxation rates somewhat relatively. The ratio of relaxation time over exchange gate time goes as $J\tau/\pi\sim 1/J^{4}$ and decreases sharply with $J$. However, even the smallest ratio in this figure, for $T_{0}\rightarrow S,\hat{\mathbf{d}}\|[001]$ at $J=2$ meV, is about 7 million, well satisfying the need for quantum error correction. We remark that for As and Sb donors, the crossover between SOC and hyperfine mechanisms occurs at smaller $J$. In this regime, our mechanism indicates that quantum computation using As Gonzales-Zalba_NanoLett14 ; Lo_APL14 or Sb Bradbury_06PRL donors has a much shorter intrinsic relaxation time compared to that of P donors. Note that while some uncertainty in the overall prefactors may raise or lower our curves in Fig. 1(a), the parametrically dominant $J^{5}$ dependence is robust and independent of quantitative details, clearly establishing the dominance of our mechanism for stronger inter-donor exchange coupling.

Finally, TR forbids spin relaxation by SOC on a single donor strictly at zero magnetic ($B$) field Hasegawa_PR60_Roth_PR60_Roth_LincolnLab60 whereas it does not affect any transition among the lowest two-electron states. As a result, the application of $B$ field has the most effect on driving one-electron spin flips, which dominate electron relaxation when $g\mu_{B}B\gtrsim J$. Indeed, we see a crossover behavior with increasing $B$ in Fig. 1(b) for an $S\rightarrow T_{-}$ relaxation. The crossover is much sooner for $T_{0}\rightarrow T_{-}$ as $g\mu_{B}B\sim I$. In fact, since the $B$ power-law dependence for single donors and our $J$ or $I$ dependence are exactly the same for any temperature, the crossover occurs at a fixed $g\mu_{B}B/J(I)$ ratio depending only on the donor type. This is an important finding independent of quantitative details.

IV Summary

In conclusion, we have established that two-electron relaxation in Si donor pairs driven by the intrinsic SOC is important, despite the decade-long perception that the SOC in this system is too weak and only hyperfine-induced relaxation needs to be taken into account Sugihara_JPCS68 ; Borhani_PRB10 . It is crucial that we identify that the SOC from the central cell region of the donor core, not from the ion Coulomb potential or bulk Si, is the dominant source. The central cell SOC stems from the fast-varying core potential (which the smooth Coulomb potential does not exhibit) and at the same time breaks the Elliott-Yafet cancellation that comes with the bulk spin-phonon interaction in Si. We unravel the intricate relaxation matrix elements which invoke inter-donor, SOC and e-ph interaction on top of the Si crystal and ionized donor potentials, by utilizing the exact symmetries possessed by the two-donor system Hamiltonian and also the approximate ones by the individual donors and bulk crystal (resulting in various $\Delta$ bands at six points of the $\Delta$ stars). In this way, we are able to finally obtain a relatively simple and robust dependence of the relaxation rates on the donor alignment and exchange interactions, phonon modes and deformation potentials, and SOC constants, which cannot be ascertained by brute-force numerical simulations. Our mechanism dominates the qubit relaxation when the two-donor separation is relatively small or the inter-qubit coupling is large, and is especially important for donor-based singlet-triplet qubits. Therefore, this study is essential for any serious quantum computing proposals based on Si-donor spin qubits.

Together with the hyperfine-driven two-electron relaxation in the coupled donor systems that has been the only known mechanism for half a century Sugihara_JPCS68 ; Borhani_PRB10 , we now have a full theory for the noise-free intrinsic relaxation limit of the two-donor spin states. When both hypefine and SOC are included, the effect should be of higher order and thus comparably negligible (i.e. our effect and hyperfine effect are additive to leading order). Future developments incorporating charge noise, realistic interface or quantum well confinement potentials can expand on our results and provide quantitative estimates for comparisons with experimental data when they become available in the strong donor-coupling regime. On the experimental side, the SOC and hyperfine related relaxation mechanisms can be disentangled and compared by varying donor separation and alignment direction, and by going over different donor species (P, As, Sb).

ACKNOWLEDGMENTS

This work is supported by Laboratory for Physical Sciences and Microsoft.

Appendix A symmetry identifications of the inter-donor coupled mixture

In this appendix, we present the allowed IRs of mixture for each of the donor-alignment directions $\hat{\mathbf{d}}$ to be [001], [111], and [110]. The leading-order mixtures of interest are such that one of the two electron states is the same as that in the unperturbed Heitler-London ground states. The mixed components and their coefficients have to obey all the spatial and time reversal (TR) symmetry of the two-donor systems.

To avoid ambiguity and for consistency, here we provide explicitly the matrices we use (i.e., our convention) for all symmetry operations upon the symmetrized states. See Tables 2-4. The physical conclusions will not depend on the specific convention (and the resulting coefficient phases) we choose as long as it is consistently carried out. We find the following most general symmetry-allowed mixtures without SOC, compatible with the forms of Eqs. (4)-(6). For $\hat{\mathbf{d}}\|[001]$ we find,

[TABLE]

the first of which is Eq. (22) of the main text. The spins are along the $z$ direction, and the small coefficients $\delta_{i}$’s (derived later in Appendix B) are real, and superscripts $c$ and $i$ denote covalence and ionic, respectively. The $\delta$’s have no relation between each other for different singlet and triplet functions unless they are connected by symmetry in a degenerate IR for a given alignment (i.e., $T_{\pm}$ in [001] or [111] alignment). Specifically, we note that $\delta^{i}_{0}=0$ for $T_{0}$ or $T_{\pm}$ due to Pauli exclusion. Regarding the restriction imposed by the symmetry, we note in particular the $C_{2z}$ symmetry operation which keeps the atom site and spin orientation. As a result, each mixed component should transform back to itself under $C_{2z}$. Another constraint is that both $\sigma_{[1\pm 10]}$ and TR connect the same pair of wave functions (i.e., $\alpha\uparrow\leftrightarrow\alpha\downarrow$, $\beta\uparrow\leftrightarrow\beta\downarrow$) up to a global phase difference. This constraint sets the phase of the coefficient for each mixed term, as well as the linear combination of degenerate mixed components.

For $\hat{\mathbf{d}}\|[111]$, we find that the first restriction (by simultaneously satisfying $C^{+}_{3[111]}$ and $C^{-}_{3[111]}$) eliminates the possibility of mixing in $E$ IR and two of the three combinations belonging to $T_{1}$ or $T_{2}$ IR.

[TABLE]

where spins are along the [111] direction, and the $\delta$’s have no relations between different $\hat{\mathbf{d}}$ alignments. Each $\delta_{i}$ for $\alpha$ and $\beta$ needs not be to the same for the $[111]n$ case as no symmetry operation in $C_{3v}$ connects atoms $\alpha$ and $\beta$, while it is the same for $[111]i$. $\delta^{i}_{0}=0$ for $T_{0}$ or $T_{\pm}$.

Finally, for $\hat{\mathbf{d}}\|[110]$, all IR components are allowed to be mixed in, since there is no ‘self-projecting’ operation as in the $\hat{\mathbf{d}}\|$ [001] or [111] case. The $\sigma_{[1-10]}$ and TR symmetry still dictates the linear combination of mutually connected components ($T_{1x}$ and $T_{1y}$; $T_{2x}$ and $T_{2y}$).

[TABLE]

where spins are along the [110] direction. $\delta^{i}_{0}=0$ for $T_{0}$ or $T_{\pm}$.

Appendix B obtaining mixture coefficients $\delta$’s

The mixture coefficients $\delta_{i}$ can be obtained perturbatively. The general procedure for the non-polar mixture is shown in the following. Take the component to be mixed as $X$ IR, which can be any of those in Eqs. (A)-(A). Together with the unperturbed (denoted by superscript ‘0’) basis states

[TABLE]

where $\chi_{A_{1}}=\langle\alpha_{A_{1}}|\beta_{A_{1}}\rangle$, the generic $S(T_{0})$ component with $A_{1}$ and $X$ electrons on different atoms has the form

[TABLE]

where $a,b,c,d$ are determined by symmetry as shown in Eqs. (A)-(A), $\chi_{A_{1}X}=\langle\alpha_{A_{1}}|\beta_{X}\rangle$. These coefficients will cancel out due to the same symmetry operations of the two-donor system when one computes $\langle S(T_{0})^{\prime}|H_{\rm int-d}|S(T_{0})^{0}\rangle$ where $H_{\rm int-d}$ is always identity under all symmetry operations. The mixtures for $T_{\pm}$ states are exactly the same as those for the $T_{0}$ state if one neglects the small anisotropic SOC effect, which accounts for the energy splitting between the triplet states (the nondegeneracy is shown in Table 1 of the main text by symmetry). Between basis states that are not exactly orthogonal, like in our cases due to the finite overlap, the leading-order approximation Herring_book66 is to replace $H_{\rm int-d}$ with $H_{ee}-\mathcal{E}_{S(T_{0})^{0}}$, where $\mathcal{E}_{S(T_{0})^{0}}=\langle S(T_{0})^{0}|H_{ee}|S(T_{0})^{0}\rangle$ and $H_{ee}$ is the total Hamiltonian for the two-electron system [Eq. (1)]. The mixture coefficient of $X$ component follows $\langle S(T_{0})^{\prime}|H_{ee}-\mathcal{E}_{S(T_{0})^{0}}|S(T_{0})^{0}\rangle/(\mathcal{E}_{S(T_{0})}-\mathcal{E}_{S(T_{0})^{\prime}})$, and straightforwardly we have

[TABLE]

where $\mathbf{r}_{1,2}$ is omitted and the first (second) wavefunction in the state belongs to the first (second) electron, and the dielectric constant $\epsilon$ is also omitted above. In Eq. (30), the second to fourth terms are the direct Coulomb energy of two neutral atoms and sum to $\sim e^{-2r_{\alpha\beta}/a^{L}_{B}}$ Heitler_London_ZP27 , where unlike the hydrogen atoms the unperturbed $S(T_{0})^{0}$ of Si donors is formed by the $1s$ state whose decay length are $a^{L(S)}_{B}$ along the larger(smaller) effective mass direction. In Eq. (31), again the third to seventh terms sum up to $\sim e^{-2r_{\alpha\beta}/a^{X}_{B}}$ where $a^{X}_{B}$ is about the longest decay length in the $X$’s envelope. In Eq. (33), the first two terms are the coupling of $\beta_{A_{1}}$ and $\beta_{X}$ by a neutral $\alpha$ atom and largely cancel out. Later we will see that the mixture coefficients due to the direct Coulomb interaction cancel out when calculating the relaxation matrix elements, leaving the physical effects induced only by exchange interaction. Quantitatively, the net effect of the inter-donor perturbation is to couple excited donor states with weights the order of $J(\alpha_{A_{1}}\beta_{A_{1}};\alpha_{A_{1}}\beta_{X})/(\mathcal{E}_{A_{1}}-\mathcal{E}_{X})$. The numerator is similar to the $S-T$ splitting involving the given envelopes in $A_{1}$ (i.e., $1s$) and $X$ (similar arguments in $J$ estimation in Mariya_PR60 ). Furthermore, for relaxation between triplet states, the mixture coefficients cancel out except for the effect of anisotropic SOC not shown above. We can also examine the mixture with polar excited states following a similar procedure. Again, the important part is the difference between mixture coefficients in the singlet and triplet states.

Appendix C Details on perturbative expansion of Eqs. (12)-(19)

Using the physical principles discussed in the main text for identifying the strongest coupling strengths among all perturbations, we obtain the essential perturbation expansions of the relaxation matrix elements derived in Eqs. (12)-(19) one by one.

We start with Eq. (12) of $\langle T_{0}|\epsilon_{xx}-\epsilon_{yy}|S\rangle_{[001]}$. As discussed, inter-donor interaction, $H_{\rm int-d}$, or overlap has to be involved to have finite relaxation. In the first term, there are two ways to have perturbation made by $H_{\rm int-d}$: the mixture of non-polar (covalent) excited states, or the mixture of the polar (ionic) excited states. For the non-polar ones, the one-electron mixture is on the same donor [see Eqs. (A)-(A)] and the coefficient is proportional to $e^{-ir_{\alpha\beta}(1/a^{L}_{B}+1/a^{X}_{B})}$. The polar ones require another overlap factor. In the second term of Eq. (12), one factor utilizes the polar mixture with coefficient $\propto e^{-r_{\alpha\beta}/a^{X}_{B}}$ while the other factor goes with overlap $\propto e^{-r_{\alpha\beta}/a^{L}_{B}}$, and is comparable to the first term (the double overlap contribution mostly cancels out due to TR symmetry).

Next we study the first term due to its non-polar coupling by $H_{\rm int-d}$. We find out that the dominant term comes with donor impurity SOC $\sim L_{z}\sigma_{z}$. $L_{z}$ belongs to the $T_{1z}$ IR of the $T_{d}$ group ($L_{z}\sim T_{1z}$). It only couples the same $\Delta$ band strongly. Then the $H_{\rm ep}$ has to couple the same $\Delta$ band (in addition to the same envelope state) as a result. We find this is possible, although $\epsilon_{xx}-\epsilon_{yy}\sim\Delta^{z}_{2}$ in the $z$ valleys, $\epsilon_{xx(yy)}\sim\Delta^{x(y)}_{1}$, and $\epsilon_{yy(xx)}=\frac{1}{2}[(\epsilon_{yy(xx)}+\epsilon_{zz})+(\epsilon_{yy(xx)}-\epsilon_{zz})]\sim\Delta^{x(y)}_{1}+\Delta^{x(y)}_{2}$ in the $x(y)$ valleys. $\epsilon_{xx}-\epsilon_{yy}$ belongs to the $E^{I}_{z}$ IR of the $T_{d}$ group, to which one of the donor state configuration out of the $\Delta_{1}$ band also belongs. So $\epsilon_{xx}-\epsilon_{yy}$ couples $\Delta_{1}-1s-A_{1}$ to $\Delta_{1}-1s-E^{I}_{z}$, which in turn is coupled by $L_{z}$ to $\Delta_{1}-1s-T_{2z}$ ($E^{I}_{z}\times T_{1z}=T_{2z}$), which is finally coupled back to $\Delta_{1}-1s-A_{1}$ by $H_{\rm int-d}$ [allowed in Eq. (A)]. For the polar mixture contribution, we need to keep everything except replacing the $T_{2z}$ mixture in the same donor to all the possible mixtures $X$’s in the other donor by $H_{\rm int-d}$ each multiplying an overlap $\langle\alpha_{T_{2z}}|\beta_{X}\rangle$. For the second term of Eq. (12), we can replace the inter-donor mixture by a simple overlap, multiplied by a polar mixture with the $T_{2z}$ IR.

As the combined interaction of $H_{\rm ep}$ with diagonal strain elements (which give rise to deformation potential within the $\Delta_{1}$ band) and SOC is used frequently, we show their combined $T_{d}$ symmetry explicitly. $T_{1}\times E=T_{1}+T_{2}$,

[TABLE]

The most important one of the perturbation expressions follows as Eq. (23) of the main text, plus the term reversing perturbation ordering. From Eqs. (30)-(32), we know that the two coefficients are differed by parts due to the exchange interaction integrals. By reversing the matrix element ordering in the numerator, it essentially takes the TR operation. The key is that $L_{i}$ is odd under TR, while $\epsilon_{ij}$ and $H_{\rm int-d}$ are even under TR. These parities yield a minus sign between these two orderings, and as a result the direct Coulomb integrals of the inter-donor mixture in $S$ and $T_{0}$ cancel out. This cancellation makes physical sense, as the singlet and triplets respond the same way to the direct Coulomb part of the interaction. The otherwise cancellation between two same $S$ (or $T$) states can also be generally reached by TR,

[TABLE]

for any $i,j,k$ components.

The other two orderings involve $T_{1z}$ states which cannot be formed out of $s$ envelopes but only with $d_{\pm 1}$ envelopes,

[TABLE]

This term is similar in magnitude to those in Eq. (23) involving $3d$ intermediate states.

[TABLE]

This term is smaller than Eq. (36) because of the two relatively large energy denominators instead of one.

For Eq. (13) of $\langle T_{0}|\epsilon_{z^{\prime\prime}y^{\prime\prime}}|S\rangle_{[110]}$, the first and third terms dominate by one less SOC. As for the $\hat{\mathbf{d}}\|$[001] case, $L_{z^{\prime\prime}}=\frac{1}{\sqrt{2}}(L_{x}+L_{y})$ is required in conjunction with $\epsilon_{z^{\prime\prime}y^{\prime\prime}}$ to couple $\langle\psi^{T_{0}}_{\alpha}|...|\psi^{S}_{\alpha(\beta)}\rangle$, due to the $\sigma_{y^{\prime\prime}}$ reflection symmetry. $\epsilon_{z^{\prime\prime}y^{\prime\prime}}\equiv\epsilon_{[110][1-10]}=\epsilon_{xx}-\epsilon_{yy}$. And we have a coupling similar to that for Eq. (12). Using Eq. (C), the combined interaction is symmetrized as,

[TABLE]

With the $T^{x}_{2}+T^{y}_{2}$ part, we can follow Eqs. (15), (36), and (37) and replace $(\epsilon_{xx}-\epsilon_{yy})L_{z}$ with $(\epsilon_{yy}-\epsilon_{zz})L_{x}$ and $(\epsilon_{zz}-\epsilon_{xx})L_{y}$. With the $T^{x}_{1}+T^{y}_{1}$ part, the perturbation expansion is similar except that there is at least one denominator energy $\sim 40$ meV, since $T_{1}$ cannot be formed from $1s$ envelopes.

For Eq. (14) of $\langle T_{+}|\epsilon_{xz}-i\epsilon_{yz}|S\rangle_{[001]}$, since $\epsilon_{xz}-i\epsilon_{yz}$ does not partially belong to the $\Delta_{1}$ in any valley, it has to go with interband coupling, which cannot be brought back to the $\Delta_{1}$ band state considerably by donor SOC or $H_{\rm int-d}$. As noted in the main text, the host SOC leads to an E-Y cancellation. In all cases, the contributing terms are substantially smaller than Eq. (23).

For Eq. (15) of $\langle T_{+}|\epsilon_{x^{\prime}z^{\prime}}+i\epsilon_{y^{\prime}z^{\prime}}|S\rangle_{[111]n}$, the first (third) term cannot be connected to the second (fourth) term due to the lack of $\alpha\leftrightarrow\beta$ symmetry operations (they are constructive in the $\hat{\mathbf{d}}\|$[001] case and destructive in the $\hat{\mathbf{d}}\|[111]i$ case). Now we go into more details about them. $\langle\uparrow|L_{x^{\prime}}\sigma_{x}+L_{y^{\prime}}\sigma_{y}|\downarrow\rangle=L_{x^{\prime}}-iL_{y^{\prime}}=\sqrt{\frac{2}{3}}[e^{i\pi/6}L_{x}-iL_{y}-e^{-i\pi/6}L_{z}]$ is necessary. $\epsilon_{x^{\prime}z^{\prime}}+i\epsilon_{y^{\prime}z^{\prime}}=\sqrt{2}/3[e^{-i\pi/6}\epsilon_{xx}+i\epsilon_{yy}-e^{i\pi/6}\epsilon_{zz}+e^{i\pi/6}\epsilon_{xy}-i\epsilon_{xz}-e^{-i\pi/6}\epsilon_{yz}]$. Since the e-ph part contains both the diagonal and off-diagonal strain tensor elements, we have respectively the perturbation within the $\Delta_{1}$ band (using diagonal $\epsilon_{ii}$ with donor SOC) and via the $\Delta_{5}$ band (using $\epsilon_{ij},i\neq j$ and host SOC), the latter of which is dropped due to E-Y cancellation.

We study for the diagonal e-ph part. The coefficients in $e^{-i\pi/6}\epsilon_{xx}+i\epsilon_{yy}-e^{i\pi/6}\epsilon_{zz}$ add up to 0, showing there is no $\epsilon_{xx}+\epsilon_{yy}+\epsilon_{zz}\sim A_{1}$ component. Therefore it contains only $E$ IR components. We know $L_{i}\sim T_{1i}$. $E\times T_{1}=T_{1}+T_{2}$. So in the first case the combined e-ph and SOC mixes in some components of $T_{1}$ or $T_{2}$. Now we check Eq. (A) to see whether and what $T_{1,2}$ components can be coupled by $H_{\rm int-d}$. Indeed, not all the components but only $T_{1x}+T_{1y}+T_{1z}$ and $T_{2x}+T_{2y}+T_{2z}$ are allowed. In order to see whether they are contained in the combined e-ph and SOC, we use Eq. (C) to make symmetrization,

[TABLE]

and find out that it does contain both $\sum_{i}T_{1i}$ and $\sum_{i}T_{2i}$ components. Basically, this is a very similar situation as shown in Eq. (C) of the [110] $S-T_{0}$ case. The major contributing terms in the first (or second) term are identified as follows,

[TABLE]

where $i,j,k$ ordered cyclicly. Polar mixture also exist in all four terms of Eq. (15) similarly to that of Eq. (12).

For Eq. (17) of $\langle T_{+}+T_{-}|\epsilon_{A_{1}}|S\rangle_{[110]}$, the first term consists of one spin-independent matrix element, $\langle\psi^{T}_{\alpha\uparrow}|\epsilon_{A_{1}}|\psi^{S}_{\alpha\uparrow}\rangle\approx\langle\psi_{A_{1}}|\epsilon_{A_{1}}|\psi_{A_{1}}\rangle=\Xi_{d}+\Xi_{u}/3$ (which are the usual dilation and shear deformation potentials Yu_Cardona_Book ), and one spin-flip overlap $\langle\psi^{T_{2}}_{\beta\uparrow^{\prime\prime}}|\psi^{S}_{\beta\downarrow^{\prime\prime}}\rangle$, which we treat similarly to the manner above:

[TABLE]

relatively small due to the $3d$ intermediate state. The second term contains $\langle\psi^{T}_{\alpha\uparrow^{\prime\prime}}|\epsilon_{A_{1}}|\psi^{S}_{\alpha\downarrow^{\prime\prime}}\rangle$, where $\epsilon_{A_{1}}$ stands for $\epsilon_{x^{\prime\prime}x^{\prime\prime}}\equiv\epsilon_{zz}$ and $\epsilon_{z^{\prime\prime}z^{\prime\prime}(y^{\prime\prime}y^{\prime\prime})}\equiv\epsilon_{[1\pm 10][1\pm 10]}=(\epsilon_{xx}+\epsilon_{yy}\pm 2\epsilon_{xy})/2$. We can similarly symmetrize the operator as before. Seeking the dominant $T_{2}$ components, we take from $\epsilon_{zz}$ or $\epsilon_{xx}+\epsilon_{yy}$ the $E^{I\!I}_{z}$ part (dropping the $A_{1}$ part), $\sim\epsilon_{xx}+\epsilon_{yy}-2\epsilon_{zz}$, and the $L_{x}-L_{y}$ part from $L_{x^{\prime\prime}}-iL_{y^{\prime\prime}}$ to obtain,

[TABLE]

similarly to Eq. (C), which the rest of the steps follow.

For Eq. (17) of $\langle T_{+}-T_{-}|\epsilon_{x^{\prime\prime}y^{\prime\prime}}|S\rangle_{[110]}$, SOC is needed such that $\langle\psi^{T_{2}}_{\alpha(\beta)}|\epsilon_{x^{\prime\prime}y^{\prime\prime}}L_{z^{\prime\prime}}|\psi^{S}_{\alpha}\rangle$ is nonvanishing (due to $\sigma_{y^{\prime\prime}}$), so the first and third terms require SOC2, much smaller than the second and fourth terms. For $\langle\psi^{T_{2}}_{\alpha\uparrow^{\prime\prime}}|\epsilon_{x^{\prime\prime}y^{\prime\prime}}|\psi^{S}_{\alpha\downarrow^{\prime\prime}}\rangle$, $\epsilon_{x^{\prime\prime}y^{\prime\prime}}=\frac{1}{\sqrt{2}}(\epsilon_{xz}-\epsilon_{yz})$. We have to go through the $\Delta_{5}$ intermediate state and host SOC. The result is suppressed by the E-Y cancellation.

The relaxation matrix elements between different triplet ($T$) states. As discussed in the main text [also see Eq. (35)], the relaxation does not vanish only when the mixtures into the two states by inter-donor interaction are not the same. This can be satisfied between $S$ and $T$ states by exchange interaction even the crystal field is completely uniform. Between the $T$ states, by the same reason, the different mixtures have to come from field anisotropy with respect to triplet spin orientations. The same origin leads to the small splitting ($I$, much smaller than the $S$-$T$ splitting, $J$) between different triplet states. We briefly discuss the perturbation expansion from these matrix elements.

For Eq. (14) of $\langle T_{+}|\epsilon_{xz}+i\epsilon_{yz}|T_{0}\rangle_{[001]}$, due to the non-diagonal e-ph element, $\epsilon_{xz}+i\epsilon_{yz}$, the host SOC is invoked and E-Y cancellation is switched on [as for Eq. (14) of $\langle T_{+}|\epsilon_{xz}-i\epsilon_{yz}|S\rangle_{[001]}$] making this matrix element small even among ones for the triplet relaxations.

For Eq. (15) of $\langle T_{+}|\epsilon_{x^{\prime}z^{\prime}}+i\epsilon_{y^{\prime}z^{\prime}}|T_{0}\rangle_{[111]n}$, the e-ph and SOC couplings follow those of $\langle T_{+}|\epsilon_{x^{\prime}z^{\prime}}+i\epsilon_{y^{\prime}z^{\prime}}|S\rangle_{[111]n}$ in Eq. (15), while the exchange part is much suppressed as discussed above. The first (third) and second (fourth) terms become equal to each other in Eq. (16) of $\langle T_{+}|\epsilon_{x^{\prime}z^{\prime}}\!+\!i\epsilon_{y^{\prime}z^{\prime}}|T_{0}\rangle_{[111]i}$.

In Eq. (18) of $\langle T_{+}+T_{-}|\epsilon_{y^{\prime\prime}z^{\prime\prime}}|T_{0}\rangle_{[110]}$, the first and third terms scale with SOC2 (without taking the inter-donor SOC effect into account) and are much smaller than the second and fourth terms. $\epsilon_{y^{\prime\prime}z^{\prime\prime}}=(\epsilon_{xx}-\epsilon_{yy})/\sqrt{2}\sim E^{I}_{z}$, and $L_{x^{\prime\prime}}-iL_{y^{\prime\prime}}=L_{z}-\frac{i}{\sqrt{2}}(L_{x}-L_{y})\sim T_{1}$. $E\times T_{1}=T_{1}+T_{2}$. So the perturbation expansion has a similar form to that in Eq. (40).

Finally, in Eq. (18) of $\langle T_{+}-T_{-}|\epsilon_{x^{\prime\prime}z^{\prime\prime}}|T_{0}\rangle_{[110]}$ and (19) of $\langle T_{+}-T_{-}|\epsilon_{x^{\prime\prime}y^{\prime\prime}}|T_{+}+T_{-}\rangle_{[110]}$, $\epsilon_{x^{\prime\prime}z^{\prime\prime}}=\frac{1}{\sqrt{2}}(\epsilon_{xz}+\epsilon_{yz})$ and $\epsilon_{x^{\prime\prime}y^{\prime\prime}}=\frac{1}{\sqrt{2}}(\epsilon_{xz}-\epsilon_{yz})$ both invoke E-Y cancellation and suppress the two matrix elements even further, as in Eq. (14) of $\langle T_{+}|\epsilon_{xz}+i\epsilon_{yz}|T_{0}\rangle_{[001]}$.