Spin-relaxation time in the impurity band of wurtzite semiconductors

TL;DR

This paper investigates the spin-relaxation time of electrons in the impurity band of wurtzite semiconductors, combining theoretical models and calculations to match experimental data and predict effects of external tuning.

Contribution

It introduces a detailed theoretical framework for calculating spin-relaxation times in wurtzite semiconductors, incorporating symmetry considerations and spin-orbit interactions.

Findings

Theoretical spin-relaxation times agree with experiments in GaN and ZnO.

Predicted significant increase in relaxation time for InN and AlN under external potential tuning.

Provided estimates for key materials' spin-relaxation times.

Abstract

The spin-relaxation time for electrons in the impurity band of semiconductors with wurtzite crystal structure is determined. The effective Dresselhaus spin-orbit interaction Hamiltonian is taken as the source of the spin relaxation at low temperature and for doping densities corresponding to the metallic side of the metal-insulator transition. The spin-flip hopping matrix elements between impurity states are calculated and used to set up a tight-binding Hamiltonian that incorporates the symmetries of wurtzite semiconductors. The spin-relaxation time is obtained from a semiclassical model of spin diffusion, as well as from a microscopic self-consistent diagrammatic theory of spin and charge diffusion in doped semiconductors. Estimates are provided for particularly important materials. The theoretical spin-relaxation times compare favorably with the corresponding low-temperature measurements in GaN and ZnO. For InN and AlN, we predict that tuning of the spin-orbit coupling constant induced by an external potential leads to a potentially dramatic increase of the spin-relaxation time related to the mechanism under study.

Tables

Table 1. Table 1: Material parameters and resulting spin-relaxation times. The effective masses m ∗ superscript 𝑚 m^{*} and the dielectric constants ϵ italic-ϵ \epsilon are taken from Ref.  Hanada, 2009 , while the spin-orbit interaction coefficients α 𝛼 \alpha , γ 𝛾 \gamma , and b 𝑏 b are extracted from the quoted references. The spin lifetime τ s subscript 𝜏 𝑠 \tau_{s} and the maximum spin lifetime τ s max superscript subscript 𝜏 𝑠 max \tau_{s}^{\mathrm{max}} are calculated within the loop-corrected self-consistent approximation (LCSCA) of Sec.  IV.2 at the density that corresponds to the Mott criterion n i 1 / 3 ​ a ≈ 0.25 superscript subscript 𝑛 i 1 3 𝑎 0.25 n_{\mathrm{i}}^{1/3}a\approx 0.25 for the metal-insulator transition. For the former, the tabulated value of α 𝛼 \alpha is used, while for the latter the optimal value α opt superscript 𝛼 opt \alpha^{\mathrm{opt}} arising from the LCSCA is adopted.

Material	$m^{\ast }$	$ϵ$	$a$ [$\AA$]	$V_0$ [meV]	$\alpha$ [meV$\AA$]	$\gamma$ [meV${\AA }^3$]	$b$	${\tau }_s$ [ns]	${\alpha }^{\mathrm{opt}}$ [meV$\AA$]	${\tau }_s^{\max }$ [ns]
GaN	0.19	10.1	29	50	4.5Stefanowicz et al. (2014)	400Yin et al. (2010)	3.959Fu and Wu (2008)	5	0.2	150
ZnO	0.22	7.8	19	98	2.2Andrearczyk et al. (2005)	320Fu and Wu (2008)	3.855Fu and Wu (2008)	15	0.35	38
InN	0.10	15.0	77	13	13.1Wang et al. (2010)	354Wang et al. (2010)	4.885Wang et al. (2010)	1	0.02	$1.8\times {10}^4$
AlN	0.28	8.5	16	104	-0.72Wang et al. (2010)	6.445Wang et al. (2010)	3.767Wang et al. (2010)	118	0.01	$4.0\times {10}^4$

Equations

$$H$$

$$H_0$$

$$H_{\text{D}}$$

$$H_{\text{extr}}$$

$$H_{\text{D},\text{1}}$$

$$H_{\text{3}}$$

$$\alpha ={\alpha }_{\mathrm{D}}+{\alpha }_{\mathrm{R}}.$$

$$V_{\text{imp}}(\mathbf{r})=\sum _m{V}_m(\mathbf{r})=-\sum _m\frac{e^2}{ϵ|\mathbf{r}-{\mathbf{R}}_m|},$$

$$\mathcal{H}=\sum _{m^{\prime }\ne m}\sum _{{\sigma }^{\prime },\sigma }|m^{\prime }{\sigma }^{\prime }⟩{\mathcal{V}}^{{\sigma }^{\prime },\sigma }({\mathbf{R}}_{m^{\prime }m})⟨m\sigma |,$$

$$\mathcal{V}(\mathbf{r})=\left(\begin{array}[]{cc}\mathcal{V}_{0}(\mathbf{r})+i\mathcal{C}_{z}(\mathbf{r})&i\mathcal{C}_{x}(\mathbf{r})+\mathcal{C}_{y}(\mathbf{r})\\ i\mathcal{C}_{x}(\mathbf{r})-\mathcal{C}_{y}(\mathbf{r})&\mathcal{V}_{0}(\mathbf{r})-i\mathcal{C}_{z}(\mathbf{r})\end{array}\right).$$

$${\mathcal{V}}_0(\mathbf{r})=-V_0\left(1+\frac{r}{a}\right)e^{-r/a}.$$

$$⟨m^{\prime }{\sigma }^{\prime }|H_1|m\sigma ⟩$$

$$⟨m^{\prime }{\sigma }^{\prime }|H_1|m\sigma ⟩=\frac{\sigma \alpha }{a}{\delta }_{{\sigma }^{\prime },\overline{\sigma }}\int \mathrm{d}\mathbf{r}{\varphi }_{m^{\prime }}^{\ast }(\mathbf{r})\frac{(x-X_m)+i\sigma (y-Y_m)}{|\mathbf{r}-{\mathbf{R}}_m|}{\varphi }_m(\mathbf{r}).$$

$$⟨m^{\prime }{\sigma }^{\prime }|H_1|m\sigma ⟩=\frac{\sigma \alpha }{3a^2}{\delta }_{{\sigma }^{\prime }\overline{\sigma }}{R}_{m^{\prime }m}\sin \theta e^{i\sigma \varphi }\left(1+\frac{R_{m^{\prime }m}}{a}\right)e^{-R_{m^{\prime }m}/a},$$

$$⟨m^{\prime }{\sigma }^{\prime }|H_3|m\sigma ⟩=\gamma ⟨m^{\prime }{\sigma }^{\prime }|(b{k}_z^2-k_x^2-k_y^2)(k_y{\sigma }_x-k_x{\sigma }_y)|m\sigma ⟩$$

$$⟨m^{\prime }{\sigma }^{\prime }|H_3|m\sigma ⟩=-\frac{\sigma \gamma }{3a^3}{\delta }_{{\sigma }^{\prime },\overline{\sigma }}\frac{R_{m^{\prime }m}}{a}\sin \theta e^{i\sigma \varphi }{e}^{-R_{m^{\prime }m}/a}\left\{5-b+\left[(b+1){\cos }^2\theta -1\right]\frac{R_{m^{\prime }m}}{a}\right\}.$$

$${\mathcal{C}}_x^{(1)}(\mathbf{r})$$

$${\mathcal{C}}_y^{(1)}(\mathbf{r})$$

$${\mathcal{C}}_x^{(3)}(\mathbf{r})$$

$${\mathcal{C}}_y^{(3)}(\mathbf{r})$$

$$\rho (\vartheta ,t)=\sum _{n=0}^{\infty }\frac{2n+1}{4\pi }\exp \left[-\frac{1}{4}n(n+1)V(t)\right]P_n(\cos \vartheta ).$$

$$\mathbf{S}(t)\cdot {\mathbf{S}}_0=2\pi {\int }_0^{\pi }\mathrm{d}\vartheta \sin \vartheta \cos \vartheta \rho (\vartheta ,t).$$

