%From 19Feb14



%\documentclass[a4paper,10pt]{article}
\documentclass[aps,prb,preprintnumbers,amsmath,amssymb]{revtex4}
%\usepackage{graphicx}% Include figure files
\usepackage{dcolumn}% Align table columns on decimal point
\usepackage{bm}% bold math
\usepackage{amsmath}
\usepackage{feynmp}

\bibliographystyle{apsrev}

\usepackage{graphicx}% Include figure files

\begin{document}

\title{Supplementary Material: Hole Spin Helix: Anomalous Spin Diffusion in Anisotropic Strained Hole Quantum Wells}

\author{Vincent E. Sacksteder IV$^{1}$}
\email{vincent@sacksteder.com}
\affiliation{Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China}
\affiliation{Division of Physics and Applied Physics, Nanyang Technological University, 21 Nanyang Link, Singapore 637371}


\author{B. Andrei Bernevig}
 \affiliation{Princeton Center for Theoretical
Science, Princeton, NJ 08544} \affiliation{$^2$Department of
Physics, Princeton University, Princeton, NJ 08544}

 %\pacs{pacs numbers}


\date{\today}% It is always \today, today,
             %  but any date may be explicitly specified

\maketitle


%compound, r in e A A, b in e V A A A, b/r in V A
%gaas, -14.62, -81.93, 5
%alas, -1.501, -33.51, 20
%inas -159.9, -50.18, 0.3
%insb -548.5, -934.8, 2
%cdte -10.79, -76.93, 7
%znse -4.099, -62.33, 16


%b kz kz =pm  E r

%b/r = pm E / (kz kz) = pm E a a = pm (delta V) a

%They use cubic terms in the 8v8v Hamiltonian that go beyond the Luttinger Hamiltonian.  One term labeled by b is caused by structural inversion asymmetry and seems to be independent of the splitting Delta.  The other term labeled by r is both caused by and proportional to the external electric field E.  

%lambdaD = 3 b/ (4 gamma2). lambdaR = 3 E r / (4 gamma2 kz kz).  These both get multiplied by kf kf kf. 

%compare to alpha =  6 e E (gamma3/ 4 gamma2 kz kz )^2 = (3/8) e E (a/2pi)^4 (gamma3/  gamma2 )^2

%So we need to compare e E (a/2pi)^4 to E r (a/2pi)^2,  1 to (r/e)/ (a/2pi)^2, 1 to (r/e) kz^2, the second wins if if a is too small and the critical (a/2pi) ranges from 2 A in znse to 25 A in insb.

%Also we need to compare e E (a/2pi)^4 to b, the second wins if a is too small and the critical (a/2 pi) ranges from 2 A / (E)^{1/4} in alas to 6 A / (E)^{1/4} in InSb.  No, one should compare (Delta V) (a/2pi)^3   to 2pi b/e.    I think that the critical values of (a/2pi) range from 6 to 14  A / (Delta V)^{1/3} but I wasn't careful.



\section{Approximate Form of the Heavy Hole Effective Hamiltonian}
The four-band Luttinger model \cite{Luttinger56, Winkler03, Lu05}, with the Bir-Pikus strain Hamiltonian, in a quantum well extended along the $x - y$ plane and grown along the $001$ ($z$) axis, is represented in the $J_z = +3/2, -3/2, +1/2, -1/2$ basis as:
\begin{eqnarray}
V(z) & + &    \begin{bmatrix} k_{HH} & 0 & S & R^*  \\ 0 & k_{HH}  & R & -S^*  \\ S^* & R^* &  k_{LL}  & 0  \\ R  & -S & 0 & k_{LL}  \end{bmatrix}
\nonumber \\
k_{HH} &=& (\gamma_1 + \gamma_2)k^2/2m + (\gamma_1 - 2 \gamma_2)k_z^2/2m  -a (\epsilon_{xx} + \epsilon_{yy} + \epsilon_{zz}) +b(\epsilon_{zz} - \epsilon_{xx}/2 - \epsilon_{yy}/2)
\nonumber \\
k_{LL} & = & (\gamma_1 - \gamma_2)k^2/2m + (\gamma_1 + 2 \gamma_2)k_z^2/2m   -a (\epsilon_{xx} + \epsilon_{yy} + \epsilon_{zz}) -b(\epsilon_{zz} - \epsilon_{xx}/2 - \epsilon_{yy}/2)
\nonumber \\
S & = & -  \sqrt{3} \gamma_3 k_- k_z/ m  -  d ( \imath \epsilon_{zy} -  \epsilon_{zx})
 \nonumber \\
