Extrinsic Spin-Orbit Coupling and Spin Relaxation in Phosphorene
S. M. Farzaneh Shaloo Rakheja


APS/123-QED

]Extrinsic Spin-Orbit Coupling and Spin Relaxation in Phosphorene


farzaneh@nyu.edu
[]

Department of Electrical and Computer Engineering, New York University, Brooklyn, NY 11201





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. 

[
    Shaloo Rakheja
    December 30, 2023
=====================





§ 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  <cit.>, was originally derived by Ohkawa and Uemura <cit.> for an inversion layer of zinc-blende crystals. Later, Vasko <cit.>, and Bychkov and Rashba <cit.>, 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 <cit.>, 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 <cit.>.
A century after black phosphorus was discovered <cit.>, phosphorene and its multilayer thin films were isolated <cit.> using mechanical exfoliation, which had been utilized earlier to isolate graphene <cit.>.
Similar to graphene, phosphorene consists of light atoms producing a spin-orbit coupling of ∼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 ∼1 ns <cit.>, 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 <cit.>, 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 ∼ 0.1 <cit.>.
This band structure anisotropy is expected to result in anisotropic extrinsic spin-orbit coupling, also confirmed through first-principles calculations <cit.>.
It is shown <cit.> 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 <cit.> 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. li2014electrons, 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 <cit.>, has been used to develop closed-form solutions of the spin-relaxation time in isotropic semiconductors <cit.>.
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 <cit.> 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_F≫ T ≫Δ E_SO. Therefore, only the electrons at the Fermi energy are taken into account and the spin-orbit coupling is treated as a perturbation. 




§ EXTRINSIC SPIN-ORBIT COUPLING


Considering the two-dimensional crystal of phosphorene with a periodic lattice potential V_0(*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_SO, and a perpendicular electric field V(z)=-eℰz is 

    ℋ = p^2/2m_0 + V_0(*r)_H_0 + ħ^2/4m_0^2c^2*p*σ×V_0(*r)_H_SO + V(z)·
 
Here, *p is the momentum operator, *σ 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_SO is neglected. 
Using envelope function approximation <cit.> and Bloch's theorem in the xy subspace, we can describe the solutions to the Schrodinger's equation as

    Ψ_n(*r) = e^i*k_∥*r_∥∑_l,σf_nlσ(z)u_lσ*0(*r)|σ⟩,
 
where *k_∥=(k_x, k_y, 0) and  *r_∥=(x, y, 0) are the in-plane wavevector and position respectively, f_nlσ(z) are the envelope functions and u_lσ*0(*r)|σ⟩ are the lattice-periodic Bloch functions at the band edge (i.e. Γ point, *k_∥=0) which provide a complete and orthonormal basis. 
Plugging Eq. <ref> and Eq. <ref> into the Schrodinger's equation, multiplying with ⟨σ'|u_l'σ'*0^*(*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),

    ∑_l,σ[(E_l(*0) + ħ^2(k_∥^2 - d^2/dz^2)/2m_0 + V(z))δ_ll'δ_σσ' + Δ_l'σ'lσ
        + ħ/m_0(*k_∥ -  id/dz*z)*P_l'σ'lσ]f_nlσ(z) = E_nf_nl'σ'(z),

where E_l(*0)=H_0l are the energies at the band edge and *P_l'σ'lσ = l'σ'*p + ħ/4m_0c^2*σ×∇V_0lσ≈l'*plδ_σσ' 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_SO are denoted by Δ_l'σ'lσ = l'σ'H_SOlσ.
We rewrite the Pauli spin-orbit coupling H_SO=(ξħ/2)*ℒ*σ, where *ℒ is the angular momentum operator, and the parameter ξ represents the strength of atomic spin-orbit coupling and includes the average radial contribution of ∇ V_0(*r).
Therefore, we rewrite Δ_l'σ'lσ = (ξħ/2)l'σ'*ℒ*σlσ.


The infinite-dimensional Hilbert space in Eq. <ref> 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. <cit.> studied the symmetry properties of phosphorene and the orbital composition of the bands at the Γ 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 Γ_i^± at the Γ point of the Brillouin zone. 
Each band at the Γ point can be labeled by one of the IRs of the space group.  
Table <ref> 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 Γ_1^+. Therefore, various products of k_x, k_y, σ_x, and σ_y could appear in the diagonalized effective Hamiltonian as long as the direct product of their corresponding IRs results in Γ_1^+. 
For instance, from Table <ref> we can see that terms such as k_zk_xσ_y and  k_zk_yσ_x will be present in the effective Hamiltonian but terms such as k_zk_xσ_x and  k_zk_yσ_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. <ref> as H = H_0 + H' where H_0 contains band-edge and free-electron energies and H' 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 <ref>. 
The resulting effective Hamiltonian is 

    H_cc = ħ^2k_x^2/2m_x + ħ^2k_y^2/2m_y + λ_xk_xσ_y + λ_yk_yσ_x,

where m_x and m_y are the in-plane effective masses written as

    1/m_x(y) = 1/m_0(1 + 2/m_0∑_l|P_x(y),cl|^2/E_c - E_l),

where the sum is over l=Γ_2^+ and l=Γ_3^+ for x and y directions, respectively. 
The coefficients of extrinsic spin-orbit coupling term, λ_x and λ_y, are of the following form 

    λ_x(y) =     (ξħ/2)dV(z)/dzħ^2/m_0^2
    ×   ∑_l,l'P_zL_y(x)P_x(y)(E_c - E_l') ± (E_c - E_l)/(E_c - E_l)^2(E_c - E_l')^2,

where P and L are matrix elements of momentum and angular momentum operators that couple c, l, and l' bands according to the invariants and their corresponding IRs. 
For instance, for l=Γ_1^+ and l'=Γ_2^+ (c=Γ_4^- denotes the conduction band), the matrix elements appearing in λ_x are P_z,cl=cp_zl, L_y,ll'=lℒ_yl', and  P_x,l'c=l'p_xc. 
Equation <ref> contains all possible terms that produce extrinsic spin-orbit coupling which is linear in spin-orbit strength ξ, 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 <cit.> 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. <ref>. 
Details of the DFT setup are provided in Table <ref>. 
Parameters of the crystal structure of phosphorene are listed in Table <ref>. 
Figure <ref> 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 <cit.>.
  
Results show that m_x=0.16m_0 and m_y=1.13m_0 which is similar to what has been reported before <cit.>.
Figure <ref>a plots λ_x and λ_y versus the external electric field. 
Quadratic terms with respect to the electric field appear only at high magnitudes, i.e. ℰ>0.4 V/Å(not shown in the figure).
The expectation value of spin, over the Fermi contour, is illustrated in Fig. <ref>b.
  
As seen from the figure, the spin is not tangential to the Fermi contour in contrast to the isotropic case.
To verify Eq. <ref> produces the same anisotropy as that of Fig. <ref>a, 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 <ref>. 
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 Γ_4^- and is the 11^th band from the bottom in Fig. <ref>. 
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 λ_x/λ_y ratio which requires the contributions from different bands to be summed up. As seen from Table <ref> this ratio is λ_x/λ_y≈ -1.4 which is close to the full first principle  result of λ_x/λ_y≈ -1.7 plotted in Fig. <ref>. 
The small discrepancy can be resolved as one adds other contributions from weakly coupled bands to the sum in Eq. <ref>.



§ SPIN RELAXATION

Utilizing the effective Hamiltonian in Eq. <ref>, we study the spin relaxation of conduction electrons. 
Decomposing the Hamiltonian H_cc = H + H', we rewrite the extrinsic spin-orbit coupling term as H'=*Ω_*k*σ where *Ω_*k=λ_yk_y*x + λ_xk_x*y is an effective magnetic field which is *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. averkiev2002spin, fabian2007semiconductor, 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 *k-dependent density matrix ρ_*k, the time evolution is given as the following kinetic equation <cit.> 

    ρ_*kt = -1/iħ[ρ_*k, *Ω_*k*σ] -
        ∑_*k'≠* kW_*kk'(ρ_*k - ρ_* k') ,

where we used [ρ_*k, H] = 0. Here W_*kk' is the probability density of transition between *k and *k' states. The first term on the right-hand side represents spin precession about *Ω_*k, and the second term represents momentum scattering between incoming wavevector *k and outgoing wavevector *k'. We assume that the density matrix can be decomposed as ρ_*k=ρ + ρ'_* k, where ρ is the average of density matrix over the Fermi contour, i.e. ρ = ℓ^-1∫ dℓρ_* k, where ℓ is the perimeter of the Fermi contour and dℓ=dθ**kθ is the differential arc length.
We assume that ρ'_*k is a small perturbation with zero average, i.e. ρ_*k'=0. 
Taking the average of Eq. <ref> over the Fermi contour, we obtain

    ρt = 
        1/iħ[ρ'_*k, *Ω_*k*σ] ,

where we used the fact that *Ω_*k is zero. The reason is that for each point *k on the Fermi contour, -*k is also on the Fermi contour. Since *Ω_*k is linear in *k and therefore an odd function of *k, i.e. *Ω_-*k=-*Ω_*k, it averages to zero over the Fermi contour.

Applying the decomposition to Eq. <ref> and dropping the terms containing product of *Ω_*k and ρ'_*k, we can find the quasistatic value of ρ'_*k, by setting *ρ'_*kt to zero, assuming that momentum relaxation is much faster than spin relaxation.
Therefore, 

    1/iħ[ρ, *Ω_*k*σ] =
        ∑_*k'≠*kW_*kk'(ρ'_*k - ρ'_* k') ·

Equations <ref> and <ref> are coupled and must be solved self-consistently. 
To do so, first we assume that the average spin polarization is in *s direction. 
Therefore, we can write ρ = (1 + *s*σ)/2.
It can be shown that [ρ, *Ω_*k*σ] = i(* s×*Ω_*k)* σ.
Using Eq. <ref>, we can solve for ρ'_*k iteratively using the following equation:

    ρ'_*k = 1/ħ(*s×*Ω_*k)* σ + ∑_*k'≠*kW_*kk'ρ'_*k'/∑_*k'≠*kW_*kk' ·

Plugging ρ'_*k into Eq. <ref>, we calculate the rate of decay * ρt or correspondingly **st=-*s/τ_s which results in the spin-relaxation time τ_s.

The collision sum in the continuum limit becomes an integral, i.e. ∑_*k'≠*kW_*kk'→ A∫ d^2*k'(2π)^-2W_*kk', where A is the area of the system. 
Using Fermi's golden rule, the probability density of transition is given as W_*kk'=2π/ħNU_*k*k'^2δ(E(* k) - E(*k')), where N is the number of scatterers and U_*k*k' is the matrix element of the scattering potential. 
For long-range scattering potential varying slowly compared to the periodic lattice potential, U_*k*k'=U(*k - *k')/A = U(*q)/A, where U(*q) is the Fourier transform of U(*r).
In two-dimensional electron systems, the effective Coulomb potential in the Fourier domain is <cit.>

    U(*q) = 2π e^2/κ(q + q_s)e^-qd,

where κ is the average relative permittivity, d is the depth of the scattering center in the substrate, and q_s≈ 2√(m_xm_y)e^2/κħ^2 is the Thomas-Fermi screening constant. 
The delta function in W_*kk' reduces the k-space integral to an integral over the Fermi contour. Therefore,

    ∑_*k'≠*kW_*kk'→n/2πħ∫ dℓ'U(* q)^2/ E(*k') ,

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

The *k-dependent momentum scattering time, τ_*k=1/∑_*k'W_*kk', is depicted in Fig. <ref>a as a function of the polar angle for a typical value of charged impurity density <cit.>, i.e. n=10^12 cm^-2. 

We assume that the monolayer is deposited on an SiO_2  substrate <cit.> with relative permittivity of ϵ_r=3.9.
As seen from the figure, the momentum scattering time shows a significant anisotropy which consequently affects the spin relaxation.

Figure <ref>b depicts the energy dependence of both average momentum scattering time, τ_*k, and spin-relaxation time for two ensembles initially polarized in the x-direction, τ_s,xx, and the y-directions, τ_s,yy.
We assume that the external electric field is fixed at ℰ=0.1 V/Å  and the spin-orbit coefficients are λ_x=-1.79 meV·Å and λ_y=1.03 meV·Å  (from Fig. <ref>). 
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. <ref> and <ref> 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/τ_s,xy=1/τ_s,yx=0.
The spin relaxation for an ensemble polarized along the *z axis is always faster than in-plane directions (not shown in the figure). 
Replacing *s with *z in Eq. <ref>, we can see that ρ'_*k obtains both σ_x and σ_y components. Therefore, the corresponding spin relaxation rate is the sum of relaxation rates along the in-plane directions, i.e. 1/τ_s,zz = 1/τ_s,xx + 1/τ_s,yy.

We note that as the energy increases, the ratio of in-plane spin-relaxation times, τ_s, yy/τ_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. <ref> which leads to a greater anisotropy in τ_*k. 
Finally, for an isotropic two-dimensional system, i.e. m_x=m_y and λ_x=λ_y, we obtain τ_s,xx=τ_s,yy=2τ_s,zz which has been reported previously in the literature <cit.>.




§ 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 <cit.> 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 ℰ∼0.1 V/Å  and a highly disordered sample with charged impurity density n=10^12 cm^-2 would be required. 





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




§ LOWDIN PARTITIONING



1.2




Lowdin partitioning <cit.>, 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. <ref> is decomposed as H=H_0 + H' where 

    H_0,ll' = (E_l + ħ^2(k_∥^2 - (d/dz)^2/2m_0)δ_ll',
    
        H_ll'' = V(z)δ_ll' + Δ_ll' + ħ/m_0(*k_∥ - id/dz)·*P_ll's·

Based on Lowdin partitioning, the leading four terms that comprise H_cc which are given as follows<cit.>. 

    
            H_cc^(0) = H_0,cc,   
    
    
            H_cc^(1) = H_cc',    
    
    
            H_cc^(2) = ∑_lH'_clH'_lc/E_c - E_l,    
    
    
        H_cc^(3) =   -1/2     ∑_lH'_ccH'_clH'_lc/(E_c - E_l)^2 -1/2∑_lH'_clH'_lcH'_cc/(E_c - E_l)^2       
+    ∑_l,l'H'_clH'_ll'H'_l'c/(E_c - E_l)(E_c - E_l')
    
    
        H_cc^(4) = 1/2    ∑_lH'_ccH'_ccH'_clH'_lc + H'_clH'_lcH'_ccH'_cc/(E_c - E_l)^3        
-1/2    ∑_l,l'(H'_clH'_ll'H'_l'cH'_cc + H'_ccH'_clH'_ll'H'_l'c)        
     ×((E_c - E_l') + (E_c - E_l)/(E_c - E_l)^2(E_c - E_l')^2)        
 -1/2    ∑_l,l'H'_clH'_lcH'_cl'H'_l'c((E_c - E_l') + (E_c - E_l)/(E_c - E_l)^2(E_c - E_l')^2)        
 +    ∑_l,l',l”H'_clH'_ll'H'_l'l”H'_l”c/(E_c - E_l)(E_c - E_l')(E_c - E_l”)
    

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 + ħ^2(k_∥^2 - (d/dz)^2)/2m_0. 
The first order term is similarly given by H' as H_cc^(1) = V(z) which vanishes upon taking the expectation value over the z direction. 
The expectation value is ⟨ H_cc^(1)⟩=∫ 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. 

    H_cc^(2) = ħ^2/m_0^2[k_x^2    ∑_l=Γ_2^+|P_x,cl|^2/E_c - E_l 
        + k_y^2∑_l=Γ_3^+|P_y,cl|^2/E_c - E_l
     
        - d^2/dz^2   ∑_l=Γ_1^+|P_z,cl|^2/E_c - E_l] + 𝒪(ξ^2)σ_0·

The first sum in the third order term vanishes upon averaging over z for it contains H'_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=Γ_1^+ are operators affecting H'_cc=V(z). In other words, since c⊗ l = Γ_4^-⊗Γ_1^+=Γ_4^-, H_cl contains k_z invariant, i.e.  H_cl=(ħ/m_0)(-id/dz)P_z,cl. Therefore, 

    H'_clH'_lcH'_cc = - ħ^2|P_z,cl|^2/m_0^2(d^2V(z)/dz^2 
        + 2dV(z)/dzd/dz
    + V(z)d^2/dz^2)·

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. ⟨ H_cc^(3)⟩=𝒪(ξ k_xk_y)σ_z + 𝒪(ξ^3)σ_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.  𝒪(k_x^4) + 𝒪(k_y^4) + 𝒪(k_x^2d^2/dz^2) + 𝒪(k_y^2d^2/dz^2) + 𝒪(d^4/dz^4) + 𝒪(ξ^2k_x^2) + 𝒪(ξ^2k_y^2)+ 𝒪(ξ^2d^2/dz^2) + 𝒪(ξ^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' or l'=l”. 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 

    ≈∑_l,l”H'_clH'_llH'_ll”H'_l”c/(E_c - E_l)^2(E_c - E_l”) + 
        ∑_l,l'H'_clH'_ll'H'_l'l'H'_l'c/(E_c - E_l)(E_c - E_l')^2·

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. 

    H_cc^(4)   ≈ (ξħ/2) dV(z)/dzħ^2/m_0^2[ 
      
        k_xσ_y    ∑_l=Γ_1^+
    l'=Γ_2^+P_z,clL_y,ll'P_x,l'c(E_c - E_l') + (E_c - E_l)/(E_c - E_l)^2(E_c - E_l')^2
    
        + k_xσ_y    ∑_l=Γ_1^+
    l'=Γ_3^-P_z,clP_x,ll'L_y,l'c(E_c - E_l') + (E_c - E_l)/(E_c - E_l)^2(E_c - E_l')^2
    
        - k_xσ_y    ∑_l=Γ_2^+
    l'=Γ_3^-P_x,clP_z,ll'L_y,l'c(E_c - E_l') - (E_c - E_l)/(E_c - E_l)^2(E_c - E_l')^2
    
        + k_yσ_x    ∑_l=Γ_1^+
    l'=Γ_3^+P_z,clL_x,ll'P_y,l'c(E_c - E_l') + (E_c - E_l)/(E_c - E_l)^2(E_c - E_l')^2
    
        + k_yσ_x    ∑_l=Γ_1^+
    l'=Γ_2^-P_z,clP_y,ll'L_x,l'c(E_c - E_l') + (E_c - E_l)/(E_c - E_l)^2(E_c - E_l')^2
    
        - k_yσ_x    ∑_l=Γ_3^+
    l'=Γ_2^-P_y,clP_z,ll'L_x,l'c(E_c - E_l') - (E_c - E_l)/(E_c - E_l)^2(E_c - E_l')^2]
        = (ξħ/2)dV(z)/dzħ^2/m_0^2∑_l,l'α_x,ll' k_xσ_y + α_y,ll' k_yσ_x·

As expected, this term is linear in spin-orbit strength ξ, electric field dV(z)/dz, and in-plane momentum k_x, k_y.  
Using first principle calculations we quantify the terms in Eq. <ref> 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 <ref> lists all significant nonzero terms and their corresponding coefficients.
The sum over all contributions (last row) determines the anisotropy  ∑α_y,ll'/∑α_x,ll'≈ -1.4.



t
1.2





    






