Extrinsic Spin-Orbit Coupling and Spin Relaxation in Phosphorene

TL;DR

This paper derives an effective Hamiltonian for monolayer phosphorene to analyze how extrinsic spin-orbit coupling, influenced by external electric fields and substrate impurities, affects spin relaxation anisotropy.

Contribution

It introduces a detailed derivation of extrinsic spin-orbit coupling in phosphorene considering symmetry and external fields, highlighting anisotropic effects on spin relaxation.

Findings

Spin splitting is linear in electric field and atomic spin-orbit strength.

Anisotropic Fermi contours influence spin relaxation times.

Charged impurity scattering affects spin dynamics.

Abstract

An effective Hamiltonian is derived to describe the conduction band of monolayer black phosphorus (phosphorene) in the presence of spin-orbit coupling and external electric field. Envelope function approximation along with symmetry arguments and Lowdin partitioning are utilized to derive extrinsic spin-orbit coupling. The resulting spin splitting appears in fourth order perturbation terms and is shown to be linear in both the magnitude of the external electric field and the strength of the atomic spin-orbit coupling, similar to the Bychkov-Rashba expression but with an in-plane anisotropy. The anisotropy depends on the coupling between conduction band and other bands both close and distant in energy. The spin relaxation of conduction electrons is then calculated within the Dyakonov-Perel mechanism where momentum scattering randomizes the polarization of a spin ensemble. We show how the anisotropic Fermi contour and the anisotropic extrinsic spin splitting contribute to the anisotropy of spin-relaxation time. Scattering centers in the substrate are considered to be charged impurities with screened Coulomb potential.

Tables

Table 1. Table 1: Character table of the space group of phosphorene ( C ​ m ​ c ​ a 𝐶 𝑚 𝑐 𝑎 Cmca : 64) for the Γ Γ \Gamma point. The plus and minus signs denote the parity under spatial inversion i 𝑖 i .

	$\left\{E\right|0\}$	$\left\{C_{2x}\right|\tau \}$	$\left\{C_{2y}\right|0\}$	$\left\{C_{2z}\right|\tau \}$	$\left\{i\right|0\}$	$\left\{R_x\right|\tau \}$	$\left\{R_y\right|0\}$	$\left\{R_z\right|\tau \}$	Basis	Invariant
${\Gamma }_1^{+}$	1	1	1	1	1	1	1	1		$k_x^2+k_y^2$
${\Gamma }_2^{+}$	1	-1	1	-1	1	-1	1	-1	$xz$	${\sigma }_y$
${\Gamma }_3^{+}$	1	1	-1	-1	1	1	-1	-1	$yz$	${\sigma }_x$
${\Gamma }_4^{+}$	1	-1	-1	1	1	-1	-1	1	$xy$	${\sigma }_z$
${\Gamma }_1^{-}$	1	1	1	1	-1	-1	-1	-1
${\Gamma }_2^{-}$	1	-1	1	-1	-1	1	-1	1	$y$	$k_y$
${\Gamma }_3^{-}$	1	1	-1	-1	-1	-1	1	1	$x$	$k_x$
${\Gamma }_4^{-}$	1	-1	-1	1	-1	1	1	-1	$z$	$k_z=-id/dz$

Table 2. Table 2: First principle calculations setup

Pseudopotential Type	Ultrasoft, Fully Relativistic
Exchange-Correlation Function	PBE
Kinetic Energy Cutoff	34.0 Ry
Charge Density Cutoff	136.0 Ry
Convergence Threshold	${10}^{-6}$ Ry
$k$-Point Grid	Monkhrost $12\times 12\times 1$
# of bands	24
Interlayer Spacing	$31.22=$ Å

Table 3. Table 3: Parameters of the crystal structure of bulk and monolayer black phosphorus based on the definition in Ref. Takao  et al. , 1981 .

	$a$ (Å)	$b$ (Å)	$d_1$ (Å)	$d_2$ (Å)	${\alpha }_1$	${\alpha }_2$
Bulk	4.376	3.314	2.224	2.244	96.34∘	102.09∘
Monolayer	4.552	3.306	2.224	2.262	93.89∘	102.98∘

Table 4. Table 4: Contributions of different bands in the extrinsic spin-orbit coupling quantified by coefficients α x , l ​ l ′ subscript 𝛼 𝑥 𝑙 superscript 𝑙 ′ \alpha_{x,ll^{\prime}} and α y , l ​ l ′ subscript 𝛼 𝑦 𝑙 superscript 𝑙 ′ \alpha_{y,ll^{\prime}} according to Eq. 18 . The values are obtained from first principle calculations and are in atomic units. The last row lists the total contributions of the bands which determines the anisotropy between extrinsic spin-orbit coefficients.

$l$	band#	$l^{\prime }$	band#	${\alpha }_{x,l{l}^{\prime }}$	${\alpha }_{y,l{l}^{\prime }}$
${\Gamma }_1^{+}$	4	${\Gamma }_2^{+}$	10	-71.24
	4		12	2.95
	9		10	10.08
	15		10	-6.39
	17		10	-5.68
${\Gamma }_1^{+}$	17	${\Gamma }_3^{-}$	13	-2.58
	21			-5.75
${\Gamma }_2^{+}$	10	${\Gamma }_3^{-}$	5	32.84
			13	-2.76
			20	46.36
			22	-8.77
${\Gamma }_1^{+}$		${\Gamma }_3^{+}$			1
${\Gamma }_1^{+}$	4	${\Gamma }_2^{-}$	14		2.10
	17				-2.76
	21				-7.02
${\Gamma }_3^{+}$		${\Gamma }_2^{-}$			1
${\sum }_{l,l^{\prime }}$				-10.95	7.67

Equations

$$\mathcal{H}=\underset{H_0}{\underbrace{\frac{p^2}{2m_0}+V_0(r)}}+\underset{H_{\text{SO}}}{\underbrace{\frac{{\hslash }^2}{4m_0^2{c}^2}p\text{dotproduct}\sigma \times \text{gradient}{V}_0(r)}}+V(z)\cdot$$