$$\mathbf{S}(t)\cdot {\mathbf{S}}_0=\exp \left(-\frac{t}{{\tau }_s}\right)$$

$$\frac{1}{{\tau }_{\mathrm{s}}}=\frac{1}{2}\frac{{\chi }^2}{{\tau }_{\mathrm{c}}}$$

$${\mathbf{S}}_{\mathrm{i}}=⟨I|\mathbf{\sigma }|I⟩.$$

$$|F⟩=\mathcal{V}({\mathbf{R}}_{m^{\prime }m})|I⟩,$$

$${\mathbf{S}}_{\mathrm{f}}=\frac{⟨F|\mathbf{\sigma }|F⟩}{⟨F|F⟩},$$

$$\cos {\alpha }_{m^{\prime }m}={\mathbf{S}}_{\mathrm{f}}\cdot {\mathbf{S}}_{\mathrm{i}}.$$

$${\alpha }_{m^{\prime }m}^2\approx \frac{4}{{\mathcal{V}}_0^2({\mathbf{R}}_{m^{\prime }m})}\left\{{\mathcal{C}}_x^2({\mathbf{R}}_{m^{\prime }m})+{\mathcal{C}}_y^2({\mathbf{R}}_{m^{\prime }m})+{\mathcal{C}}_z^2({\mathbf{R}}_{m^{\prime }m})-{\left[{\mathcal{C}}_x({\mathbf{R}}_{m^{\prime }m})S_{\mathrm{i},\mathrm{x}}+{\mathcal{C}}_y({\mathbf{R}}_{m^{\prime }m})S_{\mathrm{i},\mathrm{y}}+{\mathcal{C}}_z({\mathbf{R}}_{m^{\prime }m})S_{\mathrm{i},\mathrm{z}}\right]}^2\right\},$$

$$⟨{\alpha }_{m^{\prime }m}^2⟩\approx \frac{8}{3{\mathcal{V}}_0^2({\mathbf{R}}_{m^{\prime }m})}\left[{\mathcal{C}}_x^2({\mathbf{R}}_{m^{\prime }m})+{\mathcal{C}}_y^2({\mathbf{R}}_{m^{\prime }m})+{\mathcal{C}}_z^2({\mathbf{R}}_{m^{\prime }m})\right].$$

$$P_{m^{\prime }m}=\frac{{\mathcal{V}}_0^2({\mathbf{R}}_{m^{\prime }m})}{{\sum }_{m^{\prime \prime }\ne m}{\mathcal{V}}_0^2({\mathbf{R}}_{m^{\prime \prime }m})}.$$

$$\overline{⟨{\alpha }^2⟩}=\sum _{m^{\prime }\ne m}{P}_{m^{\prime }m}⟨{\alpha }_{m^{\prime },m}^2⟩=\frac{8}{3}\frac{{\sum }_{m^{\prime }\ne m}\left[{\mathcal{C}}_x^2({\mathbf{R}}_{m^{\prime }m})+{\mathcal{C}}_y^2({\mathbf{R}}_{m^{\prime },m})+{\mathcal{C}}_z^2({\mathbf{R}}_{m^{\prime },m})\right]}{{\sum }_{m^{\prime }\ne m}{\mathcal{V}}_0^2({\mathbf{R}}_{m^{\prime }m})}.$$

$${\chi }^2=\frac{8}{3}\frac{\int \mathrm{d}\mathbf{r}{n}_{\mathrm{i}}\left[{\mathcal{C}}_x^2(\mathbf{r})+{\mathcal{C}}_y^2(\mathbf{r})+{\mathcal{C}}_z^2(\mathbf{r})\right]}{\int \mathrm{d}\mathbf{r}{n}_{\mathrm{i}}{\mathcal{V}}_0^2(\mathbf{r})}.$$

Full text

I Introduction

The group-III nitride and the group-II oxide semiconductors have direct bandgaps, which cover the ultraviolet to infrared energy range.Hanada (2009) In particular, the wide bandgaps of AlN, GaN, and ZnO make these materials extremely important for optoelectronic and high-power applications. In the case of ZnO, the present interest is also fueled by the large exciton binding energy that allows for high efficiency light emission up to high temperatures, by the prediction that manganese doping induces room-temperature ferromagnetism,Dietl et al. (2000); Pearton et al. (2003) as well as by the relatively simple crystal-growth technology involved.Özgür et al. (2005); Janotti and de Walle (2009) A key characteristic of this class of materials, linked with very special physical properties, is that their stable crystalline structure is of hexagonal wurtzite (WZ) type. The light nuclei result in a small valence band splitting of otherwise degenerate spin states, leading to a relatively weak spin-orbit coupling, based on which long spin coherence times have been anticipated. Krishnamurthy et al. (2003) A long spin lifetime is a necessary ingredient for spintronics and quantum technology applications,Soumyanarayanan et al. (2016); Awschalom et al. (2013); Bader and Parkin (2010); Fabian et al. (2007) paving the way to all-semiconductor spintronic devices based on this group of materials.Awschalom et al. (2002); Žutić et al. (2004)

Despite the interesting physical properties and the technological promise of the group-III nitride and the group-II oxide semiconductors (and their alloys) with WZ structure, the corresponding spin properties have been less extensively studied, at the experimental and theoretical levels, than in the case of III-V cubic zincblende (ZB) semiconductors. In the latter context, detailed measurements have been performed for n-doped GaAs, yielding very long spin lifetimes (of the order of hundreds of nanoseconds) at low temperatures.Kikkawa and Awschalom (1998); Dzhioev et al. (2002); Oestreich et al. (2005) In particular, density-dependent measurements yielded the longest spin-relaxation times for doping concentrations in the neighborhood of the critical one for the metal-insulator transition (MIT). This intriguing behavior has contributed to the sustained theoretical interest that ZB materials have received.Wu et al. (2010) In the regime of longest spin-relaxation times, for doping densities between the critical one and that where the impurity band hybridizes with the conduction band, a theoretical description based in mechanisms relating spin relaxation to momentum scattering, like Dyakonov-Perel and Elliot-Yafet, cannot be applied from a conceptual point of view. This regime has been addressed, at low temperature and in the absence of magnetic field, in terms of a tight-binding Hamiltonian Tamborenea et al. (2007) which incorporates spin-orbit coupling into the Matsubara-Toyozawa model.Matsubara and Toyozawa (1961) The inclusion of the bulk Dresselhaus spin-orbit interaction resulted in a good agreement with the observed spin lifetimes for several ZB materials. Intronati et al. (2012); Wellens and Jalabert (2016)

The existing low-temperature spin-relaxation experiments in WZ semiconductors point to a scenario that is similar to the ZB case. For Si-doped GaN, the dependence of the relaxation time on both magnetic field and temperature was found to be qualitatively similar to previous studies in n-type GaAs, suggesting a common origin for spin relaxation in these systems.Beschoten et al. (2001) Indeed, at $T=5\,\mathrm{K}$, out of three studied samples with different doping densities, the one having a doping density in the vicinity of the MIT yielded the longest relaxation time (of about 20 ns). These values of the spin-relaxation time were confirmed in the framework of a detailed study of the properties of the MIT in GaN.Wolos et al. (2011a) The spin-relaxation time was found to exhibit a temperature-dependent maximum as a function of doping density.Buß et al. (2010, 2011) In the temperature interval of these studies ($T=80\,\mathrm{K}$ to room temperature) the electronic conductance is dominated by conduction-band properties and therefore the spin relaxation is consistent with the Dyakonov-Perel mechanism.

Spin-relaxation times of up to 20 ns have been measured in doped ZnO at 30 K. Bulk and epilayer samples of different doping densities were investigated and the longest relaxation time was found for the bulk sample.Ghosh et al. (2005) The temperature dependence of the relaxation rate was found to be consistent with the Dyakonov-Perel mechanism.Prestgard et al. (2015) It was noticed that a static in-plane electric field can enhance the spin lifetime in a quantum well geometry.Ghosh et al. (2008)