R & = &   \sqrt{3} \gamma_2 (k_y^2 - k_x^2) / 2m - \imath \sqrt{3} \gamma_3 k_x k_y / m  + \sqrt{3}b (\epsilon_{xx} - \epsilon_{yy}) /2+ \imath d \epsilon_{xy}
 \nonumber \\
 & = &   -\frac{\sqrt{3}}{4m} (k_+^2 (\gamma_2 + \gamma_3) + k_-^2 (\gamma_2 - \gamma_3) -  2 \gamma_3 \beta^2 e^{\imath 2 \theta}), \;  \beta^2 e^{\imath 2 \theta} = \frac{ m b}{ \gamma_3} (\epsilon_{xx} - \epsilon_{yy})   + \imath  \frac{2 m d} {\sqrt{3} \gamma_3 }\epsilon_{xy}  
 \nonumber \\
 k_+ &=& k_x + \imath k_y, \; k_- = k_x - \imath k_y
\end{eqnarray}
$m$ is the electron mass, $\gamma_1, \gamma_2, \gamma_3$ are the material-specific Luttinger Hamiltonian parameters, $a,b,d$ are the material-specific strain deformation potentials, and the $\epsilon$ parameters describe the strain on the sample.  $V(z) = V_c + V_E $ is the sum of the quantum well's confinement potential $V_c(z)$ (which is symmetric under inversions) and a symmetry-breaking term $ V_E  = - e E z$.   $k_{HH}$ and $k_{LL}$  are kinetic operators for the heavy and light holes respectively.  $S$ couples  states whose spin have the same sign, while $R$ couples states whose spin have opposite sign.

%\begin{eqnarray}
%H^r_{8v 8v} &=& r^{8v 8v}_{41} E \begin{bmatrix} 0 & \frac{\imath \sqrt{3}}{2} k_- & 0  & 0 \\  -\frac{\imath \sqrt{3}}{2} k_+ & 0 & \imath k_- & 0 \\  0 & - \imath k_+ & 0 &  \frac{\imath \sqrt{3}}{2} k_- \\ 0 & 0 & - \frac{\imath \sqrt{3}}{2} k_+ & 0 \end{bmatrix}
%\\ \nonumber
%H^r_{8v 8v} &=& r^{8v 8v}_{41} E \begin{bmatrix} 0 & 0 & \frac{\imath \sqrt{3}}{2} k_- & 0   \\ 0 & 0 & 0 & - \frac{\imath \sqrt{3}}{2} k_+  \\  -\frac{\imath \sqrt{3}}{2} k_+ & 0  & 0 & \imath k_- \\  0 &  \frac{\imath \sqrt{3}}{2} k_-  & - \imath k_+ & 0  \end{bmatrix}
%\\ \nonumber
%H^b_{8v 8v} &=& b^{8v 8v}_{41} (k_x (k_y^2 - k_z^2) J_x + k_y (k_z^2 - k_x^2) J_y + k_z (k_x^2 - k_y^2) J_z)
%\nonumber \\
 %&=& b^{8v 8v}_{41} (k_x (k_y^2 - k_z^2) /2 \begin{bmatrix} 0 & \sqrt{3} & 0 & 0 \\  \sqrt{3}  & 0 & 2 & 0 \\ 0 & 2 &  0 & \sqrt{3} \\ 0 & 0 & \sqrt{3}  & 0 \end{bmatrix} 
 %\nonumber \\ 
% &+& \imath k_y (k_z^2 - k_x^2) /2 \begin{bmatrix} 0 & -\sqrt{3} & 0 & 0 \\  \sqrt{3}  & 0 & -2 & 0 \\ 0 & 2 &  0 & -\sqrt{3} \\ 0 & 0 & \sqrt{3}  & 0 \end{bmatrix}  
 %\nonumber \\ 
 %&+& k_z (k_x^2 - k_y^2) / 2 \begin{bmatrix} 3 & 0 & 0 & 0 \\  0  & 1 & 0 & 0 \\ 0 & 0 & -1 & 0 \\ 0 & 0 &0   & -3 \end{bmatrix} )