$${\Psi }_n(r)=e^{i{k}_{\parallel }\text{dotproduct}{r}_{\parallel }}\sum _{l,\sigma }{f}_{nl\sigma }(z)u_{l\sigma 0}(r)|\sigma ⟩,$$

$$\begin{split}&\sum_{l,\sigma}\bigg{[}\quantity(E_{l}({\bf\it 0})+\frac{\hbar^{2}(k_{\parallel}^{2}-\frac{d^{2}}{dz^{2}})}{2m_{0}}+V(z))\delta_{ll^{\prime}}\delta_{\sigma\sigma^{\prime}}+\Delta_{l^{\prime}\sigma^{\prime}l\sigma}\\ &+\frac{\hbar}{m_{0}}({\bf\it k}_{\parallel}-i\frac{d}{dz}\hat{{\bf\it z}})\dotproduct{\bf\it P}_{l^{\prime}\sigma^{\prime}l\sigma}\bigg{]}f_{nl\sigma}(z)=E_{n}f_{nl^{\prime}\sigma^{\prime}}(z),\end{split}$$

$$H_{cc}=\frac{{\hslash }^2{k}_x^2}{2m_x}+\frac{{\hslash }^2{k}_y^2}{2m_y}+{\lambda }_x{k}_x{\sigma }_y+{\lambda }_y{k}_y{\sigma }_x,$$

$$\frac{1}{m_{x(y)}}=\frac{1}{m_0}\text{quantity}(1+\frac{2}{m_0}\sum _l\frac{|P_{x(y),cl}{|}^2}{E_c-E_l}),$$

$$\begin{array}{rl}{\lambda }_{x(y)}=& (\xi \hslash /2)\frac{dV(z)}{dz}\frac{{\hslash }^2}{m_0^2}\\ \times & \sum _{l,l^{\prime }}\Im P_z{L}_{y(x)}{P}_{x(y)}\frac{(E_c-E_l^{\prime })\pm (E_c-E_l)}{(E_c-E_l{)}^2(E_c-E_l^{\prime }{)}^2},\end{array}$$

$$\text{partialderivative}{\rho }_{k}t=-\frac{1}{i\hslash }[{\rho }_k,{\Omega }_k\text{dotproduct}\sigma ]-\sum _{k^{\prime }\ne k}{W}_{k{k}^{\prime }}({\rho }_k-{\rho }_{k^{\prime }}),$$

$$\text{partialderivative}\overline{\rho }t=\frac{1}{i\hslash }\overline{[{\rho }_k^{\prime },{\Omega }_k\text{dotproduct}\sigma ]},$$

$$\frac{1}{i\hslash }[\overline{\rho },{\Omega }_k\text{dotproduct}\sigma ]=\sum _{k^{\prime }\ne k}{W}_{k{k}^{\prime }}({\rho }_k^{\prime }-{\rho }_{k^{\prime }}^{\prime })\cdot$$

$${\rho }_k^{\prime }=\frac{\frac{1}{\hslash }(\hat{s}\times {\Omega }_k)\text{dotproduct}\sigma +{\sum }_{k^{\prime }\ne k}{W}_{k{k}^{\prime }}{\rho }_{k^{\prime }}^{\prime }}{{\sum }_{k^{\prime }\ne k}{W}_{k{k}^{\prime }}}\cdot$$

$$U(q)=\frac{2\pi e^2}{\kappa (q+q_s)}{e}^{-qd},$$

$$\sum _{k^{\prime }\ne k}{W}_{k{k}^{\prime }}\to \frac{n}{2\pi \hslash }\int d{\ell }^{\prime }\frac{\text{absolutevalue}{U\text{quantity}(q)}^2}{\text{absolutevalue}\text{gradient}E\text{quantity}(k^{\prime })},$$