Aside from Refs. [Beschoten et al., 2001], [Wolos et al., 2011a], and [Ghosh et al., 2005], where samples with different impurity concentrations have been measured, there is no systematic study of the dependence on doping density of the low-temperature spin lifetime close to the metal-insulator transition for the WZ materials. Most experimental data are not for the low-temperature regime, which is where the impurity-band physics is dominant and where the longest lifetimes are expected.

From the theoretical point of view, the ZB and the WZ crystal structures have the same tetrahedral nearest-neighbor atomic coordination number and share the lack of inversion symmetry that results in the splitting of spin-degenerate states. The main qualitative difference between them is that the latter presents a uniaxial anisotropy which is absent in the former. At the level of the envelope-function approximation this anisotropy affects the Dresselhaus spin-orbit effective interaction in the conduction band: its cubic-in-$k$ term becomes anisotropic and a linear-in-$k$ term (formally identical to the Rashba interaction in quantum wells) appears. The linear term has an intrinsic contribution that may depend on substrate-dependent strain in epilayers, and an external one that can be controlled by applying electric fields. The resulting spin relaxation emerges from an interplay of the effects of the cubic and the linear spin-orbit couplings, and can then be influenced by tuning the linear coupling strength.Ganichev and Golub (2014); Harmon et al. (2011)

An important quantitative difference between the WZ and the ZB structures is the much smaller valence-band splitting of the former as compared with the latter. Such feature would generically point to a much larger spin-relaxation time for WZ structures than in the ZB case. The existing experiments do not validate this simplistic conclusion, indicating the important influence on the spin-relaxation mechanisms of the uniaxial anisotropy as well as other physical parameters.

Spin relaxation in the conduction band of WZ structures has been theoretically addressed, using the Dyakonov-Perel formalism,Harmon et al. (2009, 2011); Wang et al. (2010) and predicting particularly long spin lifetimes for the case of AlN, for a narrow range of values of the linear coupling. The temperature- and magnetic field dependence of the spin lifetime was investigated in Ref. [Kang, 2015]. The anisotropy of the spin-orbit coupling was found to yield a dependence of the Dyakonov-Perel spin-relaxation rate on the initial orientation,Harmon et al. (2009); Wang et al. (2010); Buß et al. (2009) with a decay twice as fast for the component along the crystal symmetry axis as compared to in-plane spin. Such an asymmetry has also been obtained in the context of ZB quantum wells with particular crystal orientations. Cartoixà et al. (2005) The experimentally observed Ghosh et al. (2005) nonmonotonic behavior of the spin lifetime as a function of temperature in ZnO has been explained invoking the spin exchange between localized and extended states.Harmon et al. (2009) Alternatively, the numerical solution of Bloch equations for n-type ZnO quantum wells yielded a maximum relaxation time as a function of temperature, and also as a function of carrier density for sufficiently low temperatures. Lü and Cheng (2009) These theoretical studies rely on material parameters and spin-orbit coupling constants that have been extracted from numerical band-structure calculations,Fu and Wu (2008); Wang et al. (2007, 2010); Xu and Ching (1993) and therefore introduce a certain level of uncertainty in the theoretical predictions.

It is important to remark that the theoretical methods used so far to study WZ structures are not expected to be applicable to the low-temperature physics on the metallic side of the metal-insulator transition in the impurity band. Matsubara and Toyozawa (1961); Tamborenea et al. (2007) The rather incomplete experimental and theoretical information on spin physics for the impurity band of WZ semiconductors calls for further studies to explore whether long spin-relaxation times are possible in the vicinity of the metal-insulator transition in the impurity band of n-doped semiconductors having WZ crystal structure.

In this article we adapt theoretical methods developed in Refs. [Intronati et al., 2012] and [Wellens and Jalabert, 2016] in the context of ZB semiconductors, in order to calculate the spin-relaxation time in the metallic side of the impurity band of WZ bulk semiconductors, and compare our results with the experimental data available in the literature for GaN and ZnO. We also study how far the spin-relaxation time can in principle be extended by tuning the linear-in-$k$ component of the Dresselhaus spin-orbit interaction. We find that this tuning is promising for GaN and potentially dramatic for AlN.

From a practical point of view, the understanding of spin-relaxation mechanisms in WZ bulk semiconductors serves a two-fold purpose. One one hand, important information can be obtained in the cases where materials having ZB crystal structure in the bulk turn to a WZ configuration when nanostructured in nanorods or quantum dots. De and Pryor (2010); Koguchi et al. (1992); Intronati et al. (2013); Intronati (2013) On the other hand, even if the spin relaxation can be considerably slowed down in semiconductor heterostructures and nanostructures, doped bulk semiconductors are always required for the contact layers of the devices. From the fundamental point of view, it is important to understand how the symmetries of the WZ structure affect the spin-relaxation mechanisms, and also to have available a tool to better characterize the metal-insulator transition in wide bandgap semiconductors.

This article is organized as follows. In Section II we introduce the Hamiltonian and other basic elements of the tight-binding model with Dresselhaus spin-orbit interaction used to describe an electron in the impurity band of bulk wurtzite semiconductors. In particular, we obtain the hopping matrix elements of the Dresselhaus spin-orbit terms of the Hamiltonian; details of this calculation are given in Appendix A. Sections III and IV present, respectively, the semiclassical and the fully quantum-mechanical self-consistent approach to the spin-relaxation time. The latter with three successive degrees of approximation. In Sec. V we apply our theory to four specific semiconductors and examine the proposed scheme for maximization of the spin lifetime. Our conclusions are given in Section VI. For completeness, in Appendix B we apply the semiclassical formula to the case of zincblende materials. Relevant Fourier transforms are given in Appendix C.

II Hamiltonian and hopping amplitude matrix

The envelope-function approximation for conduction-band electrons in WZ semiconductors incorporates the crystal lattice-scale physics into the effective one-body Hamiltonian Nozières and Lewiner (1973); Engel et al. (2007)

[TABLE]

The spin-independent part $H_{0}$ is determined by the effective mass ($m^{*}$) and the electrostatic potential $V(\mathbf{r})$ including all potentials aside from that of the crystal lattice. We note $\mathbf{p}$ the momentum operator, $\mathbf{k}=\mathbf{p}/\hbar$, and $\bm{\sigma}$ the vector of Pauli matrices. The Dresselhaus (intrinsic) term $H_{\text{D}}$, enabled by the bulk inversion asymmetry, comprises in the case of WZ semiconductors, the linear- Lew Yan Voon et al. (1996) and the cubic-in-$k$ Wang et al. (2007); Fu and Wu (2008) components, respectively given by

[TABLE]

and therefore depends on three material-dependent constants, namely, $\alpha_{\mathrm{D}}$, $\gamma$, and $b$. These parameters are in some cases only partially known, as will be discussed in our treatment of specific materials in Sec. V. The extrinsic term $H_{\text{extr}}$ of Eq. (1d) has the same form as the spin-orbit interaction in vacuum, but with a material-dependent renormalized coupling constant $\lambda$. External potentials translating into an electric field in $z$-direction result in a Rashba contribution to $H_{\text{extr}}$ that has the same linear-in-$k$ form of $H_{\text{D},\text{1}}$, with a prefactor $\alpha_{\mathrm{R}}$ that depends on $\lambda$. While such a field cannot be significant in the bulk case, perpendicular electric fields in thin films might play a role. In the latter case, the strain induced on the film by the substrate could be an additional mechanism to modify the linear-in-$k$ coupling. The linear-in-$k$ contributions can be jointly treated by defining a Hamiltonian $H_{1}$ given by (2a), but with a coupling coefficient

[TABLE]