%\nonumber \\
% &=& b^{8v 8v}_{41}/2 ( \begin{bmatrix} 0 & \sqrt{3} k_x (k_y^2 - k_z^2) & 0 & 0 \\  \sqrt{3}  k_x (k_y^2 - k_z^2)  & 0 & 2k_x (k_y^2 - k_z^2)  & 0 \\ 0 & 2 k_x (k_y^2 - k_z^2)  &  0 & \sqrt{3}  k_x (k_y^2 - k_z^2) \\ 0 & 0 & \sqrt{3} k_x (k_y^2 - k_z^2)  & 0 \end{bmatrix} 
 %\nonumber \\ 
% &+&  \begin{bmatrix} 0 & -\sqrt{3} \imath k_y (k_z^2 - k_x^2) & 0 & 0 \\  \sqrt{3} \imath k_y (k_z^2 - k_x^2)  & 0 & -2 \imath k_y (k_z^2 - k_x^2) & 0 \\ 0 & 2 \imath k_y (k_z^2 - k_x^2) &  0 & -\sqrt{3} \imath k_y (k_z^2 - k_x^2) \\ 0 & 0 & \sqrt{3} \imath k_y (k_z^2 - k_x^2)  & 0 \end{bmatrix}  
 %\nonumber \\ 
% &+&   \begin{bmatrix} 3 k_z (k_x^2 - k_y^2) & 0 & 0 & 0 \\  0  & k_z (k_x^2 - k_y^2) & 0 & 0 \\ 0 & 0 & -k_z (k_x^2 - k_y^2) & 0 \\ 0 & 0 &0   & -3 k_z (k_x^2 - k_y^2) \end{bmatrix} )
%\nonumber \\
% &=& ( \begin{bmatrix} 0 & \sqrt{3} k_x (k_y^2 - k_z^2)   -\sqrt{3} \imath k_y (k_z^2 - k_x^2) & 0 & 0 \\  \sqrt{3}  k_x (k_y^2 - k_z^2)  + \sqrt{3} \imath k_y (k_z^2 - k_x^2)& 0 & 2k_x (k_y^2 - k_z^2)  -2 \imath k_y (k_z^2 - k_x^2)  & 0 \\ 0 & 2 k_x (k_y^2 - k_z^2)  +  2 \imath k_y (k_z^2 - k_x^2)  &  0 & \sqrt{3}  k_x (k_y^2 - k_z^2)   -\sqrt{3} \imath k_y (k_z^2 - k_x^2)  \\ 0 & 0 & \sqrt{3} k_x (k_y^2 - k_z^2)  +  \sqrt{3} \imath k_y (k_z^2 - k_x^2)& 0 \end{bmatrix} 
 %\nonumber \\ 
 %&+&   \begin{bmatrix} 3 k_z (k_x^2 - k_y^2) & 0 & 0 & 0 \\  0  & k_z (k_x^2 - k_y^2) & 0 & 0 \\ 0 & 0 & -k_z (k_x^2 - k_y^2) & 0 \\ 0 & 0 &0   & -3 k_z (k_x^2 - k_y^2) \end{bmatrix} )b^{8v 8v}_{41}/2 
