Spin Orbital Coupling matrix in $p$-orbital basis So I have the following Hamiltonian inherited from atomic Physics:
$$H_\mathrm{SOC}=\alpha \vec{L}\cdot \vec{S}=\frac{\alpha}{2}(L^{+}\sigma^{+}+L^{-}\sigma^{-}+ L^{z}\sigma^{z})$$
Where $L$ is the angular momentum, $S$ is spin, and $L^{\pm} (\sigma^{\pm})$ is the angular momentum (spin) step up/step down operator.
Now in the basis of $p$ orbitals and spin: $\{|p_x{\uparrow}\rangle, |p_x{\downarrow}\rangle, |p_y{\uparrow}\rangle, |p_y{\downarrow}\rangle, |p_z{\uparrow}\rangle, |p_z{\downarrow}\rangle\}$, we get the following $6\times 6$ matrix:
$$H_\mathrm{SOC}=\frac{\alpha}{2}\begin{pmatrix}
 0&-i  &0  &0  &0  &1 \\ 
 i& 0 &0  &0  &0  &1 \\ 
 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{pmatrix}$$
So excuse my ignorance, but how exactly are these matrix elements calculated?
I understand the first matrix element is the energy of the spin orbital coupling of the $\langle p_x{\uparrow}|$ electron acting on the $|p_x{\uparrow}\rangle$ electron, since they are they same orientation the energy, sure, should be zero. But now we have the  $\langle p_y{\uparrow}|$ electron acting on the $|p_x{\uparrow}\rangle$ orbital and we get $i$. How is this calculated? Can someone show the steps to calculate one matrix element so I can see how this is done. And no this is not homework or anything, just my personal curiosity.
I found something similar online, but they used Clebsch–Gordan coefficients which confused me more.
 A: First, I think your $H_{SOC}$ is incorrect and should read:
$$
H_{SOC}=\frac{\alpha}{2}\left(L_+\sigma_- + L_-\sigma_+ + L_z\sigma_z\right)\, .
\tag{1}
$$
The twist with your problem is that it uses the real spherical harmonics, rather than the complex exponential form more common in physics and better adapted to the evaluation of matrix elements.  The conversion between the two can be found on this wiki page.
Thus:
\begin{align}
p_y&=\frac{i}{\sqrt{2}}\left(Y_{1}^{-1}+Y_1^1\right)\quad\Rightarrow \quad \vert p_y\rangle = \frac{i}{\sqrt{2}}\left(\vert 1,-1\rangle + \vert 1,1\rangle\right)\, ,\tag{2}\\
p_x&=\frac{1}{\sqrt{2}}\left(Y_{1}^{-1}-Y_1^1\right)
\quad\Rightarrow \quad \vert p_y\rangle = \frac{i}{\sqrt{2}}\left(\vert 1,-1\rangle + \vert 1,1\rangle\right)\, ,\tag{3}\\
p_z&=Y_1^0\qquad\qquad \qquad\quad\Rightarrow \quad \vert p_z\rangle = \vert 1,0\rangle\, , \tag{4}
\end{align}
where the arguments $(\theta,\varphi)$ of the $Y_{\ell}^{m}$ have been omitted for clarity.
The next step is to recall the action of the various operators on $Y_{\ell m}$:
\begin{align}
\hat L_\pm Y_{\ell}^{m}&=\sqrt{(\ell\mp m)(\ell\pm m+1)}Y_{\ell}^{m\pm 1}\, ,\tag{5}\\
\hat L_z Y_{\ell}^{m}&= mY_{\ell}^{m}\tag{6}
\end{align}
from which one obtains for instance (if I did not make obvious mistakes), the matrix representation of $\hat L_+$ in the basis $\{\vert p_x\rangle, \vert p_y\rangle, \vert p_z\rangle\}$
\begin{align}
\hat L_+\vert p_x\rangle &= \vert p_z\rangle \qquad 
\hat L_+\vert p_y\rangle = i \vert p_z\rangle \qquad 
\hat L_+\vert p_z\rangle =-\vert p_x\rangle - i\vert p_y\rangle\,  , \tag{7}
\end{align}
so that 
\begin{align}
\hat L_+\to \left(\begin{array}{ccc}
0&0&-1\\
0&0&-i\\
1&i&0\end{array}\right)\, ,\qquad \hat L_-\to 
\left(\begin{array}{ccc}
0&0&1\\
0&0&-i\\
-1&i&0\end{array}\right)\, . \tag{8}
\end{align}
since the matrix representation for $\hat L_-$ will be transpose conjugate of the matrix representation for $\hat L_+$.  Likewise:
\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\, . \tag{9}
\end{align}
Hence:
\begin{align}
\hat L_+\hat \sigma_- \vert p_x;\uparrow\rangle &= \Bigl[\hat L_+\vert p_x\rangle \Bigr]
\Bigl[\hat \sigma_-\vert \uparrow\rangle\Bigr]= \vert p_z\rangle \vert\downarrow\rangle= \vert p_z;\downarrow\rangle\, , \tag{10}\\
\hat L_-\hat \sigma_+ \vert p_x;\uparrow\rangle &= 0 \, ,\tag{11} \\
\hat L_z\hat\sigma_z\vert p_x;\uparrow\rangle &=
 \Bigl[\hat L_z\vert p_x\rangle \Bigr]