Separating between Dresselhaus and Rashba contributions with the same functional form is generically a challenge, that has been experimentally addressed in the case of ZB quantum wells with the application of an external magnetic field.Meier et al. (2007) Given the uncertainties in the material parameters, we will not attempt such a separation in our theoretical treatment of the WZ structure and we will express our results as a function of the tunable parameter $\alpha$. In the literature, $H_{1}$ is often simply referred to as the ”Rashba” Hamiltonian Fu and Wu (2008); Stefanowicz et al. (2014), while the cubic-in-$k$ contribution, $H_{3}$, is labeled as the ”Dresselhaus” term. We will avoid here this somewhat misleading nomenclature.

The potential arising from the ionized donor impurities is given by

[TABLE]

where $\epsilon$ is the dielectric constant of the semiconductor and $\mathbf{R}_{m}$ represents the impurity positions. Its contribution to $H_{\text{extr}}$ has been shown to be extremely weak in the case of ZB structures, Tamborenea et al. (2007) and we expect the same considerations to also hold in the WZ case. Therefore, it will not be considered in our theoretical analysis. On the other hand, the potential $V_{\text{imp}}(\mathbf{r})$ strongly affects the orbital motion of the electrons, and will be the only contribution to $V(\mathbf{r})$ that we will keep in Eq. (1b).

The potential $V_{m}(\mathbf{r})$ gives rise to hydrogenic states centered at the impurity $m$, with a ground-state wave-function $\phi_{m}(\mathbf{r})=\phi(|\mathbf{r}-\mathbf{R}_{m}|)$. The anisotropy of the WZ crystal lattice induces a uniaxial anisotropy, of about ten percent, in the effective masses and the dielectric constants. Hanada (2009); Persson et al. (2001); Ferreira da Silva et al. (2001) We will neglect the resulting small deformation of the hydrogenic impurity wave functions, Rodina et al. (2001) adopting the standard isotropic ground-state wave-function $\phi(\mathbf{r})=(1/\sqrt{\pi a^{3}})\exp{(-r/a)}$, where $a$ stands for the effective Bohr radius. The electronic ground states of the different sites $m$ provide a restricted basis $\{|m\sigma\rangle\}$ to describe the electron jumping between impurity centers ($\sigma=\pm$ corresponds to a spin projection of $\pm\hbar/2$ in the $z$-direction). Choosing as energy origin the ground-state energy, the Hamiltonian in this restricted space can be expressed as

[TABLE]

where $\mathbf{R}_{m^{\prime}m}=\mathbf{R}_{m^{\prime}}-\mathbf{R}_{m}$. In the following, we neglect the overlaps between different states $m\neq m^{\prime}$ and thereby treat $\{|m\sigma\rangle\}$ as an orthonormal basis, which is justified in the regime that we are interested in (i.e., for moderate doping densities and energies close to the center of the impurity band) Majlis and Anda (1983); Figueira et al. (1984). Like in the ZB case, the hopping amplitudes can be generically expressed through a $2\times 2$ matrix in the spin subspace(Wellens and Jalabert, 2016)

[TABLE]

The spin-independent hopping amplitudes are those of the Matsubara-Toyozawa model Matsubara and Toyozawa (1961)

[TABLE]

We note ${\bf r}=(x,y,z)$ and $r=|{\bf r}|$, while $V_{0}=2E_{\mathrm{Ry}}^{(0)}m^{*}/\epsilon^{2}$ corresponds to twice the binding energy of the impurity sites ($E_{\mathrm{Ry}}^{(0)}=13.6\,\mathrm{meV}$ is the Rydberg energy).

The matrix elements of the linear-in-$k$ spin-orbit Hamiltonian, $H_{1}$, are

[TABLE]

We use the notation $\overline{\sigma}=-\sigma$, and thus $\sigma_{x}|\sigma\rangle=|\overline{\sigma}\rangle$ and $\sigma_{y}|\sigma\rangle=i\sigma|\overline{\sigma}\rangle$. Applying the operators $k_{x}=-i\partial_{x}$ and $k_{y}=-i\partial_{y}$ we obtain

[TABLE]

In Appendix A we perform the analytical integration which results in

[TABLE]

where $\mathbf{R}_{m^{\prime}m}=R_{m^{\prime}m}(\sin\theta\cos\phi\,\mathbf{x}+\sin\theta\sin\phi\,\mathbf{y}+\cos\theta\,\mathbf{z})$.

As shown in Appendix A, the matrix elements of the cubic-in-$k$ term

[TABLE]

can be cast in the form

[TABLE]

According to Eqs. (10) and (12), in the case of the WZ crystal structure we have $\mathcal{C}_{z}(r)=0$, and $\mathcal{C}_{j}(r)=\mathcal{C}^{(1)}_{j}(r)+\mathcal{C}^{(3)}_{j}(r)$ for $j=x,y$, with

[TABLE]

and

[TABLE]

In what follows we calculate the spin-relaxation rate in the above-defined model, using two different theoretical approaches.

III Semiclassical approach to the spin lifetime

The spatial diffusion of electrons through the network of impurities is accompanied by a corresponding dynamics of the electronic spin. The spin-orbit coupling implies that the spin is not conserved in a hopping event between two impurities. In this Section we present the derivation of the spin-relaxation rate based on the concept of spin diffusion in the Bloch sphere that occurs when the electron undergoes diffusive motion in real space. Identifying the electron spin with a continuous vector of fixed norm in three-dimensional space amounts to a semiclassical description of the spin diffusion.

III.1 Diffusion on a sphere

The evolution of the electron spin can be assimilated to a random walk on a sphere with a mean squared rotation angle $\chi^{2}$ in each step of the random walk. From an initial spin orientation $\mathbf{S}_{0}$ (the electron spin in units of $\hbar/2$) corresponding to a distribution of absolute certainty that the point is at a given position on the sphere, the resulting distribution after a random walk of $N(t)$ steps after a time $t$ is given by Roberts and Ursell (1960)

[TABLE]

Here, $\vartheta$ is the angular distance from the initial orientation, $P_{n}$ are the Legendre polynomials, and $V(t)=N(t)\chi^{2}$ is the variance of the corresponding plane motion. For steps being associated with hopping events occurring with a rate $\tau_{\mathrm{c}}^{-1}$, one has $N(t)=t/\tau_{\mathrm{c}}$ and then $V(t)=(t/\tau_{\mathrm{c}})\chi^{2}$.

The mean projection of the resulting spin orientation $\mathbf{S}(t)$ on the initial orientation is given by

[TABLE]

With the expression (15) and using the orthogonality of the Legendre polynomials, one finds an exponential decay of the mean spin projection

[TABLE]

with the spin-relaxation rate

[TABLE]

that depends on the hopping rate and the mean-square spin-rotation angle.

III.2 Spin rotation in a hopping event

We consider an initial state where the electron is localized on the impurity $m$ with a spin state $|I\rangle$ (not necessarily a state $|\sigma\rangle$ oriented along the $z$-axis). The initial spin expectation value writes

[TABLE]

The final spin state $|F\rangle$ after a hop to impurity $m^{\prime}$ is obtained as

[TABLE]

where $\mathcal{V}(\mathbf{R}_{m^{\prime}m})$ is the hopping matrix given by Eq. (6). The final spin expectation value is computed as

[TABLE]

and the spin rotation angle $\alpha_{m^{\prime}m}$ related to the hopping event can be extracted from the scalar product with the initial spin orientation as

[TABLE]

Using (20) and the general form of the matrix elements (6), we calculate the spin rotation. To lowest order in the spin rotation components of $\mathcal{V}(\mathbf{r})$, that we assume to be small as compared to the spin-independent hopping terms, $\mathcal{C}_{j}(\mathbf{r})\ll\mathcal{V}_{0}(\mathbf{r})$, and therefore having small rotation angles with $\cos\alpha\approx 1-\alpha^{2}/2$, we get the result

[TABLE]

which depends on the initial spin orientation $\mathbf{S}_{\mathrm{i}}$. An average over the initial spin orientation based on a homogeneous distribution over the surface of the Bloch sphere yields the mean square rotation angle occurring in a hopping event

[TABLE]

III.3 Impurity averaged rotation angle