%\nonumber \\
 %&=& ( \begin{bmatrix} 0  & 0 & \sqrt{3} k_x (k_y^2 - k_z^2)   -\sqrt{3} \imath k_y (k_z^2 - k_x^2) & 0    \\ 0 & 0 & 0 & \sqrt{3} k_x (k_y^2 - k_z^2)  +  \sqrt{3} \imath k_y (k_z^2 - k_x^2)\\  \sqrt{3}  k_x (k_y^2 - k_z^2)  + \sqrt{3} \imath k_y (k_z^2 - k_x^2) & 0 & 0 & 2k_x (k_y^2 - k_z^2)  -2 \imath k_y (k_z^2 - k_x^2)   \\ 0 & \sqrt{3}  k_x (k_y^2 - k_z^2)   -\sqrt{3} \imath k_y (k_z^2 - k_x^2) & 2 k_x (k_y^2 - k_z^2)  +  2 \imath k_y (k_z^2 - k_x^2)  &  0   \end{bmatrix} 
% \nonumber \\ 
% &+&   \begin{bmatrix} 3 k_z (k_x^2 - k_y^2) & 0  & 0 & 0 \\ 0 & -3 k_z (k_x^2 - k_y^2)  & 0 &0     \\  0  & 0 & k_z (k_x^2 - k_y^2) & 0  \\ 0 & 0  & 0 & -k_z (k_x^2 - k_y^2) \end{bmatrix} )b^{8v 8v}_{41}/2 
%\end{eqnarray}
%Maybe the above terms should be multiplied by $-1$ because I also multiplied Winkler's version of the Luttinger hamiltonian by $-1$.  Winkler's Dresselhaus term must come from $H^r_{8v 8v} $ multiplied by the $R$ matrix element and by $G_{LL}$.  $H^b_{8v 8v} \propto b^{8v 8v}_{41} k_x k_z^2$ so  Winkler's Rashba term must come from $H^b_{8v 8v} $ multiplied by the $R$ matrix element and by $G_{LL}$.  So both of Winkler's terms are propotional to $R$.   We were at second order because the only way to get splitting was by interaction of $S$'s $k_z$ with the light hole  potential $V(z)$ which is small and therefore has only a second order effect on the light holes.

Hooke's law implies that in-plane strain $\epsilon_{xx}, \epsilon_{yy}, \epsilon_{xy}$ applied to the quantum well will cause out-of-plane strain as well. \footnote{Many thanks to Roland Winkler for explaining these points about Hooke's law.}  When the crystal growth is along the 001 direction this physics simplifies and produces only one extra strain component, $\epsilon_{zz}$.  The only effect of $\epsilon_{zz}$ is a shift in the splitting between the light and heavy holes.  For other growth directions $\epsilon_{zx}, \epsilon_{zy}$ contribute to the $S$ matrix element, but we will show that these terms make no contribution to the commutator which controls the spin orbit interaction.

Since we are only interested in observables constructed from heavy holes, we can apply a unitary transformation to the light holes: $ | +1/2 \rangle  \rightarrow (k_+ / k )  | +1/2 \rangle, \;    |  -1/2 \rangle  \rightarrow (k_- / k )  |  -1/2 \rangle$.  As a result, the Hamiltonian transforms:
$ S \rightarrow  (k_+ / k )  S=- \sqrt{3} \gamma_3 k k_z/ m  -   (k_+ / k )d (\imath \epsilon_{zy} - \epsilon_{zx}), \  R \rightarrow  (k_+ / k ) R =  -\frac{\sqrt{3}}{4m} ((k_+^3 / k) (\gamma_2 + \gamma_3) + k_- k (\gamma_2 - \gamma_3) -( k_+ /k) 2 \gamma_3  \beta^2 e^{\imath 2 \theta})$.  