$$\begin{array}{r}{H}_{0,l{l}^{\prime }}=(E_l+\frac{{\hslash }^2(k_{\parallel }^2-(d/dz{)}^2}{2m_0}){\delta }_{l{l}^{\prime }},\\ H_{l{l}^{\prime }}^{\prime }=V(z){\delta }_{l{l}^{\prime }}+{\Delta }_{l{l}^{\prime }}+\frac{\hslash }{m_0}(k_{\parallel }-i\frac{d}{dz})\cdot P_{l{l}^{\prime }s}\cdot \end{array}$$

$$H_{cc}^{(0)}=H_{0,cc},$$

$$H_{cc}^{(1)}=H_{cc}^{\prime },$$

$$H_{cc}^{(2)}=\sum _l\frac{H_{cl}^{\prime }{H}_{lc}^{\prime }}{E_c-E_l},$$

$$H_{cc}^{(3)}=-\frac{1}{2}$$

$$+$$

$$H_{cc}^{(4)}=\frac{1}{2}$$

$$-\frac{1}{2}$$

$$\times \text{quantity}(\frac{(E_c-E_{l^{\prime }})+(E_c-E_l)}{(E_c-E_l{)}^2(E_c-E_{l^{\prime }}{)}^2})$$

$$-\frac{1}{2}$$

$$+$$

$$\begin{split}H_{cc}^{(2)}=\frac{\hbar^{2}}{m_{0}^{2}}\Bigg{[}k_{x}^{2}&\sum_{l=\Gamma_{2}^{+}}\frac{|P_{x,cl}|^{2}}{E_{c}-E_{l}}+k_{y}^{2}\sum_{l=\Gamma_{3}^{+}}\frac{|P_{y,cl}|^{2}}{E_{c}-E_{l}}\\ -\frac{d^{2}}{dz^{2}}&\sum_{l=\Gamma_{1}^{+}}\frac{|P_{z,cl}|^{2}}{E_{c}-E_{l}}\Bigg{]}+\mathcal{O}(\xi^{2})\sigma_{0}\cdot\end{split}$$

$$\begin{split}H^{\prime}_{cl}H^{\prime}_{lc}H^{\prime}_{cc}=-\frac{\hbar^{2}|P_{z,cl}|^{2}}{m_{0}^{2}}\Bigg{(}\frac{d^{2}V(z)}{dz^{2}}+2\frac{dV(z)}{dz}\frac{d}{dz}\\ +V(z)\frac{d^{2}}{dz^{2}}\Bigg{)}\cdot\end{split}$$

$$\approx \sum _{l,l^{\prime \prime }}\frac{H_{cl}^{\prime }{H}_{ll}^{\prime }{H}_{l{l}^{\prime \prime }}^{\prime }{H}_{l^{\prime \prime }c}^{\prime }}{(E_c-E_l{)}^2(E_c-E_{l^{\prime \prime }})}+\sum _{l,l^{\prime }}\frac{H_{cl}^{\prime }{H}_{l{l}^{\prime }}^{\prime }{H}_{l^{\prime }{l}^{\prime }}^{\prime }{H}_{l^{\prime }c}^{\prime }}{(E_c-E_l)(E_c-E_{l^{\prime }}{)}^2}\cdot$$

$$\begin{split}\expectationvalue{H_{cc}^{(4)}}&\approx(\xi\hbar/2)\frac{dV(z)}{dz}\frac{\hbar^{2}}{m_{0}^{2}}\Bigg{[}\\ k_{x}\sigma_{y}&\sum_{\begin{subarray}{c}l=\Gamma_{1}^{+}\\ l^{\prime}=\Gamma_{2}^{+}\end{subarray}}\Im{P_{z,cl}L_{y,ll^{\prime}}P_{x,l^{\prime}c}}\frac{(E_{c}-E_{l^{\prime}})+(E_{c}-E_{l})}{(E_{c}-E_{l})^{2}(E_{c}-E_{l^{\prime}})^{2}}\\ +k_{x}\sigma_{y}&\sum_{\begin{subarray}{c}l=\Gamma_{1}^{+}\\ l^{\prime}=\Gamma_{3}^{-}\end{subarray}}\Im{P_{z,cl}P_{x,ll^{\prime}}L_{y,l^{\prime}c}}\frac{(E_{c}-E_{l^{\prime}})+(E_{c}-E_{l})}{(E_{c}-E_{l})^{2}(E_{c}-E_{l^{\prime}})^{2}}\\ -k_{x}\sigma_{y}&\sum_{\begin{subarray}{c}l=\Gamma_{2}^{+}\\ l^{\prime}=\Gamma_{3}^{-}\end{subarray}}\Im{P_{x,cl}P_{z,ll^{\prime}}L_{y,l^{\prime}c}}\frac{(E_{c}-E_{l^{\prime}})-(E_{c}-E_{l})}{(E_{c}-E_{l})^{2}(E_{c}-E_{l^{\prime}})^{2}}\\ +k_{y}\sigma_{x}&\sum_{\begin{subarray}{c}l=\Gamma_{1}^{+}\\ l^{\prime}=\Gamma_{3}^{+}\end{subarray}}\Im{P_{z,cl}L_{x,ll^{\prime}}P_{y,l^{\prime}c}}\frac{(E_{c}-E_{l^{\prime}})+(E_{c}-E_{l})}{(E_{c}-E_{l})^{2}(E_{c}-E_{l^{\prime}})^{2}}\\ +k_{y}\sigma_{x}&\sum_{\begin{subarray}{c}l=\Gamma_{1}^{+}\\ l^{\prime}=\Gamma_{2}^{-}\end{subarray}}\Im{P_{z,cl}P_{y,ll^{\prime}}L_{x,l^{\prime}c}}\frac{(E_{c}-E_{l^{\prime}})+(E_{c}-E_{l})}{(E_{c}-E_{l})^{2}(E_{c}-E_{l^{\prime}})^{2}}\\ -k_{y}\sigma_{x}&\sum_{\begin{subarray}{c}l=\Gamma_{3}^{+}\\ l^{\prime}=\Gamma_{2}^{-}\end{subarray}}\Im{P_{y,cl}P_{z,ll^{\prime}}L_{x,l^{\prime}c}}\frac{(E_{c}-E_{l^{\prime}})-(E_{c}-E_{l})}{(E_{c}-E_{l})^{2}(E_{c}-E_{l^{\prime}})^{2}}\Bigg{]}\\ &=(\xi\hbar/2)\frac{dV(z)}{dz}\frac{\hbar^{2}}{m_{0}^{2}}\sum_{l,l^{\prime}}\alpha_{x,ll^{\prime}}k_{x}\sigma_{y}+\alpha_{y,ll^{\prime}}k_{y}\sigma_{x}\cdot\end{split}$$

Full text

††preprint: APS/123-QED

I introduction

Extrinsic spin-orbit coupling induced by an external electric field in two-dimensional electron systems lifts the spin degeneracy while it preserves time-reversal symmetry. The induced spin splitting, which is proportional to the magnitude of the field and the crystal wavevector, enables control of spin through movement of charge and vice versa. This effect, which has enabled several phenomena and ideas in spintronics and beyond Manchon et al. (2015), was originally derived by Ohkawa and Uemura Ohkawa and Uemura (1974) for an inversion layer of zinc-blende crystals. Later, Vasko Vas’ko (1979), and Bychkov and Rashba Bychkov and Rashba (1984), generalized the spin splitting for a two-dimensional electron system with an isotropic in-plane effective mass. Unlike Ohkawa and Uemura’s derivation based on the Kane’s model of zinc-blende crystals Kane (1957), the Vasko and Bychkov-Rashba spin splittings are phenomenological. In this paper, we utilize envelope function approximation and symmetry arguments to derive the spin splitting for monolayer black phosphorus, which demonstrates a highly anisotropic in-plane effective mass and, therefore, makes the phenomenological description inapplicable.

Black phosphorus, the most stable allotrope of phosphorus, is a layered material similar to graphite where van der Waals interaction binds individual layers together. Each monolayer, dubbed phosphorene, is a two-dimensional crystal with a puckered honeycomb structure which shares the symmetry properties of its bulk form denoted by the orthorhombic space group $Cmca$ Hultgren et al. (1935); Brown and Rundqvist (1965). A century after black phosphorus was discovered Bridgman (1914), phosphorene and its multilayer thin films were isolated Li et al. (2014); Koenig et al. (2014); Liu et al. (2014); Buscema et al. (2014); Castellanos-Gomez et al. (2014); Xia et al. (2014) using mechanical exfoliation, which had been utilized earlier to isolate graphene Novoselov et al. (2004). Similar to graphene, phosphorene consists of light atoms producing a spin-orbit coupling of $\sim$1 meV, which is weaker than that of conventional zinc-blende crystals. Hence, both graphene and phosphorene are expected to have a long spin-relaxation time $\sim$1 ns Tombros et al. (2007); Avsar et al. (2017), which could allow spin-polarized currents to flow macroscopic distances in these materials. Unlike graphene, which is gapless, phosphorene is a semiconductor with a direct band gap of 1.73 eV Li et al. (2017), which enables a wide control over its carrier density. Phosphorene also exhibits a large anisotropy in its band structure: the ratio of the in-plane effective mass of carriers along the armchair and zigzag directions is $\sim 0.1$ Liu et al. (2014). This band structure anisotropy is expected to result in anisotropic extrinsic spin-orbit coupling, also confirmed through first-principles calculations Popović et al. (2015); Kurpas et al. (2016). It is shown Ohkawa and Uemura (1974) that in conventional semiconductors with a zinc-blende crystal the three upper valence bands, which are made of only $p$ orbitals at the band edge, induce the spin splitting in the conduction band made of $s$ orbital. However, in the case of phosphorene, it is not clear which bands couple to the conduction band to produce spin splitting. While previous works Popović et al. (2015); Kurpas et al. (2016) on extrinsic spin-orbit coupling in phosphorene are solely based on first-principles calculations, this work focuses on analytic derivation of the spin-orbit splitting and demonstration of the bands involved from a group theoretical aspect. Utilizing the symmetry analysis of phosphorene, developed in Ref. Li and Appelbaum, 2014, we derive an effective Hamiltonian for the extrinsic spin-orbit coupling in the conduction band by Lowdin partitioning and then quantify the contributions of different bands in the coefficients of spin-orbit coupling in zigzag and armchair directions.

Using the effective Hamiltonian we then study how the anisotropy impacts spin lifetime, a key measure of spin transport properties. Spin-relaxation time characterizes the decay of the polarization of a non-equilibrium spin ensemble due to random fluctuations of a magnetic field. Extrinsic spin-orbit coupling acts as an effective momentum-dependent magnetic field and causes the spins of conduction electrons undergoing momentum scattering to relax. This mechanism, introduced by Dyakonov and Perel Dyakonov and Perel (1972), has been used to develop closed-form solutions of the spin-relaxation time in isotropic semiconductors Fabian et al. (2007); Averkiev and Golub (1999). However, when the effective mass of carriers is anisotropic, momentum scattering becomes anisotropic as well and must be accounted for numerically to determine spin-relaxation time accurately. The anisotropic spin relaxation time has been calculated before Kurpas et al. (2016) considering a constant momentum scattering. Here, we assume that the momentum scattering is anisotropic and, therefore, angle dependent. Generalizing the Dyakonov-Perel mechanism in the case of phosphorene, we account for the anisotropy of momentum scattering and calculate the spin-relaxation time for spin ensembles with different initial polarization. Our calculations assume that the temperature is much lower than the Fermi energy but much greater than the spin splitting i.e. $E_{\text{F}}\gg T\gg\Delta E_{\text{SO}}$. Therefore, only the electrons at the Fermi energy are taken into account and the spin-orbit coupling is treated as a perturbation.

II Extrinsic Spin-Orbit Coupling

Considering the two-dimensional crystal of phosphorene with a periodic lattice potential $V_{0}({\bf\it r})$ lying on the $xy$ plane with the armchair edge along the $x$-direction and the zigzag edge along the $y$-direction, the Hamiltonian in the presence of Pauli spin-orbit coupling, $H_{\text{SO}}$, and a perpendicular electric field $V(z)=-e\mathcal{E}z$ is

[TABLE]

Here, ${\bf\it p}$ is the momentum operator, ${\bf\it\sigma}$ is the vector of Pauli matrices, and $m_{0}$ is the mass of free electron . We note that the contribution of $V(z)$ in $H_{\text{SO}}$ is neglected. Using envelope function approximation Winkler (2003) and Bloch’s theorem in the $xy$ subspace, we can describe the solutions to the Schrodinger’s equation as

[TABLE]

where ${\bf\it k}_{\parallel}=(k_{x},k_{y},0)$ and ${\bf\it r}_{\parallel}=(x,y,0)$ are the in-plane wavevector and position respectively, $f_{nl\sigma}(z)$ are the envelope functions and $u_{l\sigma{\bf\it 0}}({\bf\it r})\ket{\sigma}$ are the lattice-periodic Bloch functions at the band edge (i.e. $\Gamma$ point, ${\bf\it k}_{\parallel}=0$) which provide a complete and orthonormal basis. Plugging Eq. 1 and Eq. 2 into the Schrodinger’s equation, multiplying with $\bra{\sigma^{\prime}}u_{l^{\prime}\sigma^{\prime}{\bf\it 0}}^{*}({\bf\it r})$, and integrating over the primitive unit cell of the lattice, we arrive at the eigenvalue equation for the envelope functions, $Hf(z)=Ef(z)$,

[TABLE]

where $E_{l}({\bf\it 0})=\expectationvalue{H_{0}}{l}$ are the energies at the band edge and ${\bf\it P}_{l^{\prime}\sigma^{\prime}l\sigma}=\matrixelement{l^{\prime}\sigma^{\prime}}{{\bf\it p}+\frac{\hbar}{4m_{0}c^{2}}{\bf\it\sigma}\times\nabla{V_{0}}}{l\sigma}\approx\matrixelement{l^{\prime}}{{\bf\it p}}{l}\delta_{\sigma\sigma^{\prime}}$ are approximated by the matrix elements of momentum operator which couple different bands at the edge. This approximation is valid for light atoms such as phosphorus where the contribution of spin-orbit coupling is orders of magnitude smaller than that of momentum operator. The matrix elements of $H_{\text{SO}}$ are denoted by $\Delta_{l^{\prime}\sigma^{\prime}l\sigma}=\matrixelement{l^{\prime}\sigma^{\prime}}{H_{\text{SO}}}{l\sigma}$. We rewrite the Pauli spin-orbit coupling $H_{\text{SO}}=(\xi\hbar/2){\bf\it\mathcal{L}}\dotproduct{\bf\it\sigma}$, where ${\bf\it\mathcal{L}}$ is the angular momentum operator, and the parameter $\xi$ represents the strength of atomic spin-orbit coupling and includes the average radial contribution of $\nabla V_{0}({\bf\it r})$. Therefore, we rewrite $\Delta_{l^{\prime}\sigma^{\prime}l\sigma}=(\xi\hbar/2)\matrixelement{l^{\prime}\sigma^{\prime}}{{\bf\it\mathcal{L}}\dotproduct{\bf\it\sigma}}{l\sigma}$.

The infinite-dimensional Hilbert space in Eq. 3 can be reduced to a finite-dimensional one by considering only the bands that are in the vicinity of the Fermi energy. These bands are usually made up of orbitals of the valence electrons namely $s$ and $p$ atomic orbitals in phosphorus. Using atomic orbitals as a basis, Li et al. Li and Appelbaum (2014) studied the symmetry properties of phosphorene and the orbital composition of the bands at the $\Gamma$ point. Following their study, we derive an effective Hamiltonian for the extrinsic spin-orbit coupling in the conduction band utilizing the theory of invariants and Lowdin partitioning. There are eight irreducible representations (IRs) in the space group of phosphorene, $Cmca$, denoted by $\Gamma_{i}^{\pm}$ at the $\Gamma$ point of the Brillouin zone. Each band at the $\Gamma$ point can be labeled by one of the IRs of the space group. Table 1 in the appendix lists the characters of these IRs. The invariants corresponding to each IR is also listed in this table. According to the theory of invariants, the energy of each band should be an invariant of $\Gamma_{1}^{+}$. Therefore, various products of $k_{x}$, $k_{y}$, $\sigma_{x}$, and $\sigma_{y}$ could appear in the diagonalized effective Hamiltonian as long as the direct product of their corresponding IRs results in $\Gamma_{1}^{+}$. For instance, from Table 1 we can see that terms such as $k_{z}k_{x}\sigma_{y}$ and $k_{z}k_{y}\sigma_{x}$ will be present in the effective Hamiltonian but terms such as $k_{z}k_{x}\sigma_{x}$ and $k_{z}k_{y}\sigma_{y}$ will not. Since the external electric field $V(z)$ breaks the inversion symmetry in the $z$-direction, it does not commute with $k_{z}$ and, therefore, after replacing $k_{z}$ with $-id/dz$, their commutator $[k_{z},V(z)]=-idV(z)/dz$ leads to the extrinsic spin-orbit coupling.

To derive the effective Hamiltonian for the conduction band $H_{cc}$, we first represent the Hamiltonian in Eq. 3 as $H=H_{0}+H^{\prime}$ where $H_{0}$ contains band-edge and free-electron energies and $H^{\prime}$ contains the electric field and the off-diagonal parts of $H$ namely k.p and spin-orbit terms. By applying Lowdin partitioning up to the fourth order perturbation, the coupling of conduction band with other bands is incorporated in the effective Hamiltonian. The $z$-dependence is averaged out afterwards to obtain the final two-dimensional effective Hamiltonian. The details of Lowdin partitioning and the derivation of the effective Hamiltonian is provided in Appendix A. The resulting effective Hamiltonian is

[TABLE]

where $m_{x}$ and $m_{y}$ are the in-plane effective masses written as

[TABLE]

where the sum is over $l=\Gamma_{2}^{+}$ and $l=\Gamma_{3}^{+}$ for $x$ and $y$ directions, respectively. The coefficients of extrinsic spin-orbit coupling term, $\lambda_{x}$ and $\lambda_{y}$, are of the following form

[TABLE]

where $P$ and $L$ are matrix elements of momentum and angular momentum operators that couple $c$, $l$, and $l^{\prime}$ bands according to the invariants and their corresponding IRs. For instance, for $l=\Gamma_{1}^{+}$ and $l^{\prime}=\Gamma_{2}^{+}$ ($c=\Gamma_{4}^{-}$ denotes the conduction band), the matrix elements appearing in $\lambda_{x}$ are $P_{z,cl}=\matrixelement{c}{p_{z}}{l}$, $L_{y,ll^{\prime}}=\matrixelement{l}{\mathcal{L}_{y}}{l^{\prime}}$, and $P_{x,l^{\prime}c}=\matrixelement{l^{\prime}}{p_{x}}{c}$. Equation 18 contains all possible terms that produce extrinsic spin-orbit coupling which is linear in spin-orbit strength $\xi$, electric field $dV(z)/dz$, and in-plane crystal momentum $k_{x(y)}$.

First principle calculations within density functional theory (DFT) using projected augmented plane wave method implemented in Quantum ESPRESSO package Giannozzi et al. (2017) were performed to verify the symmetry of the bands and also calculate the k.p parameters namely different $P$ and $L$ matrix elements appearing in Eq. 18. Details of the DFT setup are provided in Table 2. Parameters of the crystal structure of phosphorene are listed in Table 3. Figure 1 illustrates the band structure of phosphorene for the lowest 24 bands for a path along the high symmetry points. The bands are labeled with their corresponding IRs which are consistent with previous calculations Takao et al. (1981); Li and Appelbaum (2014).

Results show that $m_{x}=0.16m_{0}$ and $m_{y}=1.13m_{0}$ which is similar to what has been reported before Popović et al. (2015); Qiao et al. (2014). Figure 2a plots $\lambda_{x}$ and $\lambda_{y}$ versus the external electric field. Quadratic terms with respect to the electric field appear only at high magnitudes, i.e. $\mathcal{E}>0.4$ V/Å(not shown in the figure). The expectation value of spin, over the Fermi contour, is illustrated in Fig. 2b.

As seen from the figure, the spin is not tangential to the Fermi contour in contrast to the isotropic case. To verify Eq. 6 produces the same anisotropy as that of Fig. 2a, the various $P$ and $L$ matrix elements for different bands were calculated. The contributions of different bands in the spin-orbit coefficients are listed in Table 4. The band number corresponding to each IR is mentioned as well to uniquely specify the first 24 diagonalized bands ordered from the lowest energy to the highest. For instance, the conduction band is labeled as $\Gamma_{4}^{-}$ and is the 11${}^{\text{th}}$ band from the bottom in Fig. 1. As seen from the table, only few of the bands contribute significantly to the spin-orbit coefficients. Also contributions from different bands possess different signs which cancel each other out. The anisotropy of spin-orbit coupling is determined by the $\lambda_{x}/\lambda_{y}$ ratio which requires the contributions from different bands to be summed up. As seen from Table 4 this ratio is $\lambda_{x}/\lambda_{y}\approx-1.4$ which is close to the full first principle result of $\lambda_{x}/\lambda_{y}\approx-1.7$ plotted in Fig. 2. The small discrepancy can be resolved as one adds other contributions from weakly coupled bands to the sum in Eq. 6.

III Spin Relaxation

Utilizing the effective Hamiltonian in Eq. 4, we study the spin relaxation of conduction electrons. Decomposing the Hamiltonian $H_{cc}=H+H^{\prime}$, we rewrite the extrinsic spin-orbit coupling term as $H^{\prime}={\bf\it\Omega}_{{\bf\it k}}\dotproduct{\bf\it\sigma}$ where ${\bf\it\Omega}_{{\bf\it k}}=\lambda_{y}k_{y}\hat{{\bf\it x}}+\lambda_{x}k_{x}\hat{{\bf\it y}}$ is an effective magnetic field which is ${\bf\it k}$-dependent. This effective magnetic field causes the electrons with different momenta to precess around different axes. Therefore, scattering between different momenta randomizes the precession of a polarized spin ensemble and consequently leads to spin relaxation. This is the aforementioned Dyakonov-Perel mechanism.

To calculate the spin-relaxation time, we follow a similar procedure as in Refs. Averkiev et al., 2002; Fabian et al., 2007, but we specifically analyze an anisotropic Fermi contour with an anisotropic extrinsic spin-orbit coupling. Considering a polarized spin ensemble which is spatially-homogeneous and is described by a ${\bf\it k}$-dependent density matrix $\rho_{{\bf\it k}}$, the time evolution is given as the following kinetic equation Averkiev and Golub (1999)

[TABLE]

where we used $[\rho_{{\bf\it k}},H]=0$. Here $W_{{\bf\it kk^{\prime}}}$ is the probability density of transition between ${\bf\it k}$ and ${\bf\it k^{\prime}}$ states. The first term on the right-hand side represents spin precession about ${\bf\it\Omega}_{{\bf\it k}}$, and the second term represents momentum scattering between incoming wavevector ${\bf\it k}$ and outgoing wavevector ${\bf\it k^{\prime}}$. We assume that the density matrix can be decomposed as $\rho_{{\bf\it k}}=\overline{\rho}+\rho^{\prime}_{{\bf\it k}}$, where $\overline{\rho}$ is the average of density matrix over the Fermi contour, i.e. $\overline{\rho}=\ell^{-1}\int d\ell\rho_{{\bf\it k}}$, where $\ell$ is the perimeter of the Fermi contour and $d\ell=d\theta\absolutevalue{\partialderivative*{{\bf\it k}}{\theta}}$ is the differential arc length. We assume that $\rho^{\prime}_{{\bf\it k}}$ is a small perturbation with zero average, i.e. $\overline{\rho}_{{\bf\it k}}^{\prime}=0$. Taking the average of Eq. 7 over the Fermi contour, we obtain

[TABLE]

where we used the fact that $\overline{{\bf\it\Omega}_{{\bf\it k}}}$ is zero. The reason is that for each point ${\bf\it k}$ on the Fermi contour, $-{\bf\it k}$ is also on the Fermi contour. Since ${\bf\it\Omega}_{{\bf\it k}}$ is linear in ${\bf\it k}$ and therefore an odd function of ${\bf\it k}$, i.e. ${\bf\it\Omega}_{-{\bf\it k}}=-{\bf\it\Omega}_{{\bf\it k}}$, it averages to zero over the Fermi contour. Applying the decomposition to Eq. 7 and dropping the terms containing product of ${\bf\it\Omega}_{{\bf\it k}}$ and $\rho^{\prime}_{{\bf\it k}}$, we can find the quasistatic value of $\rho^{\prime}_{{\bf\it k}}$, by setting $\partialderivative*{\rho^{\prime}_{{\bf\it k}}}{t}$ to zero, assuming that momentum relaxation is much faster than spin relaxation. Therefore,

[TABLE]

Equations 8 and 9 are coupled and must be solved self-consistently. To do so, first we assume that the average spin polarization is in $\hat{{\bf\it s}}$ direction. Therefore, we can write $\overline{\rho}=(1+\hat{{\bf\it s}}\dotproduct{\bf\it\sigma})/2$. It can be shown that $[\overline{\rho},{\bf\it\Omega}_{{\bf\it k}}\dotproduct{\bf\it\sigma}]=i(\hat{{\bf\it s}}\times{\bf\it\Omega}_{{\bf\it k}})\dotproduct{{\bf\it\sigma}}$. Using Eq. 9, we can solve for $\rho^{\prime}_{{\bf\it k}}$ iteratively using the following equation:

[TABLE]

Plugging $\rho^{\prime}_{{\bf\it k}}$ into Eq. 8, we calculate the rate of decay $\partialderivative*{\overline{\rho}}{t}$ or correspondingly $\derivative*{\hat{{\bf\it s}}}{t}=-\hat{{\bf\it s}}/\tau_{s}$ which results in the spin-relaxation time $\tau_{s}$.

The collision sum in the continuum limit becomes an integral, i.e. $\sum_{{\bf\it k^{\prime}}\not={\bf\it k}}W_{{\bf\it kk^{\prime}}}\rightarrow A\int d^{2}{\bf\it k^{\prime}}(2\pi)^{-2}W_{{\bf\it kk^{\prime}}}$, where $A$ is the area of the system. Using Fermi’s golden rule, the probability density of transition is given as $W_{{\bf\it kk^{\prime}}}=\frac{2\pi}{\hbar}N\absolutevalue{U_{{\bf\it k}{\bf\it k^{\prime}}}}^{2}\delta(E({{\bf\it k}})-E({{\bf\it k^{\prime}}}))$, where $N$ is the number of scatterers and $U_{{\bf\it k}{\bf\it k^{\prime}}}$ is the matrix element of the scattering potential. For long-range scattering potential varying slowly compared to the periodic lattice potential, $U_{{\bf\it k}{\bf\it k^{\prime}}}=U({\bf\it k}-{\bf\it k^{\prime}})/A=U({\bf\it q})/A$, where $U({\bf\it q})$ is the Fourier transform of $U({\bf\it r})$. In two-dimensional electron systems, the effective Coulomb potential in the Fourier domain is Ando et al. (1982)

[TABLE]

where $\kappa$ is the average relative permittivity, $d$ is the depth of the scattering center in the substrate, and $q_{s}\approx 2\sqrt{m_{x}m_{y}}e^{2}/\kappa\hbar^{2}$ is the Thomas-Fermi screening constant. The delta function in $W_{{\bf\it kk^{\prime}}}$ reduces the $k$-space integral to an integral over the Fermi contour. Therefore,

[TABLE]

where $n=N/A$ is the density of scatterers.

The ${\bf\it k}$-dependent momentum scattering time, $\tau_{{\bf\it k}}=1/\sum_{{\bf\it k^{\prime}}}W_{{\bf\it kk^{\prime}}}$, is depicted in Fig. 3a as a function of the polar angle for a typical value of charged impurity density Yuan et al. (2015); Liu et al. (2016), i.e. $n=10^{12}$ cm${}^{\text{-2}}$.

We assume that the monolayer is deposited on an SiO2 substrate Wang et al. (2015); Li et al. (2014) with relative permittivity of $\epsilon_{\text{r}}=3.9$. As seen from the figure, the momentum scattering time shows a significant anisotropy which consequently affects the spin relaxation.

Figure 3b depicts the energy dependence of both average momentum scattering time, $\overline{\tau_{{\bf\it k}}}$, and spin-relaxation time for two ensembles initially polarized in the x-direction, $\tau_{s,xx}$, and the y-directions, $\tau_{s,yy}$. We assume that the external electric field is fixed at $\mathcal{E}=0.1$ V/Å  and the spin-orbit coefficients are $\lambda_{x}=-1.79$ meV$\cdot$Å and $\lambda_{y}=1.03$ meV$\cdot$Å  (from Fig. 2). Spin lifetime of the x-polarized ensemble is about 2-3 times longer than that of the y-polarized ensemble. This anisotropy increases with an increase in the Fermi energy. The spin-relaxation time decreases with Fermi energy, while the average momentum scattering time increases with Fermi energy. This opposite energy dependence is a signature of the Dyakonov-Perel mechanism. We note that our assumptions in deriving Eqs. 8 and 9 are valid as long as the average momentum scattering occurs on a faster timescale compared to spin relaxation, i.e. for $E<0.4$ eV. We also note that $1/\tau_{s,xy}=1/\tau_{s,yx}=0$. The spin relaxation for an ensemble polarized along the $\hat{{\bf\it z}}$ axis is always faster than in-plane directions (not shown in the figure). Replacing $\hat{{\bf\it s}}$ with $\hat{{\bf\it z}}$ in Eq. 10, we can see that $\rho^{\prime}_{{\bf\it k}}$ obtains both $\sigma_{x}$ and $\sigma_{y}$ components. Therefore, the corresponding spin relaxation rate is the sum of relaxation rates along the in-plane directions, i.e. $1/\tau_{s,zz}=1/\tau_{s,xx}+1/\tau_{s,yy}$. We note that as the energy increases, the ratio of in-plane spin-relaxation times, $\tau_{s,yy}/\tau_{s,xx}$, increases as well. The reason is that at higher energies screening becomes less effective as we can see from the screened Coulomb potential in Eq. 11 which leads to a greater anisotropy in $\tau_{{\bf\it k}}$. Finally, for an isotropic two-dimensional system, i.e. $m_{x}=m_{y}$ and $\lambda_{x}=\lambda_{y}$, we obtain $\tau_{s,xx}=\tau_{s,yy}=2\tau_{s,zz}$ which has been reported previously in the literature Fabian et al. (2007).

IV Summary

Using envelope function approximation along with symmetry arguments and Lowdin partitioning an effective Hamiltonian was derived to describe the extrinsic spin-orbit coupling for the conduction electrons in phosphorene. Based on the theory of invariants, we determined the bands that are involved in generating extrinsic spin splitting. In contrast to the isotropic Bychkov-Rashba and Vasko spin splittings, phosphorene shows an anisotropic spin splitting which is characterized by two coefficients. The coefficients are determined by the coupling between the conduction band and other bands. First-principles calculations were performed to verify the symmetry of the bands and also quantify the contribution of different bands in spin splitting of conduction electrons by calculating the parameters of the effective Hamiltonian. Given the effective Hamiltonian in the conduction band, we calculated the spin-relaxation time for a homogeneous polarized spin ensemble within a generalized Dyakonov-Perel mechanism. Our results show that spin-relaxation time is highly anisotropic in the plane of phosphorene. A spin ensemble polarized in the armchair direction (x-direction) relaxes about 2-3 times longer than a spin ensemble polarized in the zigzag direction (y-direction). The calculated spin lifetimes are comparable in magnitude to the recent experiment on phosphorene Avsar et al. (2017) which is shown to be dominated by Elliott-Yafet mechanism. However, in order for this anisotropy to be detected experimentally, the Dyakonov-Perel mechanism needs to be dominant. Therefore, a significant electric field $\mathcal{E}\sim 0.1$ V/Å  and a highly disordered sample with charged impurity density $n=10^{12}$ cm*-2* would be required.

Acknowledgements.

The authors acknowledge the funding support from the MRSEC Program of the National Science Foundation under Award Number DMR-1420073.

Appendix A Lowdin Partitioning

Lowdin partitioning Löwdin (1951), also known as quasi-degenerate perturbation theory, is a method for block diagonalizing a Hamiltonian with two sets of bands that weakly interact with each other. Here we are interested in the conduction band only and therefore calculate the effect of all other bands perturbatively as a separate set interacting with the conduction band. We assume that the Hamiltonian of the conduction band is denoted by $H_{cc}$ and the Hamiltonian in Eq. 3 is decomposed as $H=H_{0}+H^{\prime}$ where

[TABLE]

Based on Lowdin partitioning, the leading four terms that comprise $H_{cc}$ which are given as followsWinkler (2003).

[TABLE]

The zeroth order term is similar to the energy dispersion of a free electron confined in two-dimensions i.e. $H_{cc}^{(0)}=H_{cc}^{0}=E_{c}+\hbar^{2}(k_{\parallel}^{2}-(d/dz)^{2})/2m_{0}$. The first order term is similarly given by $H^{\prime}$ as $H_{cc}^{(1)}=V(z)$ which vanishes upon taking the expectation value over the $z$ direction. The expectation value is $\langle H_{cc}^{(1)}\rangle=\int dzf^{*}_{c}(z)V(z)f_{c}(z)=0$ which is a direct consequence of inversion symmetry of the unperturbed envelope function $f_{c}(z)$. More generally, higher order terms that contain odd powers of $V(z)$ or $k_{z}=-id/dz$ vanish upon averaging over $z$. This makes the calculation of the higher order terms more convenient as many terms become zero. The second order term contains the effective mass terms, subband energy in the $z$ direction, and other terms quadratic in spin-orbit strength which we discard.

[TABLE]

The first sum in the third order term vanishes upon averaging over $z$ for it contains $H^{\prime}_{cc}=V(z)$ in all of its terms. The same is true for the second sum although it is more subtle as the terms with $l=\Gamma_{1}^{+}$ are operators affecting $H^{\prime}_{cc}=V(z)$. In other words, since $c\otimes l=\Gamma_{4}^{-}\otimes\Gamma_{1}^{+}=\Gamma_{4}^{-}$, $H_{cl}$ contains $k_{z}$ invariant, i.e. $H_{cl}=(\hbar/m_{0})(-id/dz)P_{z,cl}$. Therefore,

[TABLE]

All the terms on the right hand side vanish upon averaging out $z$. Similarly, the third sum contains only odd terms in $V(z)$ and $d/dz$ which vanish as well. It also produces some high order terms in momentum and spin-orbit strength, i.e. $\langle H_{cc}^{(3)}\rangle=\mathcal{O}(\xi k_{x}k_{y})\sigma_{z}+\mathcal{O}(\xi^{3})\sigma_{0}$, which are negligible. The extrinsic spin orbit coupling, which couples electric field, in-plane momentum, and spin-orbit together, emerges from fourth order terms. The first sum in the fourth order term produces terms that are quadratic in electric field $(dV(z)/dz)^{2}$ which are negligible and appear only at high electric fields. The third sum produces other negligible high order terms i.e. $\mathcal{O}(k_{x}^{4})+\mathcal{O}(k_{y}^{4})+\mathcal{O}(k_{x}^{2}d^{2}/dz^{2})+\mathcal{O}(k_{y}^{2}d^{2}/dz^{2})+\mathcal{O}(d^{4}/dz^{4})+\mathcal{O}(\xi^{2}k_{x}^{2})+\mathcal{O}(\xi^{2}k_{y}^{2})+\mathcal{O}(\xi^{2}d^{2}/dz^{2})+\mathcal{O}(\xi^{4})$. In the second sum, only terms containing both $V(z)$ and $d/dz$ survive. Other terms are either quadratic in $V(z)$ or zero after averaging out $z$. First order terms in $V(z)$ can be found in the fourth sum where either $l=l^{\prime}$ or $l^{\prime}=l^{\prime\prime}$. Other terms stemming from the fourth sum are negligible for they contain higher orders of electric field, spin-orbit, or momentum. Therefore, we can approximate the fourth sum as

[TABLE]

Since $H_{cc}=H_{ll}=V(z)$ for all $l$, we can combine the second and the fourth sum together. It might seem that they would cancel each other completely. However, since $V(z)$ and $d/dz$ do not commute with each other, their commutator, i.e. $[d/dz,V(z)]=dV(z)/dz$, survives the summation. We can see that this is the direct consequence of the broken inversion symmetry by $V(z)$ which leads to the extrinsic spin-orbit term. Finally, based on the multiplication of group elements we obtain the non-zero terms that couple electric field, momentum in the $z$ direction, an in-plane momentum, and the angular momentum perpendicular to the rest. Therefore, we obtain the effective Hamiltonian for extrinsic spin-orbit coupling averaged over $z$ direction.

[TABLE]

As expected, this term is linear in spin-orbit strength $\xi$, electric field $dV(z)/dz$, and in-plane momentum $k_{x}$, $k_{y}$. Using first principle calculations we quantify the terms in Eq. 18 by calculating the matrix elements of momentum and angular momentum operator. The importance of each term is determined by the strength of coupling between the bands involved and also the energy difference between them. Table 4 lists all significant nonzero terms and their corresponding coefficients. The sum over all contributions (last row) determines the anisotropy $\sum\alpha_{y,ll^{\prime}}/\sum\alpha_{x,ll^{\prime}}\approx-1.4~{}$.

t