In order to use the result of the diffusion on a sphere (18), we need the average square spin rotation angle for hops between any pair of impurities. Towards the evaluation of that average, we start with the probability that a hop from the initial impurity $m$ ends at a given impurity $m^{\prime}$

[TABLE]

We have neglected the small corrections due to the spin-dependent hopping terms $\mathcal{C}_{j}(\mathbf{r})$. Weighting the mean square spin rotation angle (24) with this probability leads to an impurity and spin-orientation averaged squared spin rotation angle

[TABLE]

We identify the squared mean rotation angle $\chi^{2}$ used in the formalism of the diffusion on a sphere with the average of (26) over impurity configurations. This impurity average is obtained through the replacement of the sums over final positions by integrals over the spatial vector between initial and final site, weighted by the constant density of impurities $n_{\mathrm{i}}$, yielding

[TABLE]

Defining the inverse hopping rate from the time $\tau_{\mathrm{c}}$ when the occupation of a site is halved, Intronati et al. (2012); Intronati (2013) we have

[TABLE]

Replacing the denominator of (27) using (28) and plugging the result in the diffusion expression (18), we get the mean spin-relaxation rate

[TABLE]

The spin-independent hopping is dominated by $\mathcal{V}_{0}(\mathbf{r})$ [see definition after Eq. (6)]. Those terms are the hopping elements of the Matsubara-Toyozawa model. Matsubara and Toyozawa (1961) Using this form in (28), one gets the hopping timeIntronati et al. (2012)

[TABLE]

and the general expression of the mean spin-relaxation rate in terms of the model parameters as

[TABLE]

with the general property that the relaxation rate is proportional to the square root of the dopant density.

III.4 Application to materials with wurtzite crystal structure

The expression (31) for the spin-relaxation rate is completely general, only relying on the form (6) of the hopping amplitudes. The dominant spin-dependent terms depend on the crystal structure, and the application of the semiclassical spin-diffusion formalism to zincblende materials is presented in Appendix B.

For wurtzite materials, the spin parts of the hopping matrix elements have significant contributions from both, the linear-in-$k$ and the cubic-in-$k$ terms in the spin-orbit coupling. Using the expressions of Eqs. (13) and (14) in the general formula for the averaged spin-relaxation rate (31) and performing the integral over hopping vectors, one gets with the hopping time (30) the mean spin-relaxation time in wurtzite structures

[TABLE]

where we introduced the material-dependent energy $E_{\rm soc}$ associated to the spin-orbit coupling and given by

[TABLE]

The first term is due to the linear-in-$k$ spin-orbit coupling and the second and third terms appear due to the presence of the cubic-in-$k$ contribution.

As discussed in Sec. II, while the cubic coupling strength $\gamma$ is a material parameter, the linear coupling constant $\alpha$ is composed of an intrinsic component $\alpha_{\mathrm{D}}$ and an extrinsic Rashba component $\alpha_{\mathrm{R}}$. Thus, $\alpha$ can be affected by external influences like strain or an electric field, and its value can be controlled and optimized in the search for maximum spin lifetime. Our general result (32) implies a minimal spin-relaxation rate of

[TABLE]

occurring for the optimal value $\alpha^{\mathrm{opt}}$ given by

[TABLE]

where we have assumed $b=4$ for the evaluation of the approximate numerical values.

Defining the spin lifetime as the inverse of the mean relaxation rate (32), one gets the Lorentzian dependence on the linear coupling strength $\alpha$ depicted in Figure 1, with a maximum spin lifetime $\tau^{\mathrm{max}}$ corresponding to the inverse of the minimal relaxation rate (34). The relaxation time has a pronounced maximum as a function of $\alpha$ with high values close to $\alpha^{\mathrm{opt}}$. The width at half maximum of the Lorentzian dependence of the lifetime on $\alpha$ is given by

[TABLE]

A similar strongly peaked dependence of the spin lifetime on the linear coupling has been obtained Harmon et al. (2011) for conduction-band electrons in wurtzite semiconductors based on the Dyakonov-Perel mechanism. The presentation of Fig. 1 in terms of dimensionless quantities is universal and can be used for different materials once their characteristic constants are determined.

III.5 Anisotropy of the spin-relaxation rate

It can however be noticed that the spin-orbit coupling is anisotropic in wurtzite structures, with $\mathcal{C}_{z}(\mathbf{r})=0$. If we fix the initial spin orientation along the $z$-axis in (23), instead of averaging over the initial spin direction, we get the squared mean rotation angle

[TABLE]

and thus a (longitudinal) spin-relaxation rate

[TABLE]

that is enhanced by a factor of 3/2. In contrast, when the initial spin orientation is in the $x$-$y$ plane, a reduced (transverse) spin-relaxation rate of

[TABLE]

is obtained. Thus, the initial relaxation is twice as fast for a spin oriented along the symmetry axis [001] of the crystal than for a spin perpendicular to that axis. Such an anisotropy of the spin relaxation does not occur in zinc blende structures where the components of the spin-orbit coupling do not have a preferential direction (see Appendix B).

The factor 2 appearing in the ratio between $\tau_{\mathrm{s},x-y}$ and $\tau_{\mathrm{s},z}$ is a general feature of the WZ structure, and it is in line with the anisotropic electron spin relaxation measured Buß et al. (2009) in bulk GaN at a temperature $T=80\,\mathrm{K}$, in the regime where the Dyakonov-Perel mechanism sets the spin lifetimes.

IV Self-consistent approach

The semiclassical approach to the spin lifetime presented in the previous Section has been shown to be extremely successful in the case of ZB materials, Intronati et al. (2012) as it gives good account of the existing measurements and the results of microscopic theories. For the low-temperature spin relaxation in the WZ materials, the experimental situation is somewhat uncertain, and microscopic theories have not been developed. The microscopic approach developed for the case of ZB structures Wellens and Jalabert (2016) uses self-consistent approximations where the locator expansion for the one- and two-particle Green functions (and their irreducible components) fulfill important constraints, like particle conservation. Vollhardt and Wölfle (1980, 1992); Kroha et al. (1990) We present in this Section different schemes of the self-consistent approximation applicable to the WZ crystal structure. While the general features of the diagrammatic perturbation theory are the same as in the ZB case, the reduced symmetry of the WZ structure leads to considerably different results.

The impurity-averaged retarded (advanced) Green function $G^{(\pm)}(\varepsilon)$ can be written in terms of the corresponding self-energy $\Sigma^{(\pm)}(\varepsilon)$ through Dyson equation

[TABLE]

We note $\text{z}_{\pm}=\varepsilon\pm i\eta$, with $\eta$ an infinitesimal positive quantity. $G^{(\pm)}$ and $\Sigma^{(\pm)}$ are understood as on-site, and the impurity average makes the choice of the site irrelevant. Moreover, these matrices are proportional to the $2\times 2$ identity matrix, and this is why the spin indices are not explicit.

The intensity propagator $\Phi$ is a two-Green function object that, according to the Bethe-Salpeter equation, writes as

[TABLE]

in terms of its irreducible component $U$. We note $\varepsilon=(\varepsilon_{1}+\varepsilon_{2})/2$ and $\hbar\omega=\varepsilon_{1}-\varepsilon_{2}$ the semi-sum and the difference of the energies $\varepsilon_{1}$ and $\varepsilon_{2}$ of the two Green functions defining $\Phi$, and $\mathbf{r}$ the vector difference between the positions of the impurities defining the Green functions. The interest of the intensity propagator is that it allows us to obtain the probability distribution governing the spatial and spin diffusion

[TABLE]

where $\rho(\varepsilon)$ denotes the impurity-averaged density of states, obtained as

[TABLE]