The second and third terms in $R$ combine to have constant phase when $\beta^2 =  k_F^2 (1 - \gamma_2 / \gamma_3) $, and the first term in $R$ can be completely eliminated by setting $\gamma_2 / \gamma_3 = -1$, rendering $R$ completely real.  In momentum space the only remaining complex term in the Luttinger Hamiltonian is  $S$'s strain term $- \imath  (k_+ / k ) d  \epsilon_{zy}$.  

We will now prove that  $S$'s strain term $- \imath  (k_+ / k ) d ( \epsilon_{zy}+ \imath \epsilon_{zx}) $ has no effect on the coupling between the heavy holes.   The effective Hamiltonian for the heavy holes is:
\begin{eqnarray}
H_{HH}&=& V(z) +  \begin{bmatrix} k_{HH}  & 0 \\ 0 & K_{HH} \end{bmatrix} - \begin{bmatrix} S  & R^* \\ R & -S^* \end{bmatrix} (k_{LL} + V(z) - E)^{-1} \begin{bmatrix} S^*  & R^* \\ R & -S \end{bmatrix}
\nonumber \\
& = &  V(z)  +   \begin{bmatrix} \acute{k}_{HH}  &   [ (k_{LL} + V(z) - E)^{-1}, S] R^* 
 \\   - [(k_{LL} + V(z) - E)^{-1} , S^*]  R & \acute{K}_{HH}
\end{bmatrix}
\nonumber \\
\acute{k}_{HH} &= & k_{HH} - S (k_{LL} + V(z) - E)^{-1} S^* - R^* (k_{LL} + V(z) - E)^{-1} R
\end{eqnarray}
The strain term $- \imath  (k_+ / k ) d ( \epsilon_{zy}+ \imath \epsilon_{zx})$ has no $z$ dependence, so it commutes with $k_{LL}$ and $V(z)$, and makes no contribution to the commutator.   In fact the commutator reduces to  $- \sqrt{3} \gamma_3 k / m [ (k_{LL} + V(z) - E)^{-1}, k_z]$.  

After transforming  to position space the commutator is $-\imath \sqrt{3} \gamma_3 k / m [ (k_{LL} + V(z) - E)^{-1}, \partial_z]$.  Its phase is manifestly constant.   This proves that  the off-diagonal  elements of the heavy hole Hamiltonian $H_{HH}$  have (up to a factor of $\imath$) the same phase as $R$.    If $R$'s phase is constant then the spin-orbit interaction also has constant phase and one component of the  spin is conserved.  This is an exact result for the full four-band Luttinger Hamiltonian.

Moreover we note that diagonal elements of $H_{HH}$ are proportional to the identity; any splitting between the heavy holes comes only from the off-diagonal elements.


The previous steps were exact.  Now we obtain a simple but approximate form for the heavy-hole Hamiltonian.  We  simplify the effective Hamiltonian $H_{HH}$ by assuming that the quantum well is thin, and that the confinement potential $V_c(z)$ along the $z$ axis  splits the spectrum into 2-D sub bands corresponding to eigenstates of $V_c$.   We assume that the Fermi level is used to regulate the  hole carrier density to a small value where only the first 2-D sub-band contributes to transport.  Therefore we can assume that the system is in the lowest eigenstate of $V_c(z)$, which can be be replaced everywhere by its lowest eigenvalue.  

 Since the quantum well is thin, $ k_z \gg k_x, k_y$ and therefore $R$'s contribution to the diagonal is negligible.  Next we assume that the charge asymmetry is small compared to the splitting between light and heavy holes  ($V_E \ll k_{LL} - E$), we expand in powers of $V_E$, and we evaluate the expectation value $\langle k_{LL} - E \rangle =  \Delta E$. 
\begin{eqnarray}
H_{HH}&=&  V(z)  +   \begin{bmatrix} \acute{k}_{HH} - S S^* (\Delta E)^{-1}   & -\imath \alpha  \frac{2m}{\sqrt{3}  \gamma_3}  k  R^*
 \\    \imath  \alpha \frac{2m}{\sqrt{3}  \gamma_3}  k  R & \acute{k}_{HH} - S S^* (\Delta E)^{-1} 
