Spin drift-diffusion for two-subband quantum wells
SSpin drift-diffusion for two-subband quantum wells
I. R. de Assis, R. Raimondi, and G. J. Ferreira Instituto de F´ısica, Universidade Federal de Uberlˆandia, Uberlˆandia, MG 38400-902, Brazil Dipartimento di Matematica e Fisica, Universit`a degli Studi Roma Tre, Via della Vasca Navale 84, 00146 Rome, Italy (Dated: February 10, 2021)Controlling the spin dynamics and spin lifetimes is one of the main challenges in spintronics. Tothis end, the study of the spin diffusion in two dimensional electron gases (2DEGs) shows that whenthe Rashba and Dresselhaus spin-orbit couplings (SOC) are balanced, a persistent spin helix regimearises. There, a striped spin pattern shows a long lifetime, limited only by the cubic DresselhausSOC, and its dynamics can be controlled by in-plane drift fields. Here, we derive a spin diffusionequation for non-degenerate two-subbands 2DEGs. We show that the intersubband scattering rate,which is defined by the overlap of the subband densities, enters as a new nob to control the spindynamics, and can be controlled by electric fields, being maximum for symmetric quantum wells.We find that for large intersubband couplings the dynamics follows an effective diffusion matrixgiven by approximately half of the subband-averaged matrices. This extra 1/2 factor arises fromthe Matthiessen’s rule summing over the intra- and intersubband scattering rates, and leads toa reduced diffusion constant and larger spin lifetimes. We illustrate our findings with numericalsolutions of the diffusion equation with parameters extracted from realistic Schr¨odinger-Poissoncalculations.
I. INTRODUCTION
The Datta-Das spin transistor [1, 2] has establishedthe full electric control of the spin dynamics as one ofthe main goals of spintronics [3–9]. To this end, one ex-plores the spin-orbit couplings (SOCs) [10], which canbe seen as an effective momentum-dependent magneticfield. For two-dimensional electron gases (2DEGs) hostedin quantum wells (QWs) of zinc blend semiconductors,the SOCs arise from the bulk and structural inversionasymmetries, yielding the Dresselhaus [11] and Rashba[12] terms, respectively. While these allow for the spincontrol via external fields, they also introduce spin relax-ation in diffusive systems, e.g. , via the Dyakonov-Perel(DP) [13] and Elliott-Yafet (EY) [14, 15] mechanisms.Nevertheless, when the linear-in-momentum Rashba andDresselhaus SOCs are balanced, the effective SOC fieldbecomes uniaxial and it emerges a conserved SU(2) sym-metry [16, 17], yielding an helical spin mode with longlifetime, named persistent spin helix (PSH). The PSH hasbeen observed experimentally by means of spin-gratingspectroscopy [18, 19], and via time- and spatially-resolvedmagneto-optical Kerr rotation [20, 21]. The first observesa spin-lifetime enhancement when the Rashba SOC istuned to match the linear Dresselhaus SOC, while the lat-ter maps the actual formation of the PSH. The PSH dy-namics is described theoretically by a spin drift-diffusionequation [17, 22], derived from the Keldysh formalism forthe kinetic equation [23–25], and an intuitive picture canbe drawn from a random walk process [26, 27]. For sin-gle subband 2DEGs, these have been extensively studied[28–36]. A detailed derivation of the spin diffusion equa-tion can be found in Ref. [34, 36, 37]. Experimentally,the PSH has also been recently exploited in many differ-ent forms. Applying in-plane fields induces a controllableprecession frequency linear in the cubic Dresselhaus SOCand drift velocity [38, 39]. The pitch of the spin helix can be controlled by tuning the QW asymmetry and electrondensity [40]. The high symmetry of the PSH uniaxialspin regime allows for a perturbative approach to mea-sure weak (anti)-localization [41]. The spin relaxationanisotropy can be enhanced by tilting the precession an-gle with magnetic fields [42]. Under intense pump fieldswith short delay, the anisotropy might also be induced byscattering with holes [43]. For a comprehensive review ofthe PSH, see Ref. [44].For two-subbands QWs [45–48], the spin drift-diffusionhas only been studied recently [27, 49–52]. Particularly,it has been observed strong spin relaxation anisotropywith long lifetimes [52]. Theoretically, it has been pre-dicted that the Rashba and Dresselhaus SOCs from dis-tinct subbands can be set to match along perpendic-ular axes [49], indicating the formation of a persistentSkyrmion lattice (PSL). Even though, the PSH is a ro-bust effect in a single-subband system, the role of in-tersubband scattering in a two-subband 2DEG has beenshown to be relevant for the evaluation of the transportproperties [53]. It is then an open problem whether thePSH and PSL are robust in the presence of intersubbandscattering [27]. The present paper aims to shed light onthis problem, by showing that the intersubband scatter-ing reveal itself as a new knob to actually control the spinlifetime in quantum well devices.In this paper we derive and analyze the spin drift-diffusion equation for two-subbands 2DEGs accountingfor both intra- and intersubband scattering. The deriva-tion follows from the the Keldysh formalism for the ki-netic equation [24, 25, 36, 37]. We show that the intersub-band scattering introduces a time scale t c (to be properlydefined later) for the subband relaxation, and its com-parison with the momentum relaxation time τ and spinlifetime τ S allows for a classification of the spin dynamicsbetween two extreme limits that match the intuitive pic-tures taken from a random walk process [27]. First, if t c a r X i v : . [ c ond - m a t . m e s - h a ll ] F e b is comparable to τ S the intersubband dynamics is slow,defining a weak-coupling regime where electrons fromeach subband diffuse and precesse in time independently.On the other hand, for t c (cid:28) τ S the fast intersubbanddynamics leads to a subband-averaged strong-couplingregime. In both limits, and for strong enough disorder,the spin relaxation time τ S follows the Dyakonov-Perel(DP) behavior τ S ∝ /τ . For intermediate regimes, theintersubband coupling breaks this proportionality, sincewe find t c ∝ τ . Consequently, the competition betweenthe DP mechanism and the EY-like subband relaxation t c leads to a non-monotonic transition between strong- andweak-coupling regimes as τ increases. Therefore, thecontrol over the diffusion regime requires a careful choiceof the disorder strength and of the degree of asymmetryof the quantum well.Applying these results to the PSL system fromRef. [49], we find that the Skyrmion checkerboard patternmay exist only in the weak-coupling regime. An intuitivephysical picture can be understood by recalling that, asalso remarked in Ref. [49], the PSL inherits its robust-ness from the PSHs in the separate subbands. From thispoint of view, it is crucial to achieve a weak-coupling orintermediate regime where the spin dynamics is faster than the intersubband dynamics. As shown in SectionIV D, we are able to demonstrate how the checkerboardspin pattern, initially killed by the switching on of theintersubband scattering, rises again upon increasing themomentum relaxation time.The rest of the paper is organized as follows. In Sec-tion II we introduce the two-subband model with bothRashba and Dresselhaus SOC, and as well as intra- andintersubband disorder scattering. In Section III we derivethe spin drift-diffusion equation for two non-degeneratesubbands. In Section IV we present our results with anextended discussion of the emergence of the PSL in atwo-subband quantum well device. We provide our con-clusions in Section V. Finally a number of technical ap-pendices give further details about our derivation andsolution of the diffusive equation. II. 2DEG MODEL
Let us start with a clean two-dimensional electron gas(2DEG) with two-subbands given by the 2D Hamiltonian[45–49], H = (cid:18) H H H † H (cid:19) ≈ (cid:18) H H (cid:19) , (1)where H j are 2 × j = { , } , which includes the Rashba and Dressel-haus spin-orbit couplings (SOC). This 2D model is de-fined in the basis of the quantum well eigenstates foreach subband, i.e. (cid:104) r | j k σ (cid:105) = ϕ j ( z )[ e i k · r / √ A ] ξ σ , where k = ( k x , k y ) is the in-plane wave-vector, ξ σ are the spineigenvectors with σ = {↑ , ↓} , A is the normalization area,and ϕ j ( z ) is the confined eigenstates of the quantum (a) (b)(c) (d) Figure 1. (a) Symmetric GaAs double quantum well with twooccupied subbands. The dashed line marks the Fermi level at ε F = 0 for a total density n T = 8 × cm − , and the ϕ j ( z ) envelope functions are shown for the first ( j = 1, blue)and second ( j = 2, red) subbands. (b) Fermi circles of bothsubbands labeled by their color with respect to panel (a). Theblurring illustrates a small impurity broadening. The insetillustrates the parabolic dispersion with a small SOC. (c) Fora quantum well tilted by an electric field F z = 0 . µ m the ϕ ( z ) [ ϕ ( z )] shifts to the right [left] edges of the well, whichincreases the energy splitting, as shown in panel (d). well, which are illustrated in Fig. 1(a). We are inter-ested in a nondegenerate regime, i.e., the subband ener-gies ε (0)1 (cid:54) = ε (0)2 and the Fermi circles do not overlap, asshown in Fig. 1(b). Therefore we can neglect off diagonalintersubband SOC in H [49]. The 2D Hamiltonian foreach subband reads H j = (cid:104) ε j ( k ) − e F · r (cid:105) σ + H soc j ( k ) , (2) H soc j = λ j, + σ x k y + λ j, − σ y k x + 2 β ,j k x − k y k ( k y σ x − k x σ y ) , (3)where ε j ( k ) = ε (0) j + (cid:126) k / m , F = ( F x , F y ,
0) is an ex-ternal in-plane electric field, σ = × and σ x | y | z are thePauli matrices in spin space. H soc j contains the Rashba( α j ), linear ( β ,j ) and cubic ( β ,j = γ D k /
4) Dresselhausintrasubband couplings, with λ j, ± = β ,j ± α j . As pre-sented later in Section IV, the SOC and other parameterscan be controlled by the z-component of the electric field F z .For the spin diffusion, in the next section, we considerscattering by random short-range impurities. Therefore,we must add it to the clean H of Eq. (1) as a per-turbation. The 3D form of the impurity potential is V ( r ) = (cid:80) i v δ ( r − R i ), where v is the intensity and R i = ( X i , Y i , Z i ) is the position of each of the N randomimpurities, labeled by i = { , N } . The 2D projection of V ( r ) into the basis | j k σ (cid:105) reads (cid:104) j (cid:48) k (cid:48) σ (cid:48) | V ( r ) | j k σ (cid:105) = δ σ,σ (cid:48) N (cid:88) i =1 ˜ v j (cid:48) ,j ( Z i ) e i ( k − k (cid:48) ) · R i , (4)˜ v j (cid:48) ,j ( Z i ) = v A ϕ † j (cid:48) ( Z i ) ϕ j ( Z i ) . (5)Notice that the impurities allow for subband scattering,since V = (cid:104) k σ | V ( r ) | k σ (cid:105) (cid:54) = 0. In practice, we areinterested in the impurity ensemble self-average [23, 24,54], which reads (cid:104) V ba (cid:105) = δ b,a n imp v ≡ , (6) (cid:104) V b,c V c,a (cid:105) = δ k b ,k a n imp v A Λ j a ,j c ≈ δ b,a n imp v A Λ j a ,j c , (7)Λ j a ,j c = (cid:90) | ϕ j a ( z ) | | ϕ j c ( z ) | dz, (8)where we have introduced a compact notation for theindices as a → ( j a , k a , σ a ), and n imp is the impurity den-sity. The first order average, (cid:104) V ba (cid:105) ≡
0, is set as theenergy reference. For more details on the derivation ofthese quantities, see the Appendix A. The second order (cid:104) V b,c V c,a (cid:105) conserves energy and momentum, but it is notblock diagonal in subband space a priori . However, inthe non-degenerate scenario the Fermi circles from dif-ferent subbands do not overlap, see Fig. 1(b), and δ k b ,k a implies subband conservation, yielding the full δ b,a in theapproximate expression in Eq. (7). Nevertheless, it stillallows for a virtual jump to a different subband j c , whichwill lead to the intersubband coupling on the collisionintegral of the kinetic equation, as shown later in Eqs.(12) and (15). The density overlap integral Λ j,j (cid:48) definesthe intensity of this coupling. III. SPIN DIFFUSION WITH TWO SUBBANDS
To derive the spin drift-diffusion equation for two sub-band quantum wells, we follow the kinetic equation ap-proach [23, 24, 34, 36] and extend the formalism to thetwo subband model. Our starting point is the well-knownquantum Boltzmann transport equation for a Wigner dis-tribution function ρ ( k , R , T ), which is a 4 × R , T ) are the Wigner co-ordinates [24, 25]. Here we take a block diagonal approx-imation with ρ = ρ ⊕ ρ , where ρ j ( k , R , T ) are 2 × j , and obeys the kinetic equation ∂ T ρ j − i (cid:126) [ H j , ρ j ] − (cid:126) {∇ R H j , ∇ k ρ j } + 12 (cid:126) {∇ k H j , ∇ R ρ j } = I coll j , (9) I coll = i π (cid:126) (cid:90) dE (cid:16) { A, Σ K } − { Γ , G K } (cid:17) . (10)The block diagonal approximation is valid in the non-degenerate regime presented in the previous section,where the Fermi circles of each subband are distinguish-able, see Fig. 1(b), and it is consistent with the self-consistent Born approximation (SCBA), as shown nextas we present the collision integral I coll , which will couple ρ and ρ through intersubband impurity scattering.The I coll above is defined in terms of the spectralfunction A ( E, k ) = i [ G R ( E, k ) − G A ( E, k )], broaden-ing Γ( E, k ) = i [Σ R ( E, k ) − Σ A ( E, k )], and the Keld-ish component of the self-energy Σ K ( E, k , R , T ) andGreen’s function G K ( E, k , R , T ). The retarded and ad-vanced Green’s functions are block diagonal G R ( A ) = G R ( A )1 ⊕ G R ( A )2 , with G R ( A ) j ≈ [ E − ε j ( k ) ∓ iη ] − for eachsubband j . For simplicity, here we neglect small SOCcontributions, while we show these perturbative correc-tions Appendix B. Thus, the spectral function becomes A = A ⊕ A , with A j ≈ πδ ( E − ε j ( k )). The G K plays acentral role in deriving Eq. (9), as shown in Refs. [24, 25],with ρ ( k , R , T ) = − (cid:82) dE πi G K ( E, k , R , T ).The last missing pieces to characterize our system arethe self-energies Σ R , Σ A and Σ K . To derive these, wefollow the standard impurity self-average within the self-consistent Born approximation (SCBA) [54] and extendit to the two-subbands scenario using the impurity aver-aged potentials from Eq. (7), yielding for each subband j Σ νj = n imp v A (cid:88) c Λ j,j c (cid:104) G νc (cid:105) , (11)which is written generically for all components ν = R, A, K . The block diagonal form here, Σ ν = Σ ν ⊕ Σ ν isa consequence of the non-degenerate scenario introducedin Eq. (7), and it justifies the block diagonal form of G K and ρ , since it allows for the decoupling of the blockdiagonal and non-diagonal components. This approxi-mation will fail if the Fermi circles from different sub-bands overlap with the range defined by the broadeningof the spectral function. In this case, the intersubbandSOC [27, 45–47, 49] in H could also play a significantrole.Collecting the approximations above, we calculate theleading order contribution to the collision integral I (0) = I (0)1 ⊕ I (0)2 , where I (0) j ( k ) = (cid:88) j c (cid:104) ρ j c ( ε j ( k ) , R , T ) (cid:105) − ρ j ( ε j ( k ) , θ, R , T ) τ j,j c , (12) k = ( k cos θ, k sin θ ) is now in polar coordinates, and (cid:104) ρ (cid:105) = (2 π ) − (cid:82) ρdθ is the θ -average over the Fermi circle.The intra- and inter-subband momentum relaxationrates are 1 τ j,j c = φ j,j c τ , (13)where τ − = Λ , n imp mv / (cid:126) is the first subband mo-mentum relaxation rate, and the overlap integral ratio φ j,j c = Λ j,j c / Λ , defines the relative intensity betweenthe intra- and intersubband scattering times ( τ , , τ , ,and τ , ). Essentially, φ j,j c measures the extension of thewave-functions of the quantum well ϕ j ( z ), Fig. 1, whichdefines how many impurities couple the initial and finalstates. In the next section, the spin diffusion will bewritten in terms of1 τ j = (cid:88) j c τ j,j c = 1 τ (cid:88) j c φ j,j c , (14)where τ j is the total relaxation rate for each subband j ,thus obeying Matthiessen’s rule.Considering I coll j ≈ I (0) j in Eq. (9) already capturesthe intersubband coupling between ρ and ρ and leadsto the Dyakonov-Perel spin relaxation. However, for con-sistency it is important to carry the perturbation expan-sion to the same order in SOC in both sides of Eq. (9).Therefore, in the Appendix B we show that consideringthe SOC corrections in G R ( A ) and in the spectral func-tion, we obtain I (1) j ( k ) ≈ (cid:88) j c τ j,j c (cid:110) H soc j ( k ) , ∂∂ε (cid:104) ρ j c ( ε, θ, R , T ) (cid:105) (cid:111) ε = ε j ( k ) . (15)Now we can proceed with I coll j ≈ I (0) j + I (1) j in Eq. (9) toderive the spin diffusion equation. A. The spin diffusion equation
The Wigner distribution in Eqs. (9), (12), and (15)can be expanded [36] into charge and spin compo-nents as ρ j ( k , R , T ) = g ij ( k , R , T ) σ i and (cid:104) ρ j ( k, R , T ) (cid:105) = (cid:10) g ij ( k, R , T ) (cid:11) σ i (Einstein’s notation for the sums overrepeated indices is implied). Thus g ij and (cid:10) g ij (cid:11) are the σ i component of ρ j for each subband j . Using these,straightforward calculation of the (anti-)commutators al-low us to cast the Boltzmann equation, Eq. (9), into thematrix form (cid:104) K j − eτ (cid:126) F · ∇ k + τ τ j (cid:105) g j = (cid:88) j (cid:48) φ j,j (cid:48) (1 + T j ) (cid:104) g j (cid:48) (cid:105) , (16)where the matrices K j and T j are presented in Ap-pendix C, and the vector g j = ( g j , g xj , g yj , g zj ). Since these are components of the distribution function, integratingover k gives (cid:80) k g j = { N j , S xj , S yj , S zj } , i.e., the chargeand spin densities of each subband j .To obtain a closed set of equations for (cid:104) g j (cid:105) , we firstrearrange Eq. (16) as [34, 36] g j = (cid:88) j (cid:48) φ j,j (cid:48) (cid:104) K j + τ τ j (cid:105) − (cid:104) (1 + T j ) (cid:104) g j (cid:48) (cid:105) (cid:105) + eτ (cid:126) (cid:104) K j + τ τ j (cid:105) − F · ∇ k g j . (17)The inverses can be expressed as the geometric series ex-pansion for ||K j || (cid:28) τ /τ j of order one. Recursivelyreplacing g j on the right hand side up to first order inthe electric field F and integrating over k = ( k, θ ), weobtain (cid:32) φ , φ , φ , φ , (cid:33) ∂ (cid:126) S ∂t = (cid:32) −D − γ γ − φ , φ , D γ − φ , φ , D −D − γ (cid:33) (cid:126) S , (18)where is the 4 × (cid:126) S =[ N , (cid:126)S , N , (cid:126)S ] T is written in a compact notation, with N j ≡ N j ( q , t ), (cid:126)S j ≡ (cid:126)S j ( q , t ), and q = ( q x , q y ) are the co-ordinates of the reciprocal space on the Fourier transform F r → q [ N j ( r , t )] = N j ( q , t ), and equivalently for (cid:126)S j ( q , t ).For simplicity, starting in Eq. (18) and hereafter, we referto the Wigner coordinates ( R , T ) simply as ( r , t ).To express the diffusion matrices for each subband j ,we split it as D j = [ D j q − i q · v j ] + (cid:18) T Θ + Θ D S j (cid:19) , (19)where the spin block D S j contains the spin precession andrelaxation terms, D S j = τ j, − iD j q x Q ∗ j, − τ j, + − iD j q y Q ∗ j, + − iD j q x Q ∗ j, − +2 iD j q y Q ∗ j, + + 1 τ j, − + 1 τ j, + + v j,x Q ∗∗ j, − − v j,y Q ∗∗ j, + − v j,x Q ∗∗ j, − + v j,y Q ∗∗ j, + , (20)while the spin-charge couplings are given by the rect-angular blocks Θ and Θ , which are presented in theAppendix D together with the spin Hall angles. Above, Q ∗ j, ± = Q j, ± − Q j, , Q ∗∗ j, ± = Q j, ± − Q j, and the otherelements are D j = τ j (cid:18) (cid:126) k F j m (cid:19) , (21) v j = eτ j m F , (22) Q j, ± = 2 m (cid:126) λ j, ± = 2 m (cid:126) ( β ,j ± α j ) , (23) Q j, = 2 m (cid:126) β ,j , (24)1 τ j, ∓ = D j (cid:16) Q j, ∓ − Q j, ∓ Q j, + 2 Q j, (cid:17) , (25) γ j = 1 τ j φ , φ j,j , (26)which refer, for each subband j , to the diffusion constant D j , the drift velocity v j = ( v x,j , v y,j ), the linear Q j, ± and cubic Q j, spin-orbit momenta, the Dyakonov-Perelrelaxation rates 1 /τ j, ∓ , and γ j is the intersubband relax-ation rate. B. Strong and weak intersubband coupling regimes
The intersubband coupling is dictated by the φ , over-lap ratio. A particularly interesting scenario arises when φ j,j (cid:48) ≈ j, j (cid:48) ), which allow us to combine theequations for [ N , (cid:126)S ] and [ N , (cid:126)S ] from Eq. (18) as ∂∂t (cid:18) N T (cid:126)S T (cid:19) = −D avg (cid:18) N T (cid:126)S T (cid:19) , (27)where N T = N + N and (cid:126)S T = (cid:126)S + (cid:126)S are the totalcharge and spin densities. In this regime the charge-spin dynamics is led by the subband averaged diffusionmatrix D avg = ( D + D ) /
2. A qualitative picture of thisscenario is clearly seen as a random walk of the particlequickly jumping between subbands [27], thus feeling onlythe subband averaged drift and diffusion forces.In the opposite limit of intersubband coupling φ , ≈ N , (cid:126)S ] and [ N , (cid:126)S ], decouples as Eq.(18) becomes block diagonal. In this case each subbandevolves in time independently [27] as ∂∂t (cid:18) N j (cid:126)S j (cid:19) = −D j (cid:18) N j (cid:126)S j (cid:19) . (28)While the φ , coupling defines the intensity of theintersubband coupling, a more insightful quantity is theintersubband relaxation time t c ≈ φ , − φ , φ , (1 + φ , )( φ , + φ , ) τ . (29)This expression is derived in Appendix E and it charac-terizes the time scale for the initial condition to relaxtowards the subband averaged dynamics given by thestrong coupling regime. Therefore, the weak couplingregime is only dominant for t (cid:28) t c and the strong cou-pling dominates for t (cid:29) t c . As shown in Fig. 2(a), t c /τ vanishes as φ , → φ , . More importantly, for small φ , the t c /τ ≈ / (2 φ , ), which shows that t c /τ (cid:29) φ , (cid:28)
1. Consequently, for typical τ = 1 ps, one would − − φ , − − t c / τ (a) φ , = 1 . φ , = 0 . F z [V/ µ m] . . . . φ j , j (b) φ , φ , φ , t c / τ Figure 2. (a) Intersubband relaxation time t c /τ fromEq. (29) as a function of φ , . For small φ , , t c /τ ≈ / (2 φ , ) shows as a straight line in the log-log scale, andquickly drops to zero as φ , approaches φ , . (b) Densityoverlap ratios φ j,j (cid:48) as a function of F z calculated from theSchr¨odinger-Poisson model with parameters from Table I. Thedashed line refers to the right axis and shows t c /τ calculatedfrom the φ j,j (cid:48) shown in this panel. get t c ≈
100 ps only for φ , ≈ − . Typical values of φ , are much larger than that, as shown in Fig. 2, whichfavors the strong coupling regime. Therefore, we expectthe weak coupling regime to only occur for extreme casesof nearly isolated quantum wells.The Eqs. (27) and (28) above are nearly identical.However, two important differences between the strongand weak coupling regimes are hidden in the details.First, let us assume that the diffusion constants D ≈ D for simplicity. Then, for α ≈ − α the ( Q , ± + Q , ± ) / α ≈ φ = 1 in Eq. (27), which typically shows cir-cular patterns on the spin maps S z ( x, y ) [27, 31], whilefor φ = 0 in Eq. (28) the subbands evolve indepen-dently, each with its α j , which could lead to the persis-tent Skyrmion lattice [49] if the limit t c /τ (cid:29) D j is proportional to τ j , thus D j ∝ / ( φ j, + φ j, ) [see Eq. (21)]. Therefore, for φ , ≈ φ , ≈
1, the D j drops by a factor of 1/2 as φ , is increased from 0 to1. Consequently, the smaller diffusion constant enhancesthe spin lifetime by a factor ∼ IV. DISCUSSION
Hereafter, we are interested in the spin diffusion dy-namics and the spin lifetime. Therefore we may neglectthe small spin-charge coupling in Eq. (19) for simplicity,which reduces the diffusion matrices from 4 × × × × × × × S z,j ( r , t = 0) ≈ δ ( r ) / j = 1 and2), which models the pump-and-probe time-resolved Kerrrotation experiments [20]. In the reciprocal q -space thiscorresponds to an uniform initial S z,j ( q , t = 0). The nu-merical solutions are obtained by propagating Eq. (18) intime independently for each q , and applying the inverseFourier transform to return to r -space.Additionally, in Section IV D we consider the quantumwell from Ref. [49], where it is proposed the persistentSkyrmion lattice (PSL) regime where the subbands areset with orthogonal PSH regimes. This system clearlyshows an insightful transition from strong to weak inter-subband coupling regimes as τ is increased. A. System parameters
Let’s consider a two-subband 2DEG defined by a GaAsquantum well confined along z (cid:107) [001], with x (cid:107) [110] and y (cid:107) [1¯10]. To achieve a two-subbands regime we use awide quantum well with a central barrier as defined inTable I and shown in Fig. 1. Both the central barrierand the lateral regions are composed by Al . Ga . As,
Table I. Typical set of parameters obtained for a wide 38 nmGaAs quantum well with a 0.5 nm central barrier, and totaldensity n D = 8 × cm − . The subband energies ε (0) j ,densities n j , and SOCs α j , β ,j , and β ,j , and overlaps φ j,j (cid:48) are listed for each subband j = { , } at transverse field F z =0. The energy reference is set at the Fermi level ε F = 0.The subband split is ∆ ε = ε (0)2 − ε (0)1 ≈ (cid:126) /τ (cid:28) ∆ ε , thus τ (cid:29) (cid:126) / ∆ ε ≈ . τ = 1 ps. To calculate β ,j = γ D πn j / γ D = 11 meVnm [38].Parameter Value Description n . cm − ] n . α α β , .
16 Linear Dresselhaus couplings [meV nm] β , . β , .
074 Cubic Dresselhaus couplings [meV nm] β , . φ , .
00 Overlap integrals [dimensionless] φ , . φ , . η , .
25 Intersubband Rashbaand Dresselhaus [meV nm]Γ , F z [V/ µ m] . . . α j , β ∗ j [ m e V n m ] (a) + α − α β ∗ β ∗ n j [ c m − ] F z [V/ µ m] . . β , j [ m e V n m ] (b) j = 1 j = 2 Figure 3. Spin-orbit couplings for each subband j = { , } as a function of the electric field F z for the quantum well pa-rameters in Table I. (a) Comparing ± α j with β ∗ j = β ,j − β ,j shows a PSH +1 (blue dot) and two PSH − (red dots) crossings.(b) The cubic Dresselhaus β ,j = γ D πn j / n j (right axis) are approximately constant for small F z and varies linearly for large F z . and symmetric doping is considered to be far away fromthe quantum well. Thus, the Schr¨odinger-Poisson equa-tions [45–48] are solved self-consistently for effective mass m ∗ = 0 . m and dielectric constant κ = 12 .
9, where m is the bare electron mass. To break the structural in-version symmetry and induce a Rasbha SOC we considera transverse electric field F z (cid:107) z and vary its intensity.The SOC coefficients are calculated following Ref. [45–48], yielding the 2D model in Eq. (2).The numerically calculated overlap integrals φ j,j (cid:48) areshown in Fig. 2(b). In the symmetric regime, F z = 0, theenvelope functions ϕ j ( z ) are spread along both quantumwells as symmetric and anti-symmetric solutions, lead-ing to a maximal φ , . As F z increases, the ϕ j ( z ) splitinto left and right quantum wells, as shown in Fig. 1(a),reducing φ , in Fig. 2(b). However, t c only reaches rea-sonably large values for large F z , where φ , ≈ .
03 yields t c ≈ τ . The spin-orbit couplings and subband densi-ties are shown in Fig. 3 as a function of F z . For uncoupledsubbands, one would expect PSH regimes for α j = ± β ∗ j ,which we label as PSH ± j . In Table I we show the in-tersubband SOC parameters η , and Γ , [45–48] as areference, but these are neglected in our model since weconsider non-overlapping Fermi circles, see Fig. 1(b). B. Spin lifetime
The numerical solution of Eq. (18) give us S z ( q , t ), anda Fourier transform yields S z ( r , t ) = F [ S z ( q , t )]. To ana-lyze these numerical data, it is interesting to look at wellknown analytical solutions at particular cases for singlesubband systems and gain insight. Here we consider driftvelocities v j = 0 for simplicity. There are two scenariosthat allow for simple analytical solutions: (i) Q ∗ j, ± = 0;and (ii) Q ∗ j, + = Q ∗ j, − , and τ j, + = τ j, − . The first occursnear PSH regimes [17], and the second is the isotropicregime [31]. As shown in Appendix F, at r = 0, theasymptotic t (cid:29) τ S solutions for these scenarios are S z ( r = 0 , t (cid:29) τ S ) ∝ e − t/τ S t n , (30)where (i) n = 1 and (ii) n = 1 / τ S are shown in Appendix F.Assuming that the numerical solutions of Eq. (18) ap-proximately follow the functional form of Eq. (30), thespin lifetime τ S can be numerically obtained from ∂∂t ln[ S z (0 , t )] = − τ S − nt . (31)Since we assume t (cid:29) τ S , the second term can be ne-glected and the numerical derivative yields τ S indepen-dently of n in the asymptotic limit. Indeed, Fig. 4 showsthat this is a good approximation.To understand the role of the intersubband coupling φ , on the spin lifetime τ S , in Fig. 4 we plot τ S as afunction of F z . In Fig. 4(a) we enforce the weak cou-pling regime by setting φ , ≡
0, while keeping other pa-rameters with their original values from the Schr¨odinger-Poisson simulation. In Fig. 4(b) we restore φ , to itsoriginal values from Fig. 2(b), which leads to the strongcoupling regime. In both Figs. 4(a) and 4(b) the solidblack line is the numerical τ S obtained from Eq. (31), F z [V/ µ m] τ S [ n s ] (a) PSH +1 PSH − PSH − numericalPSH +1 PSH − . . . F z [V/ µ m] . . . τ S [ n s ] (b) PSH +1 numericalPSH +avg t [ns] − t S z ( r = , t ) (c) F z =0.1 V/ µ m F z =0.3 V/ µ m F z =1.6 V/ µ mAnalytical t [ns] − − √ t S z ( r = , t ) (d) Analytical F z =0.0 V/ µ m F z =0.1 V/ µ m F z =0.3 V/ µ m Figure 4. (a)-(b) Spin lifetime τ S as a function of the electricfield F z . In (a) we force φ = 0, while in (b) we use φ (cid:54) = 0from Fig. 2(b). Gray-dotted vertical lines mark the PSH ± j regimes from Fig. 3(a), and the colored dashed lines are takenfrom Eq. (F9). Representative t n S z ( r = 0 , t ) are shown for(c) φ , = 0 and (d) φ , (cid:54) = 0, with a color code matching thepanels above. Dashed lines are analytical solutions from (c)Eq. (F15) and (d) Eq. (F13). The numerical τ S in panels (a)and (b) are taken in the asymptotic t (cid:29) τ S limit, where weexpect for (c) S z ∝ e − t/τ S /t and for (d) S z ∝ e − t/τ S / √ t . while the dashed lines are the analytical expressions de-rived in Appendix F, which generically reads1 τ S ≈ τ ∓ + 1 τ ± − Q ∗± τ ∓ ) D − ( Q ∗± ) D. (32)For the weak coupling regime of Fig. 4(a) we use τ ± → τ j, ± , Q ∗± → Q ∗ j, ± , and D → D j for j = 1 and 2, while inthe strong coupling regime of Fig. 4(b) these parametersare taken as the subband averages, which follows fromthe approximate Eq. (27). In both cases the agreementbetween the numerical τ S and the analytical approxima-tions is patent.Complementary, Figs. 4(c) and 4(d) show tS z (0 , t ) and √ tS z (0 , t ), respectively, for illustrative F z points. InFig. 4(c) all lines match the sum of two independent ex-ponentials [see Eq. (F15)] for the first and second sub-bands, but for simplicity we highlight the analytical solu-tion only for F z = 1 . µ m, which is the one that showsthe largest deviation from a simple exponential. For thestrong coupling we only have approximate analytical so-lutions [Eq. (F13)] for the isotropic regime of F z = 0,which follows from Ref. [31]. In this case, the small de-viation seen in Fig. 4(d) is due to the small imprecisionof τ S in Fig. 4(b).A striking feature seen in Fig. 4(a) is that the τ S peaksare shifted from the PSH conditions ( α j = ± β ∗ j ) fromFig. 3(a). The shift occurs due to the strong dependenceof β ,j upon F z at large fields, as shown in Fig. 3(b).Since the subbands are decoupled for φ , = 0, the spinlifetime of each subband simplifies to τ S ∝ β j, + ( β ∗ j ± α j ) . If β j, were constant, the peak would occur at theexpected PSH condition α j = ± β ∗ j . Indeed for small F z this is approximately true and the peaks in Fig. 4(a)show only a small deviation. However, for large fields thepeak occurs when dτ S /dF z = 0, and the F z dependenceof β j, leads to the shift. This feature is not limited tothe two-subbands system, but can also happen in singlesubband quantum wells if the 2DEG is grounded and thetotal density changes with F z .For realistic φ , (cid:54) = 0, the τ S always shows a peak nearthe isotropic case of F z = 0, as in Fig. 4(b). This is dueto the Matthiessen’s rule, Eq. (14), and the strong depen-dence of φ , with F z shown in Fig. 2. Since 1 /τ j, ± ∝ D j [see Eq. (25)], we find that τ S ∝ (1 + φ , ), which arisesfrom Eq. (14) with φ , = φ , ≈
1. This is a good qual-itative approximation, as seen from Fig. 2(b). Thus, as φ , drops from 1 to 0, τ S falls by a factor of 2 as F z increases. In practice, the peak seen in Fig. 4(b) changesby a factor even larger than 2, which is due to changesin the prefactors beyond this qualitative picture. Never-theless, this shows φ , plays the role of an efficient knobto control the spin lifetime around the isotropic F z = 0configuration and does not require the PSH fine tuningof the SOC. − − q y [ µ m ] (a ) τ S = φ , = 0 −
25 0 25 − y [ µ m ] (a ) t = τ S / −
25 0 25 − (a ) t = 2 τ S − q x [1/ µ m] − q y [ µ m ] (b ) τ S = φ , = 0 −
25 0 25 x [ µ m] − y [ µ m ] (b ) t = τ S / −
25 0 25 x [ µ m] − (b ) t = 2 τ S Figure 5. Comparing the spin maps evolution for (a n ) φ , =0 and (b n ) φ , (cid:54) = 0 in q- ( n = 1) and r-spaces ( n > F z ≈ .
09 V/ µ m (blue dots in Fig. 4). (a n ) Enforcing theweak coupling regime by setting φ , = 0, the S z ( q , t ) showpoles at ± q y ( ± q x ) for the first (second) subbands, whichleads checkboard patterns in (a , a ) S z ( r , t ). (b n ) For a finite φ , the strong coupling regime dominates and the subbandaveraged dynamics show broaden poles in S z ( q , t ) at ± q y ,which gives the nearly striped patterns in (b , b ) S z ( r , t ). C. Spin patterns
Illustrative examples of the spin patterns, S z ( q , t ) and S z ( r , t ), are shown in Figs. 5, 6 and 7 for F z = 0 .
09, 0 . µ m, respectively. These correspond to the blue,orange and purple dots in Figs. 4(a) and 4(b). In thesefigures we compare the spin patterns for φ , = 0 (toppanels) and the realistic φ , (cid:54) = 0 (bottom panels).For F z = 0 .
09 V/ µ m in Fig. 5, the spin lifetimesfor φ , = 0 is similar for both subbands, as shown inFig. 4(a). Moreover, since in this case α and α haveopposite signs, the S z ( q , t ) in Fig. 5(a ) shows four peaksin ± q x and ± q y , corresponding to crossed spin excita-tions, which leads to the checkerboard pattern seen inFigs. 5(a ) and 5(a ). However, for F z = 0 .
09 V/ µ m, theactual subband overlap ratio is φ , ≈ .
9, which gives t c ∼ − ns, and the strong coupling regime dominates.In this case the dynamics follows the subband averageddiffusion matrix, as in Eq. (27). Qualitatively, this leadsto a partial cancellation of D α + D α (since α and α have opposite signs) in the average of the non-diagonalterms of ( D S + D S ) / α j = 0). Since this cancellation is not ex-act, it gives Q ∗ + > Q ∗− and favors S z ( q , t ) peaks in ± q y , asseen in Fig. 5(b ). This leads to the intermediate patternbetween circular and vertical stripes seen in Figs. 5(b )and 5(b ). Numerical experiments show that to recoverthe checkerboard pattern of the weak coupling regime onewould need to reduce φ , by three orders of magnitudeand obtain a large t c = 0 . F z = 0 . µ m shown in Fig. 6 aresimilar to those discussed above for Fig. 5. However, the τ S for φ , = 0 in this case are significantly different for − − q y [ µ m ] (a ) τ S = φ , = 0 −
25 0 25 − y [ µ m ] (a ) t = τ S / −
25 0 25 − (a ) t = 2 τ S − q x [1/ µ m] − q y [ µ m ] (b ) τ S = φ , = 0 −
25 0 25 x [ µ m] − y [ µ m ] (b ) t = τ S / −
25 0 25 x [ µ m] − (b ) t = 2 τ S Figure 6. Equivalent to Fig. 5, but for F z = 0 . µ m (orangedots in Fig. 4). (a n ) For the weak coupling regime ( φ , = 0)the checkboard pattern is transient due to the large differencebetween the subband lifetimes in Fig. 4(a). For large t wesee only the stripes due to the second subband. (b n ) In thestrong coupling regime the subband averaged dynamics leadsto stripped patterns. − − q y [ µ m ] (a ) τ S = φ , = 0 −
25 0 25 − y [ µ m ] (a ) t = τ S / −
25 0 25 − (a ) t = 2 τ S − q x [1/ µ m] − q y [ µ m ] (b ) τ S = φ , = 0 −
25 0 25 x [ µ m] − y [ µ m ] (b ) t = τ S / −
25 0 25 x [ µ m] − (b ) t = 2 τ S Figure 7. Equivalent to Fig. 5, but for F z = 0, which isan isotropic limit due to the symmetric quantum well. In thiscase the circular pattern is present in both (a n ) weak couplingregime with φ , = 0, and (b n ) strong coupling regime with φ , (cid:54) = 0. the first and second subbands, see Fig. (4)(a). Conse-quently, the checkerboard pattern is transient, while forlarge t (cid:29) τ S the second subband spin excitation dom-inates and forms the vertical stripped pattern. Never-theless, for realistic φ , (cid:54) = 0, the strong coupling dom-inates and leads to horizontal stripped pattern seen inFig. 6(b n ).The strong coupling regime is maximum at F z = 0,where φ , peaks. As discussed in the previous section,this φ , peak yields an enhancement of τ S around F z =0. However, since this is the isotropic limit, there’s noquantitative difference between the circular spin patternsin both φ , = 0 and finite φ , (cid:54) = 0 shown in Fig. 7. − .
25 0 .
00 0 . F z [V/ µ m] ∆ E [ m e V ] (b) ¯ h / ∆ E [ p s ] . . . φ j , j (a) φ , φ , φ , t c / τ − .
25 0 .
00 0 . F z [V/ µ m] − . . . . α j , β ∗ j [ m e V n m ] (c) + α − α β ∗ β ∗ Figure 8. Data from the Schr¨odinger-Poisson solution for thequantum well from Ref. [49]. (a) Density overlap ratios φ j,j (cid:48) as a function of the electric field F z and the intersubbandrelaxation time t c /τ ≈ / (2 φ , ) (dashed line, right axes).(b) The subband energy splitting ∆ E = ε − ε (blue line, leftaxes) defines the scale of the momentum scattering time thatsatisfies the non-degeneracy condition τ (cid:29) (cid:126) / ∆ E (dashedline, right axes). (c) The SOC as a function of F z with thePSH conditions β ∗ j = ± α j marked by the circles. In all panels,the vertical dashed gray line mark the point F z = 0 .
084 V/ µ mwhere we consider the PSL condition with α ≈ + β ∗ and α ≈ − β ∗ . D. Persistent Skyrmion Lattice
The checkerboard pattern seen previously in Figs. 5and 6 in the weak coupling limit ( φ , = 0) were pre-dicted in Ref. [49] as a novel topological spin texture,and dubbed persistent Skyrmion lattice (PSL). It ariseswhen the subbands are set with orthogonal PSH regimes, e.g., with α = + β ∗ and α = − β ∗ . Consequently, onesubband shows vertical stripes and the other horizontalones, such that their sum yields the checkerboard pat-tern. However, while the usual single subband PSH isrobust against impurity and electron-electron scattering[17], the two subbands PSL is not robust against inter-subband scattering. To see this, let’s consider the totalHamiltonian composed by H D from Eq. (1) and the pro-jected impurity potential V j (cid:48) j from Eq. (4), which reads H T = (cid:18) H + V V V † H + V (cid:19) , (33)where we have neglected the intersubband SOC H forsimplicity. Now, let’s assume that each single subbandblock commutes with a spin operator as [ H j + V jj , σ j ] = 0.For instance, the PSL conditions α = + β ∗ and α = − β ∗ imply that one subband commutes with σ = σ x and the other is orthogonal with σ = σ y . For a genericspin operator σ T = σ ⊕ σ , one finds[ H T , σ T ] = (cid:18) V ( σ − σ ) V † ( σ − σ ) 0 (cid:19) , (34) where we have used that the V potential is scalar. Con-sequently, [ H T , σ T ] = 0 only if (i) the subbands are inparallel PSH regimes with σ = σ ; or (ii) σ (cid:54) = σ ,but the intersubband impurity coupling V ≈
0. Thefirst case is a trivial superposition of two identical PSHregimes. The second case corresponds to the PSL regimewith σ ⊥ σ . Next, the numerical results of our spin dif-fusion equation show that the PSL checkerboard patternis killed in the strong-coupling regime and it’s recoveredas we increase t c to approach the weak-coupling regime,for which V is indeed negligible.
1. Transition from strong to weak coupling regime
The quantum well discussed in the previous sectionsapproximately satisfy the PSL criteria ( α ≈ + β ∗ and α ≈ − β ∗ ). However, it shows an overall small t c /τ in Fig. 2(b), which favors the strong coupling and thespin pattern becomes nearly circular (Figs. 5 and 6). Tocontrast these results, in this section we now consider thesystem proposed in Ref. [49]. It consists of a 45 nm widequantum well with a wider central barrier of 3 nm, andtotal density n = 4 . × cm − . The data extractedfrom the self-consistent Schr¨odinger-Poisson calculationfor this system are shown in Fig. 8. Here, the widercentral barrier (3 nm instead of 0.5 nm) allows for anenhanced separation of the quantum wells with a reducedoverlap φ , , yielding larger t c /τ , as seen in Fig. 8(a).This helps to favor the weak coupling limit and reachthe PSL regime. However, at F z ≈
0, it also leads tothe small subband energy split shown in Fig. 8(b). Thenon-degenerate regime considered in this paper requires∆ E (cid:29) (cid:126) /τ , which implies τ (cid:29) (cid:126) / ∆ E ≈ . F z ≈
0. For such small ∆ E the subband correlationsneglected by our block-diagonal approximations of H , G ,and Σ might fail. Therefore, hereafter we shall rely onlyon the data for finite F z , such that ∆ E is large enough tojustify the block-diagonal approximations, see Fig. 8(b).From the SOC data in Figs. 8(c) we see that the PSLcondition ( α ≈ β ∗ and α ≈ − β ∗ ) is approximatelysatisfied near F z ∼ .
084 V/ µ m. Around this electricfield we find ∆ E ≈ .
75 meV, which satisfies the non-degenerate condition for a reasonable τ (cid:29) . t c ≈ . τ that is largeenough to allow us to analyse a transition from the strongto the weak coupling regime as we increase τ from 1 psto 10 ps.In Fig. 9(a) we show the spin lifetime τ S multipliedby τ , i.e., τ τ S . This composition is interesting becauseboth in the weak and strong coupling regimes we expect τ S ∝ /τ , due to the DP spin relaxation mechanism.Therefore the product τ τ S becomes an universal quan-tity to characterize the extreme weak and strong couplinglimits. However, in between these extreme regimes, theintersubband couplings γ j breaks this proportionality inEq. (18). Consequently, in Fig. 9(a) we see a transitionfrom the strong to the weak coupling regime as we in-0 − . . . F z [V/ µ m] τ τ S [ p s n s ] (a) τ = 1 ps τ = 5 ps τ = 10 psweakstrongweakstrong . . . /τ [ps − ] . . . . τ S [ n s ] (b) F z [V/ µ m] . (PSL) − . Figure 9. (a) The product τ τ S shows the transition fromthe strong to the weak-coupling regime as τ increases. At F z ≈ t c [see Fig. 8], while the transition emergesfor finite F z due to the larger t c . (b) τ S as a function of 1 /τ for selected F z from panel (a), showing the linear DP trend forboth large and small 1 /τ , with a non-monotonic transitionfor the intermediate regime. In panel (b) the dashed linescorrespond to the strong-coupling limit. In both panels thedashed lines are taken from Eq. (32), where both strong andweak-coupling limits follow the DP relaxation with τ S ∝ /τ .Solid lines are extracted from the numerical solution of thespin diffusion using Eq. (31). crease τ . In Fig. 9(b) we show τ S as a function of 1 /τ for selected values of F z . For F z ≈ t c [seeFig. 8(a)], thus it shows the linear DP trend τ S ∝ /τ over the full range in Fig. 9(b). On the other hand, for F z = − . .
084 V/ µ m we see a non-monotonicevolution of τ S with 1 /τ due to the larger t c for finite F z . For large 1 /τ all cases fall into the DP linear trendof the strong coupling regime (dashed lines). Similarly,for small 1 /τ they approach the weak coupling regime,which also shows the DP linear trend, but it only con-verges to the weak-coupling τ S in for extremely small1 /τ . Notice in Fig. 9(a) that even for τ = 10 ps the τ τ S lines are approaching, but still far from the weak-coupling limit. In between these limits, for intermediate1 /τ , we see the “N”-shaped transition, which is morepronounced for cases with larger t c .The transition from strong to weak coupling regimes isalso clearly seen in the spin density patterns in Fig. 10.These were calculated near the PSL condition for F z ∼ .
084 V/ µ m. First, for τ = 1 ps the system is at thestrong coupling regime with t c = 8 . τ to 5 ps and 10 ps in Figs. 10(b-c), thecharacteristic checkerboard pattern of the PSL regimeemerges as the system approaches the weak couplingregime with t c = 44 ps and 88 ps, respectively. The ex-treme limit of the weak coupling with φ , = 0 is shownin Fig. 10(d) for comparison.This system clearly shows how t c defines the diffu-sion regime between weak and strong subband couplings.However, one must also pay attention to the validity ofdiffusion picture, which requires τ S (cid:29) τ , i.e., the spin −
20 0 20 − y [ µ m ] (a) φ , = 0 , τ = 1 ps −
20 0 20 − (b) φ , = 0 , τ = 5 ps −
20 0 20 x [ µ m] − y [ µ m ] (c) φ , = 0 , τ = 10 ps −
20 0 20 x [ µ m] − (d)weak coupling limit Figure 10. Spin patterns S z ( r , t ) for F z = 0 .
084 V/ µ mcalculated at large t = 2 τ S . As τ increases from (a) 1 ps, (b)5 ps, and (c) 10 ps, the spin dynamics go from (a) the strongcoupling with a nearly circular pattern and approaches (d)the weak coupling regime, where the checkerboard pattern ofthe PSL [49] emerges. For each panel, τ S can be extractedfrom Fig. 9, yielding (a) τ S = 0 .
377 ns, (b) 0.115 ns, and (c)0.137 ns. relaxation must be slow in comparison with the momen-tum scattering time. From Fig. 9, for F z ∼ .
084 V/ µ m,we see that this criterion is certainly valid for τ = 1 psand τ S ≈
377 ps, while for τ = 5 and 10 ps we get τ S ≈
115 ps and 137 ps, respectively. Corresponding tofactors of τ S /τ ∼
23 and 13 .
7, respectively. Therefore,for large τ one might fall towards a (quasi-) ballistic pic-ture and deviations from our model would be expected. V. CONCLUSIONS
We have extended the spin drift-diffusion model tothe case of two-subbands 2DEGs following the kineticequation approach from the Keldysh formalism and semi-classical approximations to achieve the quantum Boltz-mann equation and its linearization towards the diffu-sion equation. Considering the intra- and intersubbandscattering by impurities, we show that the subband dy-namics is controlled by the time scale t c ≈ τ / (2 φ , ),which defines a new knob to control the spin lifetime.We find that for small t c the subbands are strong-coupledand the spin follows a subband-averaged dynamics, whilefor large t c the subbands are nearly independent in theweak-coupling regime. Applying our model to the persis-tent Skyrmion lattice setup [49], we show that its char-acteristic checkerboard pattern is killed by strong inter-subband coupling and reemerges for large t c as it ap-proaches the weak-coupling regime. Moreover, for anynon-degenerate two-subband system we find that the spinlifetime peaks around the symmetric well configurationdue to the Matthiessen’s rule, which follows the corre-1sponding peak in intersubband scattering coupling.Our results apply for non-degenerate subbands, mean-ing that we assume that the broadening of the spec-tral function is small compared to the subband splitting,which justifies the quasi-particle approximation. Thisleads to the block-diagonal approximations in the Green’sfunctions and self-energy. In the opposite limit, for near-degenerate subbands the block-diagonal approximationfail and a model for this scenario shall remain a chal-lenge for future works. VI. ACKNOWLEDGEMENTS
G.J.F. thanks Felix G. G. Hernandez for useful discus-sions. I.R.A. and G.J.F. acknowledge financial supportfrom the Brazilian funding agencies CNPq, CAPES, andFAPEMIG.
Appendix A: Self-energy and impurity self-average
Here we derive the self-energy Σ within the self-consistent Born approximation, considering the weakSOC limit, and neglecting intersubband SOC. We followa similar derivation as in Refs. [24, 54]. However, dueto the two-subband structure, one must consider a basisset by a discrete subband index j = { , } and the usualin-plane plane-wave continuum momentum k = ( k x , k y ),i.e. (cid:104) r | j, k (cid:105) = A − / e i k · r ϕ j ( z ), where A is the normal-ization area for the plane-wave and ϕ j ( z ) is the quantumwell solution of each subband j . Since we will be consid-ering only scalar impurities, we neglect the spin quantumnumber.Within the | j, k (cid:105) basis, the Dyson equation reads G b,a = δ b,a G b + (cid:88) c G b V b,c G c,a , (A1)= δ b,a G b + G b V b,a G a + (cid:88) c G b V b,c G c V c,a G a + · · · such that the random impurity potential matrix elementis V b,c = V j b ,j c ( k b , k c ) = (cid:104) j b , k b | V | j c , k c (cid:105) . As previouslymentioned in the main text, we have introduced a com-pact notation for the indices, e.g., a → j a , k a , σ a .Here we consider random scalar short-range impuritiesset by V ( r ) = (cid:80) i v δ ( r − R i ), where R i = ( x i , y i , z i ) isthe position of each impurity i , and v is the intensity (inunits of energy × volume). The matrix element V b,c fora single set of N random impurities reads V b,c = N (cid:88) i =1 ˜ v j b ,j c ( z i ) e − i ( k b − k c ) · r i , (A2)˜ v j b ,j c ( z i ) = v A ϕ † j b ( z i ) ϕ j c ( z i ) . (A3)To consider an ensemble of random impurities set of iden-tical impurity concentrations n imp = N/ Ω (where Ω is the volume), one defines the ensemble average of anyquantity Q as [24, 54] (cid:104)Q(cid:105) ≡ (cid:90) (cid:104) N (cid:89) i =1 d r i Ω (cid:105) Q ( r , r , ..., r N ) . (A4)To apply the ensemble average on the Dyson equa-tion from Eq. (A1), it is sufficient to consider (cid:104) V b,a (cid:105) and (cid:104) V b,c V c,a (cid:105) . The ensemble average on (cid:104) V b,a (cid:105) straightfor-ward yields (cid:104) V b,a (cid:105) = δ b,a n imp v = , (A5)which is a trivial constant term that can be incorporatedinto the chemical potential. We shall neglect it from nowon. The second-order term, (cid:104) V b,c V c,a (cid:105) reads, (cid:104) V b,c V c,a (cid:105) = N (cid:88) n,p =1 (cid:90) (cid:104) N (cid:89) m =1 d r m Ω (cid:105) × ˜ v j b ,j c ( z n ) e − i ( k b − k c ) · r n ˜ v j c ,j a ( z p ) e − i ( k c − k a ) · r p . (A6)Here one might have two scattering events with the sameimpurity ( n = p ), or with distinct impurities ( n (cid:54) = p ).These two cases have to be considered separately.For n (cid:54) = p the integrals in Eq. (A6) simplify into pairsof products of two integrals equivalent to the one on thefirst-order term, i.e. (cid:104) V b,c V c,a (cid:105) n (cid:54) = p = (cid:104) V b,c (cid:105) N (cid:104) V c,a (cid:105) N − ≈(cid:104) V b,c (cid:105) (cid:104) V c,a (cid:105) . The indexes N and N − N −
1, since we are remov-ing the n = p contribution. Moreover N (cid:29)
1, thus,we can neglect this detail and assume the approximateexpression above. Since this term corresponds to a re-ducible sequence of two first-order scattering events, itcan be neglected as it vanishes by a chemical potentialrenormalization.The n = p contributions from Eq. (A6) simplify into N equivalent integrals, yielding (cid:104) V b,c V c,a (cid:105) n = p = n imp v A Λ b,c,a δ k b , k a , (A7)Λ b,c,a = (cid:90) dzϕ † j b ( z ) | ϕ j c ( z ) | ϕ j a ( z ) . (A8)Later we’ll also use a short notation for Λ a,c ≡ Λ a,c,a .Here, the momentum conservation implies subband con-servation, i.e. δ k b , k a → δ a,b , since we are considering thatthe subbands Fermi circles do not overlap, as shown inFig. 1.Back on the Dyson equation, we get (cid:104) G a (cid:105) = G a + G a (cid:104) n imp v A (cid:88) c Λ a,c G c (cid:105) G a + · · · , (A9)where we have already used (cid:104) G b,a (cid:105) ≈ δ b,a (cid:104) G a (cid:105) . Thisblock-diagonal form is only valid due to the approxima-tion δ k b , k a → δ a,b above.2The term above between square-brackets [ · · · ] is theself-energy in the 1st Born approximation,Σ a = n imp v A (cid:88) c Λ a,c G c = . (A10)Finally, for the self-consistent Born approximation, theusual replacement G c → (cid:104) G c (cid:105) yieldsΣ SCBA a = n imp v A (cid:88) c Λ a,c (cid:104) G c (cid:105) = . (A11) Appendix B: Collision integral
Since H has the block-diagonal form of Eq. (1), hereit is sufficient to analyze each block as H j = ε j ( k ) σ + H soc j ( k ). Consequently, the equilibrium G R ( A ) j are G R ( A ) j ( E, k ) = 1 E − ε j ( k ) σ − H soc j ( k ) ± iη . (B1)The spectral function A j ( E, k ) = i ( G Rj − G Aj ) might takea complicated form due H soc j . However, for small SOC( ε F (cid:29) | H soc j | ), we can consider a series expansion, yield-ing A j ( E, k ) = 2 π exp (cid:104) − H soc j ( k ) ∂∂E (cid:105) δ ( E − ε j ( k )) , (B2) ≈ π (cid:104) δ ( E − ε j ( k )) − H soc k ( k ) δ (cid:48) ( E − ε j ( k )) (cid:105) , where the second expression considers the approximationup to linear order in H soc j . We shall consider this expres-sion for A j ( E, k ) for the quasi-classical approximation.Namely, we take the Kadanoff-Baym ansatz as [23] G Kj ( E, k , R , T ) = − i (cid:110) A j ( E, k ) , ρ j ( k , R , T ) (cid:111) , (B3)where ρ j ( k , R , T ) = 1 − F j ( k , R , T ) is the non-equilibrium distribution function. At equilibrium F j ( k , R , T ) → f j ( k ) becomes the Fermi-Dirac distribu-tion f j ( k ) ≈ (1 + exp[ β ( ε j ( k ) − µ )]) − . Assuming we arenot far from equilibrium, it is reasonable to consider alinear response regime, with ρ j ( k , R , T ) ≈ ρ (0) j ( k ) − N f (cid:48) j ( k ) ρ (1) j ( k , R , T ) , (B4) ≈ ρ (0) j ( k ) + 1 N δ ( ε j ( k ) − ε F ) ρ (1) j ( k , R , T ) , where the second form is taken in the zero temperaturelimit ( β → ∞ , µ → ε F ). Here, the 2DEG density ofstates N = Am/ π (cid:126) is introduced to simplify futureexpressions. The equilibrium contribution ρ (0) j ( k ) van-ishes in Eq. (9), therefore, to focus on the relevant contri-butions, we may write our non-equilibrium distribution function as ρ j ( k , R , T ) ≈ N δ ( ε j ( k ) − ε F ) ρ (cid:48) j ( k , R , T ) . (B5)The spectral function A in Eq. (B2) carries two terms.Therefore, we split the collision integral in Eq. (10) intozeroth and first order terms I coll j ( k ) = I (0) j ( k ) + I (1) j ( k ).The zeroth order term I (0) j ( k ), accounts the contri-butions from the leading order of A j ( E, k ) ≈ πδ ( E − ε j ( k )). As a result both A j and Γ j become scalar matri-ces and the anticommutators simplify, e.g. { A j , Σ Kj } =2 A j Σ Kj . Namely, Γ j = ( n imp /v A ) (cid:80) d Λ j,d A d ( E, k d ) . Therefore, the collision integral in Eq.(12) reads I (0) j ( k ) = 2 πn imp v A (cid:126) (cid:88) d Λ j,d δ ( ε j ( k ) − ε d ( k d )) × (cid:16) ρ d ( k d , R , T ) − ρ j ( k , R , T ) (cid:17) , (B6)where we have expressed G Kj in terms of ρ j ( k , R , T ) usingEq. (B3).The sum in Eq. (B6) can be written as (cid:80) d = (cid:80) j d , k d = N π (cid:80) j d (cid:82) dθ d (cid:82) dε d . The energy conservation implied bythe Dirac delta above enforces ε j ( k ) = ε d ( k d ), whichconstrains the absolute values k and k d . Since ε j ( k )is monotonic in k , we may change variables to cast ρ j ( k , R , T ) → ρ j ( ε j ( k ) , θ, R , T ), leading to I (0) j ( k ) = (cid:88) j d (cid:104) ρ j d ( ε j ( k ) , R , T ) (cid:105) − ρ j ( ε j ( k ) , θ, R , T ) τ j,j d . Finally, I (1) j ( k ) accounts for first order correction in H soc j from Eq. (B2). Following the same procedures ap-plied on I (0) j ( k ), we find I (1) j ( k ) = (cid:88) j d τ j,j d (cid:34)(cid:110) H soc j ( k ) , ∂∂ε (cid:104) ρ j d ( ε, θ, R , T ) (cid:105) (cid:111) − ∂∂ε (cid:68)(cid:110) H soc j d ( ε, θ d ) , ρ j d ( ε, θ d , R , T ) (cid:111)(cid:69) θ d (cid:35) ε = ε j ( k ) , (B7)where we have applied the change of variables k → ( ε j ( k ) , θ ), and on the second term the average is over θ d . Since (cid:104) H soc (cid:105) = 0, the second term yields a negligiblecoupling to the “p-wave” part of ρ j d ( ε, θ d , R , T ). Here-after we neglect this term, yielding I (1) j ( k ) ≈ (cid:88) j d τ j,j d (cid:110) H soc j ( k ) , ∂∂ε (cid:104) ρ j d ( ε, θ, R , T ) (cid:105) (cid:111) ε = ε j ( k ) . (B8)3 Appendix C: Boltzmann equation in the matrix form
Here we present the matrices in Eq. (16). The K j originates by replacing the Hamiltonian in Eq. (2) on the lefthand side of Eq. (9). Similarly, the T j arises by replacing the SOC Hamiltonian in Eq.(15). The K j and T j matricesare given respectively by K j = τ (cid:126) (cid:126) Ω τ iq y λ j, + iq x λ j, − iq y λ j, + (cid:126) Ω τ − k x λ j, − iq x λ j, − (cid:126) Ω τ +2 k y λ j, + k x λ j, − − k y λ j, + (cid:126) Ω τ ++ τ γ D (cid:126) i [2 k x k y q x + q y ( k x − k y )] 2 i [2 k x k y q y + q x ( k y − k x )] 02 i [2 k x k y q x + q y ( k x − k y )] 0 0 4 k x ( k x − k y )2 i [2 k x k y q y + q x ( k y − k x )] 0 0 4 k y ( k x − k y )0 − k x ( k x − k y ) − k y ( k x − k y ) 0 , (C1) T j = λ j, + + 2( k x − k y ) γ D ] k y ( λ j, − − ( k x − k y ) γ D ) k x λ j, + + 2( k x − k y ) γ D ] k y λ j, − − ( k x − k y ) γ D ] k x ∂∂ε j ( k ) , (C2)where Ω = τ ∂ t + i τ (cid:126) m k · q , and γ D is the cubic Dresselhaus SOC coefficient that defines β ,j = γ D k F,j / Appendix D: Spin-charge coupling and the spin Hall angle
Here we present the theta matrices in Eq. (19). These matrices are responsible for the spin-charge coupling andare given by Θ = − iD j q y (cid:104) θ j + − ( Q j, − − Q j, ) + 2 θ j , [ Q j, + − Q j, − ] (cid:105) − iD j q x (cid:104) θ j + − ( Q j, + − Q j, ) + 2 θ j , [ Q j, − − Q j, + ] (cid:105) − i (cid:104) ( θ j + − − θ j , + − θ j , − − θ j , )( v × q ) z − θ j , + − θ j , − )( v x q y + v y q x ) (cid:105) , (D1)Θ = − v y (cid:104) θ j + − ( Q j, − − Q j, ) + 6 θ j , [ Q j, + − Q j, − ] (cid:105) − v x (cid:104) θ j + − ( Q j, + − Q j, ) + 6 θ j , [ Q j, − − Q j, + ] (cid:105) − i ( θ j − θ j − )( v x q y + v y q x ) , (D2)where θ jµ,ν = ( (cid:126) τ j / m ) Q j,µ Q j,ν is the spin-hall angle [36]. Appendix E: Intersubband relaxation rate
To establish a criterion to distinguish between thestrong and weak coupling regimes, we look at the sub-band relaxation rates γ j . In Eq. (18), the γ j enter asdiagonal blocks that dictate the coupling between thesubbands. If we neglect the matrices D j for simplicity, the equation splits into 2 × (cid:32) φ , φ , φ , φ , (cid:33) ∂ (cid:126)F∂t = (cid:18) − γ + γ + γ − γ (cid:19) (cid:126)F , (E1)where (cid:126)F is any pair of coupled charge or spin components,i.e. ( N , N ), or ( (cid:126)S , (cid:126)S ). The general solution is (cid:126)F = c (cid:126)F + c e − t/t c (cid:126)F , (E2)4where (cid:126)F = [1 , T , and (cid:126)F ≈ [1 , − T for φ , ≈
1. Therelaxation time is t c = φ , − φ , φ , (1 + φ , )( φ , + φ , ) τ (E3)For φ , → t c → ∞ and the linear combination of (cid:126)F ± (cid:126)F in Eq. (E2) indicates that the dynamics of the firstand second subbands can be decoupled as in the weakcoupling regime. On the other hand, for φ , → (cid:112) φ , we see that t c →
0, which tell us that the dynamicsquickly falls into the subband-averaged strong couplingregime. Therefore, t c defines how fast the dynamics re-laxes towards the strong coupling regime, as discussed inthe main text. Appendix F: Approximate expressions for the spinrelaxation times
To find analytical solutions for the spin lifetimes, weneglect the small spin-charge coupling and the drift veloc-ity in Eq. (19). This simplifies the spin diffusion equa-tions for the weak and strong coupling, Eqs. (27)-(28),which take the 3 × ∂ (cid:126) S ∂t = −D (cid:126) S , (F1)where D = Dq + D S , with D S = τ − iDQ ∗− q x τ + − iDQ ∗ + q y − iDQ ∗− q x +2 iDQ ∗ + q y τ − + 1 τ + . (F2)Notice that here these are written in a generic form.Namely, for the weak subband regime one must considerthe independent subbands by adding the j index to thequantities above: (cid:126) S j , D j , Q ∗ j, ± , τ j, ± . On the other hand,for the strong coupling regime, D must be taken as thesubband averaged diffusion matrix, yielding (cid:126) S → (cid:126) S + (cid:126) S , (F3) D → D + D , (F4) Q ∗± → D Q ∗ , ± + D Q ∗ , ± D + D , (F5)1 τ ± → (cid:18) τ , ± + 1 τ , ± (cid:19) . (F6)In this section we proceed with the generic form ofEqs. (F1)-(F2).A generic solution of Eq. (F1) can be written in termsof its eigenvalues as (cid:126) S ( q , t ) = (cid:88) n c n e − ω n t φ n , (F7) where ω n ≡ ω n ( q ) and φ n ≡ φ n ( q ) are the eigenvaluesand eigenvectors of D , and c n ≡ c n ( q ) are the linearcombination coefficients to be set by the initial condi-tions. For the generic matrix D in Eq. (F2) the eigenval-ues/eigenvectors expressions are cumbersome, but simpleanalytical solutions exist for the cases where: (i) Q ∗− = 0or Q ∗ + = 0; and (ii) Q ∗− = Q ∗ + ≡ Q ∗ and τ + = τ − ≡ τ .These solutions are well known [17, 31], and below wesimply revise and write them in a convenient and genericform.To obtain the spin lifetimes in each of these casesabove, we seek for the the q that extremizes Re[ ω n ( q )].For instance, let’s consider the first case of Q ∗− = 0 andfocus on the decoupled 2 × D . Its eigenvaluesare ω ± ( q ) = − τ − − τ + − Dq ± (cid:112) Q ∗ + τ − Dq y ) τ − , (F8)which are maximum (least negative) at q = (0 , ± q y ),with q y = (cid:112) Q ∗ + ) ( Dτ − ) − / (4 Q ∗ + Dτ − ). Fromthese, the spin lifetime is defined as τ S = − /ω + ( q ).Namely,1 τ ( i ) S ≈ τ ∓ + 1 τ ± − Q ∗± τ ∓ ) D − ( Q ∗± ) D, (F9)where the top and bottom signs refer to the cases of Q ∗− =0 or Q ∗ + = 0, respectively. Similarly, for the case (ii) wefind 1 τ ( ii ) S ≈ τ − Q ∗ τ ) D − ( Q ∗ ) D. (F10)Notice that the expression for τ ( i ) S reduces to τ ( ii ) S in theappropriate limit of Q ∗− = Q ∗ + ≡ Q ∗ and τ + = τ − ≡ τ .Therefore, it is sufficient to consider only Eq. (F9) as ageneric expression for the spin lifetime that incorporatesboth cases.In the single subband case, or in the weak couplingregime, case (i) corresponds to the PSH regime and τ (i) S simplifies to1 τ PSH ± j S = 2 m D j (cid:126) (cid:104) β j, + ( β ∗ j ± α j ) (cid:105) . (F11)
1. Expressions for the spin dynamics
The full expressions for S z ( q , t ) or its Fourier transform S z ( r , t ) = F [ S z ( q , t )] are also well known [17, 31]. Forthe case (i) with Q ∗− = 0, the solution in q-space have twopoles at q = ± (0 , q y ), which leads to the stripped patternin the S z ( r , t ) as these poles form a cosine function in theFourier transform, which reads S z ( r , t ) ≈ e − r Γ e − t/τ S cos( κy ) , (F12)5where the diffusion broadening Γ = √ Dt , the spin life-time τ S is given by Eq. (F9), and κ ≈ Q ∗ + . Similarly, for Q ∗ + = 0 the solution follows from the one above replacingthe last term by cos( κx ), with κ ≈ Q ∗− .For the isotropic limit of case (ii) above, with Q ∗− = Q ∗ + ≡ Q ∗ and τ + = τ − ≡ τ , an approximate solutionfor S z ( r , t ) in the asymptotic limit of large times t (cid:29) τ S is shown in Ref. [31] (see their Eqs. 28-30). It’s alsopossible to obtain an exact solution for S z ( r = 0 , t ) bydirect integration of their Eq. 29 with r = 0, which reads S z (0 , t ) = e − t/τ + e − t/τ πt/τ + ζ ( t ) e − t/τ S (cid:112) πt/τ , (F13) ζ ( t ) = ( γ − γ (cid:88) ± erf (cid:32) (cid:0) γ ± (cid:1) γ (cid:114) tτ (cid:33) , (F14)where τ S ≡ τ ( ii ) S from Eq. (F10), and γ = 2 Q ∗ √ Dτ .Notice that at r = 0 the denominator of Eq. (F12) isproportional to t , while in Eq. (F13) it has t and √ t con-tributions at short times, and asymptotically approaches √ t as the second term dominates for large t . These fea-tures are shown and discussed in Fig. 4 in the main text.Particularly, for the weak coupling case of Fig. 4(c), thesolution is the sum of the contributions for each subband,each in the form of Eq. (F12), which at r = 0 reads S z (0 , t ) = e − t/τ D t + e − t/τ D t . (F15) [1] S. Datta and B. Das, Electronic analog of the electro-optic modulator, Appl. Phys. Lett. , 665 (1990).[2] P. Chuang, S.-C. Ho, L. W. Smith, F. Sfigakis, M. Pep-per, C.-H. Chen, J.-C. Fan, J. P. Griffiths, I. Farrer, H. E.Beere, G. A. C. Jones, D. A. Ritchie, and T.-M. Chen,All-electric all-semiconductor spin field-effect transistors,Nat. Nanotechnol. , 35 (2014).[3] J. M. Kikkawa and D. D. Awschalom, Lateral drag of spincoherence in gallium arsenide, Nature , 139 (1999).[4] S. A. Wolf, Spintronics: A spin-based electronics visionfor the future, Science , 1488 (2001).[5] I. ˇZuti´c, J. Fabian, and S. Das Sarma, Spintronics: Fun-damentals and applications, Rev. Mod. Phys. , 323(2004).[6] J. Fabian, A. Matos-Abiague, C. Ertler, P. Stano, andI. ˇZuti´c, Semiconductor spintronics, Acta Phys. Slovaca , 565 (2007).[7] D. D. Awschalom and M. E. Flatt´e, Challenges for semi-conductor spintronics, Nat. Phys. , 153 (2007).[8] M. Wu, J. Jiang, and M. Weng, Spin dynamics in semi-conductors, Phys. Rep. , 61 (2010).[9] A. Manchon, H. C. Koo, J. Nitta, S. M. Frolov, and R. A.Duine, New perspectives for Rashba spin–orbit coupling,Nat. Mater. , 871 (2015).[10] R. Winkler, Spin–Orbit Coupling Effects in Two-Dimensional Electron and Hole Systems (Springer BerlinHeidelberg, 2003).[11] G. Dresselhaus, Spin-Orbit Coupling Effects in ZincBlende Structures, Phys. Rev. , 580 (1955).[12] Y. A. Bychkov and E. I. Rashba, Oscillatory effects andthe magnetic susceptibility of carriers in inversion layers,J. Phys. C , 6039 (1984).[13] M. I. D’yakonov and V. I. Perel’, Possibility of orientingelectron spins with current, JETP Lett. , 467 (1971).[14] R. J. Elliott, Theory of the effect of spin-orbit coupling onmagnetic resonance in some semiconductors, Phys. Rev. , 266 (1954).[15] Y. Yafet, g factors and spin-lattice relaxation of conduc-tion electrons, in Solid State Physics (Elsevier, 1963) pp.1–98.[16] J. Schliemann, J. C. Egues, and D. Loss, Nonballisticspin-field-effect transistor, Phys. Rev. Lett. , 146801 (2003).[17] B. A. Bernevig, J. Orenstein, and S.-C. Zhang, ExactSU(2) Symmetry and Persistent Spin Helix in a Spin-Orbit Coupled System, Phys. Rev. Lett. , 236601(2006).[18] J. D. Koralek, C. P. Weber, J. Orenstein, B. A. Bernevig,S.-C. Zhang, S. Mack, and D. Awschalom, Emergence ofthe persistent spin helix in semiconductor quantum wells,Nature , 610 (2009).[19] C. P. Weber, J. Orenstein, B. A. Bernevig, S.-C. Zhang,J. Stephens, and D. D. Awschalom, Nondiffusive spindynamics in a two-dimensional electron gas, Phys. Rev.Lett. , 076604 (2007).[20] M. Walser, C. Reichl, W. Wegscheider, and G. Salis, Di-rect mapping of the formation of a persistent spin helix,Nat. Phys. , 757 (2012).[21] J. Ishihara, Y. Ohno, and H. Ohno, Direct imag-ing of gate-controlled persistent spin helix state in amodulation-doped GaAs/AlGaAs quantum well, Appl.Phys. Express , 013001 (2013).[22] E. G. Mishchenko, A. V. Shytov, and B. I. Halperin, Spincurrent and polarization in impure two-dimensional elec-tron systems with spin-orbit coupling, Phys. Rev. Lett. , 226602 (2004).[23] J. Rammer and H. Smith, Quantum field-theoreticalmethods in transport theory of metals, Rev. Mod. Phys. , 323 (1986).[24] J. Rammer, Quantum field theory of non-equilibriumstates (Cambridge University Press, 2007).[25] H. Haug and A.-P. Jauho,