From Eqs. (13) and (14) (together with $C_{z}=0$) it follows that the hopping amplitude matrix $\mathcal{V}(\mathbf{r})$ defined in Eq. (6) fulfills ${\mathcal{V}}(x,y,z)=e^{i\frac{\phi}{2}\sigma_{z}}{\mathcal{V}}(x^{\prime},y^{\prime},z)e^{-i\frac{\phi}{2}\sigma_{z}}$ with $x^{\prime}=\cos(\phi)x-\sin(\phi)y$ and $y^{\prime}=\sin(\phi)x+\cos(\phi)y$ and ${\mathcal{V}}(x,y,z)=\sigma_{y}{\mathcal{V}}(x,-y,z)\sigma_{y}$. These transformation properties dictate that the Fourier transform of the matrix $U$ has the following form for ${\bf q}=0$:

[TABLE]

Restricting ourselves to the two-dimensional subspace of diagonal spin density operators (entries $++$ and $--$ of the matrix $\tilde{U}$), we obtain the probability-conserving condition, as well as the longitudinal relaxation rate (in the $z$ direction) Wellens and Jalabert (2016)

[TABLE]

The remaining subspace (entries $+-$ and $-+$) provides the damping of coherences in the chosen representation, corresponding to the transverse relaxation rate (in the $x$ and $y$ directions)

[TABLE]

The access to the energy-dependent spin-relaxation rates provided by the self-consistent approximation is important in order to be able to address not only the case of uncompensated semiconductors (with a half-filled impurity band), but also that of weak compensation. In addition, it is a necessary information for treating the hot-electron condition that might arise from the initial spin injection. Wu et al. (2010)

IV.1 Simplest self-consistent approximation

In the simplest self-consistent approximation (SSCA) the locator expansion for the self-energy is restricted to the only term that represents a processes where the electron hops from site $m$ to another site $m^{\prime\prime}\neq m$, and then back to $m$. Therefore,

[TABLE]

with

[TABLE]

where $E_{\rm soc}$ has been defined in Eq. (33).

The self-consistent retarded self-energy $\Sigma^{(+)}(\varepsilon)$ with negative imaginary part reads

[TABLE]

and, according to Eq. (43), the resulting density of states is given by the semicircle law

[TABLE]

where $\Theta$ stands for the Heaviside door function. Since, within this approximation, the impurity band is symmetric around the energy origin, the Fermi energy is $\varepsilon_{\mathrm{F}}=0$

The irreducible component $U$ of the intensity propagator compatible with the simplest approximation (51) for the self-energy is such that

[TABLE]

and it can thus be expressed as the tensor product $U(\varepsilon,\omega,{\bf r})=n_{\mathrm{i}}{\mathcal{V}}({\bf r})\otimes{\mathcal{V}}^{*}({\bf r})$. In Fourier space we have

[TABLE]

where $\tilde{\mathcal{V}}({\bf k})$ stands for the Fourier transform of the hopping amplitude matrix given in Appendix C, and ${\bf k}_{\pm}={\bf k}\pm{\bf q}/2$.

According to Eqs. (48) and (55),

[TABLE]

From (49), (52), and (54) the longitudinal relaxation rate for electrons at the Fermi energy $\varepsilon_{\mathrm{F}}$ is

[TABLE]

and the transverse relaxation rate is

[TABLE]

The term that multiplies $1/V_{0}^{2}$ in the denominator of Eq. (58) stems from the small correction of the density of states due to the spin-orbit coupling. Neglecting such a term, and taking into account the numerical factors of Eqs. (38)-(39) relating the averaged spin-relaxation rate with the longitudinal and transverse ones, we see that the semiclassical and the simplest self-consistent approaches yield the same functional form for the relaxation rates, with only a difference of $1/\sqrt{2}$ between the prefactors. Such a difference is not surprising, since in the semiclassical approach the exact definition of the relaxation time is to some extent arbitrary. For instance, the hopping time $\tau_{\mathrm{c}}$ is defined in (28) as the time when the occupation of a site is halved while a somewhat different criterion would have been equally adequate. Within the diagrammatic approach, the hopping time can be expressed in terms of the imaginary part of the self energy as follows:Wellens and Jalabert (2016)

[TABLE]

which, when taking into account Eqs. (52) and (53), differs from (28) again by a factor $\sqrt{2}$. In total, the simplest self-consistent approximation exactly reproduces the semiclassical approach with hopping time $\tau_{\mathrm{c}}$ defined according to (60).

IV.2 Loop-corrected self-consistent approximation

The simplest self-consistent approximation developed in Subsec. IV.1 needs to be complemented by adding the terms of the locator expansion that represent processes in which the electron visits more than one impurity before hopping back to the starting one. Such an improvement constitutes the so-called loop-corrected self-consistent approximation (LCSCA), Wellens and Jalabert (2016) which for the spinless case, amounts to the approach used by Matsubara and Toyozawa Matsubara and Toyozawa (1961) in order to obtain the density of states in the impurity band, as well as the conductivity within the diffusion approximation.

Since the self-energy is now expressed as a geometrical series representing hopping events, Eq. (40) can be written as a self-consistent equation for $\Sigma^{(\pm)}(\varepsilon)$

[TABLE]

where we note $\tilde{\mathcal{V}}^{(+)}({\bf k})=\tilde{\mathcal{V}}({\bf k})$, $\tilde{\mathcal{V}}^{(-)}({\bf k})=\tilde{\mathcal{V}}^{*}({\bf k})$. This equation needs to be numerically solved. Since the spin-orbit coupling only gives a small correction to the density of states, we solve for $\Sigma^{(\pm)}(\varepsilon)$ in Eq. (61) trading $\tilde{\mathcal{V}}^{(\pm)}({\bf k})$ by $\tilde{\mathcal{V}}_{0}({\bf k})$.

The loop-corrected sequence translates into a renormalized hopping amplitude matrix

[TABLE]

with $\mathbb{I}_{2}$ the $2\times 2$ unit matrix. Thus, the irreducible component of the intensity propagator takes the form (56), but with the hopping amplitude matrix now replaced by the effective one. Wellens and Jalabert (2016) That is,

[TABLE]

Neglecting the spin-orbit flipping terms in the denominator of (62) results in

[TABLE]

Performing the angular integrals

[TABLE]

and the longitudinal spin-relaxation rate follows from Eq. (49) after the $k$-integration. From (65) we readily see that, similarly to the SSCA, the spin-relaxation time within the LCSCA has a Lorentzian dependence on the variable $\alpha a^{2}/\gamma$. However, contrary to the simpler approximation, the density dependence is not universal in the parameter $n_{i}a^{3}$. As in the SSCA approximation, we have $\tilde{u}_{1}=\tilde{u}_{3}$, and therefore the same relationship (59) between the longitudinal and the transversal relaxation times.

In Fig. 1 we present the LCSCA spin-relaxation times in dimensionless variables for three different impurity densities (solid lines). For low values of $\alpha a^{2}/\gamma$, the relaxation time decreases as the doping density increases. Such a behavior is in line with that of the ZB case in the metallic side of the MIT. Intronati et al. (2012); Wellens and Jalabert (2016) The optimal values of $\alpha a^{2}/\gamma$ that maximize the spin-relaxation time are close to those of the simpler approximations. Interestingly, the linear-in-$k$ term of the WZ structure induces a cross-over value beyond which the relaxation times increase with impurity density.

IV.3 Repeated-scattering-corrected self-consistent approximation

The LCSCA can be improved by the inclusion of cross diagrams in the locator expansion for the self-energy that describe the repeated scattering from selected impurities. The so-called repeated-scattering-corrected self-consistent approximation (RSCSCA) Wellens and Jalabert (2016) restricts the repeated scattering to just a pair of impurities, allowing for arbitrary loops between them [represented by the renormalized hopping amplitude (62)]. The irreducible component of the intensity propagator in the RSCSCA has an expression considerably more complicated than that of Eq. (63) since, even restricting to a pair of repeatedly visited impurities, there is an important proliferation of contributing diagrams.

The RSCSCA results (not shown) are very close to those of the LCSCA, leaving aside the high-energy tail of the impurity band, as well as the very low densities, for which the repeated scattering is more relevant. Elyutin (1981) Since our approach is restricted to impurity densities larger than the critical one for the metal-insulator transition, the latter case does not apply to our study. Since we work with uncompensated or weakly compensated samples, the first case is also not relevant.