\end{bmatrix}
\nonumber  \\
 \alpha  &=&  \imath 6 \gamma_3^2 [V_E, k_z]/(2 m \Delta E)^2, \; [V_E, k_z ] = -\imath E_z 
\end{eqnarray}
$\alpha$ is the strength of the heavy hole spin-orbit interaction.  Lastly we approximate the diagonal elements of the effective Hamiltonian by  ignoring terms which are independent of $k_+, k_-$ and finding the effective mass $m_H$ which controls the diagonal's in-plane kinetic energy $k^2 / 2 m_H$.  This requires calculation of the expectation value of $ \langle k_z^2 \rangle$, which we perform using  $S S^* \approx 3 \gamma_3^2  k^2 \langle k_z^2 \rangle / m^2$.   We  use $ \Delta E = \langle k_{LH} - k_{HH} \rangle \approx \frac{ 2 \gamma_2}{m} \langle k_z^2 \rangle  $ and obtain the final effective Hamiltonian of the heavy holes:
\begin{eqnarray}
H_{HH}&=&   \begin{bmatrix} k^2 / 2 m_H  & -\imath \alpha   \frac{2m}{\sqrt{3}  \gamma_3}  k  R^*
 \\    \imath \alpha   \frac{2m}{\sqrt{3}  \gamma_3}  k  R & k^2 / 2 m_H
\end{bmatrix}
\\ \nonumber 
&=&  \begin{bmatrix} k^2 / 2 m_H  & \imath \alpha  (k_-^3 (1+\gamma_2/  \gamma_3 )/2 - k^2 k_+ (1-\gamma_2 / \gamma_3 )/2 -  k_- \beta^2 e^{-\imath 2 \theta})
 \\  - \imath \alpha  (k_+^3 (1+\gamma_2/  \gamma_3 )/2 - k^2 k_- (1-\gamma_2 / \gamma_3 )/2 -  k_+ \beta^2 e^{\imath 2 \theta})& k^2 / 2 m_H
\end{bmatrix}
\nonumber \\
 &=& k^2/2 m_H + a_x \sigma_x + a_y \sigma_y, \; 
 \nonumber \\
  a_x   &=& - \alpha  k_y(- 2 k_x^2  (1+\gamma_2/  \gamma_3 ) + k^2   (\gamma_2/  \gamma_3 )) -  \alpha  \beta^2  (k_x \sin 2 \theta + k_y \cos 2 \theta)
 \nonumber \\
            a_y    &=&  -  \alpha  k_x(- 2  k_y^2 (1+\gamma_2/  \gamma_3 ) + k^2  (\gamma_2/  \gamma_3 ) )+  \alpha  \beta^2  (k_x \cos 2 \theta - k_y \sin 2 \theta)
\end{eqnarray}
The  renormalized in-plane mass is $ m_H= m/(\gamma_1 + \gamma_2 - 3 \gamma_3^2/\gamma_2)$.




\section{Derivation of The Spin Diffusion Equations}


We study the equations of motion governing diffusion of the heavy holes, neglecting excitation and diffusion of light holes.  The heavy holes form  a doublet with spin $\pm \frac{1}{2}$ and total angular momentum $j_z = \pm \frac{3}{2}$; we write the    charge density $N$ and the spin densities $S_i$ as a $4-$vector $\vec{\rho} = \begin{bmatrix} N, & S_x, & S_y, & S_z \end{bmatrix} $.     At time scales larger than the elastic scattering time $\tau$ their  diffusion and coupling to each other are controlled by the partial differential equation $\mathcal{D}_{ij}^{-1} \vec{\rho} = 0$,   where the  $4 \times 4$ matrix  $\mathcal{D}_{ij}$ is called the diffuson and is determined by  $ \mathcal{D}_{ij}(\vec{q}, \omega)  = (1 - I)^{-1}$. The joint scattering operator $I_{ij}$ is  is given by the integral 
\begin{equation} 
I_{ij} = \frac{\hbar}{4 \pi \nu \tau} \int d\vec{k} \, {Tr}(\, G^A( \vec{k} - \vec{q}/2, E_F)  \\ \sigma_i \, G^R(  \vec{k} + \vec{q}/2,  E_F + \hbar \omega) \, \sigma_j \, ) 
\label{JointScatteringOperator}
\end{equation}
$G^A$ and $G^R$  are the disorder-averaged single-particle Green's functions which express uncorrelated movements of $\psi$ and $\psi^\dagger$, while $\vec{q}$ is the diffuson momentum.  The trace is taken over the spin indices of $G^A, G^R, \sigma_i,$ and $\sigma_j$, which are all $2 \times 2$ matrices in spin space.  