\Bigl[\hat \sigma_z\vert \uparrow\rangle\Bigr]=i\vert p_y\rangle \vert \uparrow\rangle=
i\vert p_y;\uparrow\rangle \tag{12}
\end{align}
This is not quite the first column of your matrix, but as I understand it this matrix is not quite correct either since it ought to be hermitian but is not: for instance the matrix element in position $(2,6)$ should be the complex conjugate of the matrix element in position $(6,2)$. 
Anyways, my calculation would give the elements in the first column as
$$
(0,0,i,0,0,1)^T \tag{13}
$$
This might be because the basis elements are ordered as 
$\{\vert p_x\uparrow\rangle,\vert p_y\uparrow\rangle,\vert p_z\uparrow\rangle,\vert p_x\downarrow\rangle,\vert p_y\downarrow\rangle,\vert p_z\downarrow\rangle\}$ rather than the ordering you give (or else it's quite possible I made an error somewhere).  With this ordering my first column would come out as $(0,i,0,0,0,1)^T$, which is identical to yours.  Note that this does not fix the hermiticity problem mentioned earlier.
The other columns are found in the same manner.
A: For example look at entry 1,2:
$(H_{SOC})_{1,2} \,=\,  <p_x\uparrow| H_{SOC}  |p_x\downarrow>\,=\, <p_x\uparrow| \frac{\alpha}{2}(L^{+}\sigma^{+}+L^{-}\sigma^{-}+ L^{z}\sigma^{z})|p_x\downarrow> \,=\,$
$=\,\frac{\alpha}{2}(<p_x\uparrow|L^{+}\sigma^{+}|p_x\downarrow> + <p_x\uparrow|L^{-}\sigma^{-}|p_x\downarrow> + <p_x\uparrow|L^{z}\sigma^{z}|p_x\downarrow>).$
Thats all getting lenghty, so lets look
only at the first term, here the spin-ladder operator only acts on the spin part of the product and ladder up makes spin-up from spin-down so we get:
$<p_x\uparrow|L^{+}\sigma^{+}|p_x\downarrow>\,=\,<\uparrow \uparrow><p_x|L^{+}|p_x>$,
in the expression to the right of the equal sign the "integral" over the functions of only the spin coordinate is on the left and the integral over space on the right, and with the above outlined stuff from AccidentialFourierTransform (its actually about how to express the spherical harmonics in a basis which can be evaluated using the ladder operators) that gives
$1\cdot<p_x \frac{1}{\sqrt{2}}|0 + |0> \,=\,\frac{1}{\sqrt{2}}<p_x|0>\,= ... \,=\, \frac{1}{\sqrt{2}}\cdot 0=0.$
This is about what you have to do to evaluate these matrix elements (up to some mistakes).
A: If we start on a basis of $\{|p_{-1}{\downarrow}\rangle, |p_{-1}{\uparrow}\rangle, |p_{0}{\downarrow}\rangle, |p_{0}{\uparrow}\rangle, |p_{1}{\downarrow}\rangle, |p_{1}{\uparrow}\rangle\}$ we get for the spin-orbit coupling interaction
$H_{lm\sigma}=\left(
\begin{array}{cccccc}
 \frac{1}{2} & 0 & 0 & 0 & 0 & 0 \\
 0 & -\frac{1}{2} & \frac{1}{\sqrt{2}} & 0 & 0 & 0 \\
 0 & \frac{1}{\sqrt{2}} & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & \frac{1}{\sqrt{2}} & 0 \\
 0 & 0 & 0 & \frac{1}{\sqrt{2}} & -\frac{1}{2} & 0 \\
 0 & 0 & 0 & 0 & 0 & \frac{1}{2} \\
\end{array}
\right)$
We can transform the basis from the complex spherical harmonics to the real spherical harmonics with the matrix $T = \left(
\begin{array}{cccccc}
 \frac{1}{\sqrt{2}} & 0 & 0 & 0 & -\frac{1}{\sqrt{2}} & 0 \\
 0 & \frac{1}{\sqrt{2}} & 0 & 0 & 0 & -\frac{1}{\sqrt{2}} \\
 \frac{i}{\sqrt{2}} & 0 & 0 & 0 & \frac{i}{\sqrt{2}} & 0 \\
 0 & \frac{i}{\sqrt{2}} & 0 & 0 & 0 & \frac{i}{\sqrt{2}} \\
 0 & 0 & 1 & 0 & 0 & 0 \\
 0 & 0 & 0 & 1 & 0 & 0 \\
\end{array}
\right)$
The transformation $T^*.H.T^T$ then yields the Hamiltonian on the basis $\{|p_{x}{\downarrow}\rangle, |p_{x}{\uparrow}\rangle, |p_{y}{\downarrow}\rangle, |p_{y}{\uparrow}\rangle, |p_{z}{\downarrow}\rangle, |p_{z}{\uparrow}\rangle\}$: $H_{\gamma,\sigma}=\left(
\begin{array}{cccccc}
 0 & 0 & \frac{i}{2} & 0 & 0 & -\frac{1}{2} \\
 0 & 0 & 0 & -\frac{i}{2} & \frac{1}{2} & 0 \\
 -\frac{i}{2} & 0 & 0 & 0 & 0 & -\frac{i}{2} \\
 0 & \frac{i}{2} & 0 & 0 & -\frac{i}{2} & 0 \\
 0 & \frac{1}{2} & 0 & \frac{i}{2} & 0 & 0 \\
 -\frac{1}{2} & 0 & \frac{i}{2} & 0 & 0 & 0 \\
\end{array}
\right)$
Note that there is no coupling between the $p_x$ orbitals with opposite spin. The spin orbit coupling interaction can only raise or lower the angular momentum by 1, the $p_x$ orbital is a linear combination f $m=\pm1$ and thus not coupled by the spin-orbit interaction. (i.e. there are several typos in your original matrix)
I calculated these matrices and transformations using the package Quanty http://www.quanty.org