It is interesting to remark that, while the quantitative improvements of the RSSCA respect to the LCSCA are generically small in the range of parameters that we work with, the magnitude of these corrections depends on the value of the linear-in-$k$ $\alpha$ coefficient. For small values of $\alpha$, where the spin-relaxation time is close to $\tau_{s}^{max}$, the corrections are of the order of 10%, while for larger values of $\alpha$ the correction practically vanishes. The latter result is a consequence of the fact that for large $\alpha$ the integral leading to $\tilde{u}_{2}$ is dominated by the small $k$-values, which in turn are associated with large values of $r$, where the repeated scattering is not relevant.

V Spin-relaxation times in GaN, ZnO, InN, and AlN

Direct comparison with the few available low-temperature data of spin-relaxation times in WZ materials is hindered by the limited knowledge of some material parameters and the uncertainty on certain experimental conditions, like the excess doping-density beyond the critical one or the tuning of the linear-in-$k$ spin-orbit coupling constant induced by an electrostatic potential. Moreover, the exact nature of the MIT in wide gap WZ doped semiconductors has rarely been experimentally investigated, other than in the case of GaN.Wolos et al. (2011a) For these materials, the standard Mott criterion, as well as more refined calculations, Ferreira da Silva and Persson (2002); Ferreira da Silva et al. (2001) only provide a qualitative estimate of the critical densities.

Despite these limitations, we verify that our theoretical model yields the appropriate order of magnitude of the spin-relaxation times for particularly important materials. We then analyze the improvements that can be made on the spin-relaxation time by adjusting the tunable linear-in-$k$ coupling constant $\alpha$, along the lines of the proposal made in Ref. [Wang et al., 2010] to reduce Dyakonov-Perel relaxation in wurtzite quantum wells. In Table 1 we summarize the used material parameters and present results for the spin-relaxation times obtained within the framework of the LCSCA. The differences with the other calculational schemes are not particularly important, especially in the regime of large $\alpha a^{2}/\gamma$, where the linear-in-$k$ term dominates the spin relaxation.

The quoted effective masses $m^{*}$ and dielectric constants $\epsilon$ are direction-averaged values, i.e. $m^{*}=(m^{\parallel}+2m^{\perp})/3$, with $m^{\parallel}$ and $m^{\perp}$ being, respectively, the longitudinal and transverse effective masses with respect to the c-axis of the WZ structure. Hanada (2009) These parameters determine the effective Bohr radius $a$ and the spin-independent hopping amplitude $V_{0}$, according to $a=a_{0}\epsilon/m^{*}$ and $V_{0}=2E_{\mathrm{R}}^{(0)}m^{*}/\epsilon^{2}$, where $a_{0}$ and $E_{\mathrm{R}}^{(0)}$ are, respectively, the bare values of the Bohr radius and the Rydberg energy of an isotropic hydrogen atom. Rodina et al. (2001)

The values of the spin-orbit interaction coefficients $\alpha$, $\gamma$, and $b$, were taken from the quoted references. Whenever known, the experimentally determined values were used. For instance, in the case of GaN we chose for the cubic-in-$k$ spin-orbit coupling constant the experimental value Yin et al. (2010) $\gamma=400\,\mathrm{meV\AA^{3}}$, instead of the prediction arising from tight-binding band-structure calculations Fu and Wu (2008) $\gamma=330\,\mathrm{meV\AA^{3}}$. Similarly, for the linear-in-$k$ spin-orbit coupling constant we use the value $\alpha=4.5\,\mathrm{meV\AA}$ determined from weak antilocalization measurements, Stefanowicz et al. (2014) rather than the result $\alpha=9\,\mathrm{meV\AA}$ arising from ab initio computations. Majewski (2005); Majewski and Vogl (2005) In the case of ZnO we used the experimental value Andrearczyk et al. (2005) $\alpha=2.2\,\mathrm{meV\AA}$ instead of $\alpha=1.1\,\mathrm{meV\AA}$ given by band-structure calculations. Lew Yan Voon et al. (1996)

The values of $\tau_{s}$ in Table 1 are extracted from the data of Fig. 1, using the material-dependent physical parameters and setting the doping density to the critical one, according to the Mott criterion $n_{\mathrm{i}}^{1/3}a\approx 0.25$. The column $\alpha^{\mathrm{opt}}$ indicates the value of the linear-in-$k$ coupling constant for which the relaxation time is maximum. The last column $\tau_{s}^{\mathrm{max}}$ gives the value of the spin-relaxation time obtained when the linear-in-$k$ coupling constant takes the value $\alpha^{\mathrm{opt}}$, and sets up an upper bound for the times than can be achieved by tuning the contribution arising from $\alpha_{\mathrm{R}}$ by the application of an external voltage in the $z$-axis direction. The electrical tuning of the Rashba spin-orbit interaction can be very large. In particular two-fold Scherübl et al. (2016) and six-fold Liang and Gao (2012) tuning of the Rashba coefficient have been reported in InAs nanowires. In the case of bulk and epilayer wurtzite materials, estimating the extent of the tuning span is not simple and remains a task for future investigations. The values of $\tau_{s}^{\rm max}$ listed in Table 1 therefore indicate potential improvements, realizable under the condition that $\alpha$ can be tuned to the required value $\alpha^{\mathrm{opt}}$. The width $\Delta\alpha$ at half maximum of the Lorentzian dependence of the lifetime, given approximately by Eq. (36), indicates how close $\alpha$ has to be to $\alpha^{\mathrm{opt}}$ in order to obtain a lifetime that is of the order of the optimal value $\tau_{s}^{\rm max}$.

For GaN, the calculated $\tau_{s}$ is of the same order as the measured maximum spin-relaxation time of around 20 ns at low temperature, Wolos et al. (2011a); Beschoten et al. (2001) and the discrepancy between them could be due to the uncertainty in the precise value of $\alpha$. Stefanowicz et al. (2014); Majewski (2005); Majewski and Vogl (2005); Wolos et al. (2011b) In particular, the large number of dislocations reported by Ref. Beschoten et al., 2001 might affect the anisotropy of the crystal and thereby the linear spin-orbit coupling. Furthermore, we see that tuning $\alpha$ from $4.5$ meV$\AA$ to $\alpha^{\rm opt}=0.2$ meV$\AA$ leads to a 30-fold increase of the spin relaxation time related to the mechanism under study.

ZnO is the other material for which low-temperature spin-relaxation experimental data is available. Values of $\tau_{s}$ around $20~{}\mbox{ns}$ at $T=30~{}\text{K}$ have been reported, Ghosh et al. (2005) albeit for a sample with a density of $1.26\times 10^{15}\text{cm}^{-3}$, which is very low compared to the critical MIT density. With this caveat, we remark the good agreement between our theoretical prediction and the measured value of the spin lifetime in ZnO. The comparatively larger theoretical estimation of $\tau_{s}$ in ZnO with respect to GaN is essentially due to the fact that, in the case of ZnO, the value of $\alpha$ is closer to the optimum value $\alpha^{\rm opt}$ already without external field and, in addition, the width $\Delta\alpha$ at half maximum of the Lorentzian dependence of the lifetime, see Eq. (36), is larger ($\Delta\alpha=1.4$ meV$\AA$ for ZnO as compared to $\Delta\alpha=0.8$ meV$\AA$ for GaN). On the other hand, the upper bound $\tau_{s}^{\rm max}$ for ZnO is considerably smaller than for GaN, mainly due to the proportionality of $\tau_{s}^{\rm max}$ to $a^{6}$ (see the universal scaling indicated in the unit of $\tau_{s}$ in Fig. 1).