We have computed the diffuson  systematically to next to leading order in $1/E_F \tau$ and to fourth order in $E_{SO}/E_F$.  This involved going to fourth order in $E_F \tau, q,$ and $ E_{SO}/E_F$.    However the dominant physics is already visible at second order.  At this order the diffuson simplifies to:
\begin{eqnarray}
 \mathcal{D}^{-1}  &  = & 1- \frac{1}{ 16 \pi  (1 - \imath \omega \tau)}\sum_{s,\acute{s},\nu= \pm 1}   \int d\theta  M \;  \frac{1}{1 - \imath \tau \nu \hbar^{-1} (E_F -  E(\vec{k} +  \nu \vec{q} , (1 - \nu)s/2 + (1 + \nu)\acute{s}/2) )}, 
 \nonumber \\
E_F - E & \approx & -2 \nu q \cos(\theta - \theta_q)  E_F/k_F 
+  \nu (s- \acute{s})\sqrt{a_x^2 + a_y^2}
   \nonumber \\
 M & = &  \begin{bmatrix}1 + s \acute{s}  & (s  + \acute{s}) \cos \phi &   (s + \acute{s} ) \sin \phi & 0  \\ (s  + \acute{s} ) \cos \phi & 1 + s \acute{s} \cos(2 \phi)  & s \acute{s} \sin(2 \phi) & -\imath (s - \acute{s}) \sin \phi \\  (s  + \acute{s} ) \sin \phi & s \acute{s} \sin(2 \phi)  &1-s \acute{s} \cos(2 \phi) & \imath (s - \acute{s}) \cos \phi  \\ 0 & \imath (s - \acute{s}) \sin \phi & - \imath ( s-  \acute{s}) \cos \phi & 1 - s \acute{s}  \end{bmatrix}
\nonumber \\
\cos \phi & = & a_x(\vec{k})/ \sqrt{a_x^2 + a_y^2}, \; \sin \phi  =  \pm \sqrt{1 - \cos^2 \phi} = a_y(\vec{k})/\sqrt{a_x^2 + a_y^2},  
\nonumber \\
\vec{k} & = & k_F(\cos \theta \hat{x} +  \sin \theta \hat{y}), \vec{q}  = q(\cos \theta_q \hat{x} +  \sin \theta_q \hat{y})
\end{eqnarray}

After some algebra we obtain: 
\begin{eqnarray}
 \mathcal{D}^{-1} \tau^{-1}  &  = &  
