Searching for a shell
Spin-orbital coupling in tight binding fomalism
Intrinsic SOC does not depend on the symmetry of the system and can be expressed by:
\[H_\mathrm{SO}=\lambda \hat{L} \cdot \hat{S}\]Where:
\[\begin{align} \hat{L} & = [\hat{L_x}, \hat{L_y}, \hat{L_z}] = -i \hbar (\vec{r} \times \nabla) \\ \hat{S} & = \frac{1}{2}[\hat{\sigma_x}, \hat{\sigma_y}, \hat{\sigma_z}] \end{align}\]To make things simple, we can re-write this with non-hermition raising and lowering operators $\boldsymbol{S_{\pm}}$ and $\boldsymbol{L_{\pm}}$:
\[H_\mathrm{SO}=\frac{1}{2} \lambda (\hat{L}_{+}\hat{\sigma}_{-} + \hat{L}_{-}\hat{\sigma}_{+} + \hat{L}_{z}\hat{\sigma}_{z})\]Where:
\[\begin{align} \hat L_+ &= \hat L_x + i \hat L_y \\ \hat L_- &= \hat L_x - i \hat L_y \end{align}\]and the raising and lowering $\hat L_{\pm}$ act on complex spherical harmonics, resulting in:
\[\begin{align} \hat L_\pm Y_{\ell}^{m}&=\sqrt{(\ell\mp m)(\ell\pm m+1)}Y_{\ell}^{m\pm 1}\\ \hat L_z Y_{\ell}^{m}&= mY_{\ell}^{m}, \end{align}\]while the raising and lowering $\hat \sigma_{\pm}$ act on spin components of the wavefunction, resulting in:
\[\begin{align} \hat \sigma_+\vert \uparrow \rangle &=0\, , \quad \hat \sigma_+\vert \downarrow \rangle =\vert \uparrow\rangle\, ,\quad \hat \sigma_-\vert \uparrow \rangle =\vert\downarrow\rangle\, , \quad \hat \sigma_-\vert \downarrow \rangle =0 . \end{align}\]Note thtat the usual real spherical harmonics like $p_x$, $p_y$ $p_z$ are made up by linear combinate spherical harmonics:
\[\begin{array}{l} p_{x}=\frac{1}{\sqrt{2}}\left(Y_{1}^{-1}-Y_{1}^{1}\right) \\ p_{y}=\frac{i}{\sqrt{2}}\left(Y_{1}^{-1}+Y_{1}^{1}\right) \\ p_{z}=Y_{1}^{0} \end{array}\]For example, if we are only considering the basis of:
\[\left\{\left|p_{x}, \uparrow\right\rangle,\left|p_{y}, \uparrow\right\rangle,\left|p_{z}, \uparrow\right\rangle,\left|p_{x}, \downarrow\right\rangle,\left|p_{y}, \downarrow\right\rangle,\left|p_{z}, \downarrow\right\rangle\right\}\]then we can calculate the matrix elements of $\hat{\boldsymbol{L}}_{\pm}$ by:
\[\begin{align} \hat{L}_{+} \vert p_x\rangle &= \hat{L}_{+} \frac{1}{\sqrt{2}}\left(Y_{1}^{-1}-Y_{1}^{1}\right)\\ &= \frac{1}{\sqrt{2}} (\sqrt{(1+1)(1-1+1)}Y_{1}^{0}-0) \\ &= Y_1^0 = \vert p_z \rangle \\ \hat{L}_{+} \vert p_y\rangle &=i \vert p_z \rangle \\ \hat{L}_{+} \vert p_z\rangle &= -\vert p_x \rangle - i \vert p_y \rangle \end{align}\]Hence, the $\hat{\boldsymbol{L}}_{\pm}$ operators can be written as:
\[\begin{array}{ccc} L_{+} = \left(\begin{array}{ccc} 0 & 0 & -1 \\ 0 & 0 & -i \\ 1 & i & 0 \end{array}\right) , & L_{-} = \left(\begin{array}{ccc} 0 & 0 & 1 \\ 0 & 0 & -i \\ -1 & i & 0 \end{array}\right), & L_{z} = \left(\begin{array}{ccc} 0 & -i & 0 \\ i & 0 & 0 \\ 0 & 0 & 0 \end{array}\right) \\ \sigma_{+} = \left(\begin{array}{cc} 0 & 1 \\ 0 & 0 \end{array}\right) , & \sigma_{-} = \left(\begin{array}{cc} 0 & 0 \\ 1 & 0 \end{array}\right), & \sigma_{z} = \left(\begin{array}{cc} 1 & 0 \\ 0 & -1 \end{array}\right). \end{array}\]And the SOC hamiltonian can be expressed in the matrix form:
\[H_{\mathrm{SO}}=\frac{\lambda}{2}\left(\begin{array}{cccccc} 0 & -i & 0 & 0 & 0 & 1 \\ i & 0 & 0 & 0 & 0 & -i \\ 0 & 0 & 0 & -1 & i & 0 \\ 0 & 0 & -1 & 0 & i & 0 \\ 0 & 0 & -i & -i & 0 & 0 \\ 1 & i & 0 & 0 & 0 & 0 \end{array}\right)\]This is what’s I’ve implemented in the 
wanSOC package. For example, to generate the above mentioned Hamiltonian, we can:
from __future__ import print_function
import numpy as np
from wanSOC.basis import *
from wanSOC.io import *
from wanSOC.helper import *
from wanSOC.hamiltonian import *
orb = 'p'
q_axis = [0,0,1]
H_SO = gen_Hsoc(orb, q_axis)
print(H_SO)