InN and AlN represent extreme cases in terms of their spin-orbit parameters. While the value of $\alpha_{\mathrm{D}}$ of InN is large, the one corresponding to AlN is very small, and close to $\alpha^{\mathrm{opt}}$. Furthermore, AlN is peculiar in the sense that the theoretical value Wang et al. (2010) for the cubic coupling $\gamma=6.445\,\mathrm{meV\AA^{3}}$ is very small in comparison with the values of other materials discussed. Note that, for both materials, the optimized spin-relaxation time $\tau_{s}^{\mathrm{max}}$ arising from the spin-orbit mechanisms described in Eqs. (1)-(2) is extremely high. In the case of InN this is mainly due to the large value of the Bohr radius $a$, while for AlN it is due to the small magnitude of $\gamma$. Thus, these two materials are promising candidates to obtain long spin-relaxation times, especially AlN which would require only a weak tuning of $\alpha$. The small values of the optimization window $\Delta\alpha$ for InN and AlN, $0.1$ and $0.04$ meV$\AA$ respectively, translate into very narrow and tall peaks for the $\alpha$-dependence of the spin-relaxation time. A similar conclusion was obtained from a theoretical study of spin relaxation in the conduction band at high temperatures, Harmon et al. (2011) where a very high and narrow peak was found in the $\alpha$-dependence of the spin lifetime, with a maximum value of up to 0.5 $\mu$s at room temperature. To the best of our knowledge, no measurements of the low-temperature spin-relaxation time near the metal-insulator transition in the impurity band have been reported for these two interesting semiconductors.

VI Conclusions

We have theoretically studied spin relaxation in the metallic regime of the impurity band in semiconductors with wurtzite crystal structure. We adapted theoretical concepts previously developed and successfully applied in the context of zinc blende semiconductors. Our basic model is solved using two approaches, namely, a semiclassical one and a fully quantum-mechanical microscopic theory. The latter is pursued at three levels of self-consistent approximation: the simplest (SSCA), the loop-corrected (LCSCA), and the repeated-scattering-corrected (RSCSCA), incorporating progressively higher-order terms of the spin-orbit interaction in a diagrammatic locator expansion.

The anisotropic nature of the wurtzite structure gives rise to a corresponding anisotropy in the spin-relaxation time. The interplay between the linear- and cubic-in-$k$ terms of the spin-orbit Hamiltonian leads to a richer scenario of physical behavior, as compared with the case of zinc blende materials which only have cubic terms. The obtained theoretical results are generically expressed in terms of material-dependent parameters and can thus be applied to different wurtzite semiconductors of practical interest. We emphasize that our theory contains no adjustable parameters, and all comparisons with experiment have been done using the most reliable material parameters available in the literature.

The theoretical estimates of the spin-relaxation time obtained with the LCSCA were shown to be of the same order of magnitude as the available experimental data on GaN and ZnO. At the quantitative level, while the theory underestimates somewhat the spin lifetime for GaN found in experiments, it gives a fairly accurate result for ZnO. The discrepancy between theory and experiment for GaN calls for further investigation of this important material. Although the scenario for ZnO looks at this point more consistent than that for GaN, experimental data for ZnO with densities closer to the MIT would be needed to draw definite conclusions. With an eye to potential spintronic applications, we also analyzed the maximization of the spin-relaxation time made possible by adjusting the linear-in-$k$ term of the spin-orbit interaction, and showed that radical improvements could be made for GaN and AlN, leading to lifetimes given by the discussed mechanism that are so long that another mechanism can be expected to dominate the spin relaxation.

Acknowledgements.

We gratefully acknowledge support from the French National Research Agency ANR (Projects ANR-11-LABX-0058_NIE and ANR-14-CE36-0007-01), the French-Argentinian collaborative project PICS 06687, Universidad de Buenos Aires (Project UBACyT 2011-2014 No. 20020100100741), and CONICET (Project PIP 11220110100091).

Appendix A Calculation of hopping matrix elements

In this Appendix we calculate analytically the matrix elements of the spin-orbit Hamiltonian between hydrogenic impurity states. For the linear-in-$k$ term we change the integration variables in Eq. (9) by shifting the origin of coordinates to $\mathbf{R}_{m}$, and we call $\mathbf{R}_{m^{\prime}m}=\mathbf{R}_{m^{\prime}}-\mathbf{R}_{m}$, obtaining

[TABLE]

We have defined

[TABLE]

We write $\mathbf{R}_{m^{\prime}m}=(X_{m^{\prime}m},Y_{m^{\prime}m},Z_{m^{\prime}m})=R_{m^{\prime}m}(\sin\theta\cos\phi\,\mathbf{x}+\sin\theta\sin\phi\,\mathbf{y}+\cos\theta\,\mathbf{z})$, and then Eq. (67a) takes the form

[TABLE]

The rotation of the coordinate system by $(\theta,\phi)$ induces a coordinate change characterized by the transformation

[TABLE]

Performing the corresponding change of coordinates, the integral in (68) becomes

[TABLE]

The cylindrical symmetry around the $z^{\prime}$-axis calls for a further change of variables in favor of the cylindrical coordinates $(\rho,\alpha,z)$:

[TABLE]

The integration over $\alpha$ only leaves the last term, and then

[TABLE]

Analogously we get

[TABLE]

and collecting both terms in Eq. (66) we get

[TABLE]

The integral $I_{1}$ can be solved exactly:

[TABLE]

leading to the matrix element of the linear Dresselhaus spin-orbit coupling of Eq. (10).

The cubic-in-$k$ matrix elements (11) can be written as

[TABLE]

We have six terms of the form $\langle n|k_{i}k_{j}^{2}|m\rangle$, where $i,j=x,y,z$ ($i$ and $j$ may be equal). After acting with the differential operators on the hydrogenic states we get

[TABLE]

The remaining calculation is lengthier than the one of the linear-in-$k$ matrix element, but it follows the same steps: i) a shift of the origin of coordinates to $\mathbf{R}_{m}$, ii) a rotation of the coordinate system, iii) a switch to cylindrical coordinates. Rather than writing all the details of this procedure, we just make explicit the analytical expression of some useful two-dimensional integrals:

[TABLE]

Using the integrals $I^{3,1}(\xi)=\frac{\xi}{3}e^{-\xi}$, $I^{1,3}(\xi)=\left(\frac{\xi}{2}+\frac{\xi^{2}}{3}\right)e^{-\xi}$, $I^{3,0}(\xi)=e^{-\xi}$, $I^{1,2}(\xi)=\frac{1}{2}(1+\xi)\,e^{-\xi}$, $I^{0}(\xi)=e^{-\xi}$, and $I^{1}(\xi)=\frac{\xi}{2}e^{-\xi}$, we obtain for the cubic matrix element (12).

Appendix B Application of the semiclassical approach to zincblende semiconductors

In the case of impurities in semiconductors with zincblende crystal structure, the spin-orbit coupling yields the symmetric form Intronati et al. (2012); Wellens and Jalabert (2016)

[TABLE]

where $\gamma$ is the Dresselhaus spin-orbit coupling strength. Using these expressions in (31) and performing the spatial integral, one gets with (30) the final result for the spin-relaxation time in zincblende structures

[TABLE]

where the value of the numerical prefactor is approximately $0.36$, correcting the one of Eq. (18) in Ref. Intronati et al., 2012.

Appendix C Fourier transform of the hopping amplitude matrix

The diagrammatic expansions are more easily performed in reciprocal space, therefore it is useful to work with the Fourier transform of the hopping amplitude matrix

[TABLE]

The spin-independent part is given by

[TABLE]

In the case of the WZ crystal structure we have $\tilde{\mathcal{C}}_{z}(r)=0$, and thus $\tilde{\mathcal{C}}_{j}(\mathbf{k})=\tilde{\mathcal{C}}^{(1)}_{j}(\mathbf{k})+\tilde{\mathcal{C}}^{(3)}_{j}(\mathbf{k})$ for $j=x,y$. These Fourier transforms can be readily done by a rotation of the integration variables that places the new $z$-axis in the direction of the wave vector $\mathbf{k}$:

[TABLE]

[TABLE]

These analytical expressions allow us to perform the relevant integrals of our work, like those of Eq. (65).