-\imath \omega  + D q^2  + \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & U  \cos 2\theta+ 1/T & U  \sin 2\theta & - \imath (C_1 q \cos(\theta_q + 2 \theta) + C_2 q_x)  \\ 0 & U  \sin 2\theta &  - U \cos 2 \theta + 1/ T &  -\imath (C_1 q \sin(\theta_q  + 2 \theta) - C_2   q_y)    \\ 0 &   \imath (C_1 q \cos(\theta_q + 2 \theta) + C_2 q_x)      &  \imath (C_1 q \sin(\theta_q  + 2 \theta) - C_2   q_y) &  2/ T \end{bmatrix}
\nonumber \\
1/T & = &   k_F^2 (\alpha  / \hbar)^2   \tau (k_F^4(1 + (\gamma_2/\gamma_3)^2) + 2 \beta^4)   \propto (E_{SO} \tau / \hbar)^2 / \tau
\nonumber \\
D & = & 2 (E_F \tau / \hbar k_F)^2 / \tau = v_F^2  \tau / 2 \propto (E_{SO} \tau / \hbar)^2 (E_F / E_{SO} k_F)^2  / \tau
\nonumber \\
U & = &  2 (\alpha  / \hbar)^2  \beta^2 k_F^4  \tau (1 - (\gamma_2/\gamma_3))   \propto (E_{SO} \tau / \hbar)^2 / \tau
\nonumber \\
 C_1 &=& 4 \alpha \beta^2   E_F  \tau /\hbar^2 \propto (E_{SO} \tau / \hbar)^2 (E_F / E_{SO} k_F)  / \tau, \;\;\;   
 \nonumber \\
 C_2 &=& (1 - (\gamma_2/\gamma_3))\, 2 \alpha k_F^2   E_F \tau /\hbar^2 \propto (E_{SO} \tau / \hbar)^2 (E_F / E_{SO} k_F) / \tau
\end{eqnarray}

Going to higher order, we found that all the non-zero elements of $ \mathcal{D}^{-1}$ have corrections.  There is also a spin-charge coupling at higher order, but it is small even compared to the corrections which we just mentioned.


%: $\imath e^{-\imath 2 \theta} = e^{\imath \xi}, \; - 2 \theta = \xi - \pi/2, \; \theta = \pi/4-\xi/2$


 
 Assuming the time-dependence $\exp{(-i \omega t -  i \vec{q}\cdot \vec{r})}$, we obtain the equations of motion:
\begin{eqnarray}
\partial_t N & = & D \nabla^2 N 
\nonumber \\ 
\partial_t S_x &=& D \nabla^2 S_x     - (U  \cos 2\theta + \frac{1}{T})S_x - U  \sin 2\theta \, S_y  -  C_2 \partial_x S_z - C_1 (\cos 2 \theta  \, \partial_x - \sin 2 \theta  \,  \partial_y )S_z
\nonumber \\ 
\partial_t S_y & = & D \nabla^2 S_y - U  \sin 2\theta \, S_x - (-U  \cos 2\theta +\frac{1}{T})S_y    +  C_2 \partial_y S_z - C_1 (\sin 2 \theta \,  \partial_x + \cos 2 \theta  \, \partial_y )S_z
\nonumber \\ 
  \partial_t S_z & =  &D \nabla^2 S_z -\frac{2}{T} S_z     +  C_2 (\partial_x S_x  - \partial_y S_y) +  C_1 (\cos 2 \theta  \, \partial_x - \sin 2 \theta  \, \partial_y )S_x + C_1 (\sin 2 \theta  \,  \partial_x + \cos 2 \theta  \,  \partial_y )S_y
\end{eqnarray}
 
When the strain is along the $x$ axis ($\theta = 0$) this simplifies to:
\begin{eqnarray}
\partial_t N & = & D \nabla^2 N 
\nonumber \\ 
\partial_t S_x &=& D \nabla^2 S_x     - (U   + \frac{1}{T})S_x   -  C_2 \partial_x S_z - C_1  \partial_x S_z
\nonumber \\ 
\partial_t S_y & = & D \nabla^2 S_y  - (-U   +\frac{1}{T})S_y    +  C_2 \partial_y S_z - C_1  \partial_y S_z
\nonumber \\ 
  \partial_t S_z & =  &D \nabla^2 S_z -\frac{2}{T} S_z     +  C_2 (\partial_x S_x  - \partial_y S_y) +  C_1  \partial_x S_x + C_1   \partial_y S_y
\end{eqnarray}

%\begin{appendix}
\bibliography{PSHHole}
\end{document}
