Can the Pauli equation be used to derive the usual $\mathbf L\cdot \mathbf S$ spin-orbit coupling? If not, what's the simplest way? Can the Pauli equation,
$$
 \left[{\frac {1}{2m}}({\boldsymbol {\sigma }}\cdot (\mathbf {p} -q\mathbf {A} ))^{2}+q\phi \right]|\psi \rangle =i\hbar {\frac {\partial }{\partial t}}|\psi \rangle,
$$
be used as a starting point to derive the spin-orbit coupling hamiltonian, of the form $\mathbf L\cdot \mathbf S$, for the hydrogen atom? If so, how is this done? If not, what's the core reason for this, and what's the cleanest starting ground for deriving this coupling directly from a more fundamental Schrödinger equation?
I'm interested in sharp, clean arguments, as well as references that contain them.
 A: As the other two replies mentioned, the answer is no, spin-orbit coupling cannot be derived from the Pauli equation. Why? In a sense you could just say the term has been neglected in an approximation already made to arrive at the Pauli equation. The starting point you need is the Dirac equation. 
To my eyes, the most convenient way to find the spin-orbit term is to perform the Foldy-Wouthuysen (FW) transformation on the Dirac equation including a central Coulomb potential. My favorite reference for this, albeit not a very commonly used one, is Chapter 3.4-3.6 of Brown & Carrington, Rotational Spectroscopy of Diatomic Molecules. I'll just roughly reproduce their derivation here, but I've heard that a pdf of that book is not hard to find... 
Start with the Dirac Hamiltonian (in old-fashioned notation)
$${\mathcal H}= \beta m c^2 - e \phi + c \vec{\alpha} \cdot \vec{\pi}$$
$$ \equiv \beta m c^2 + V + O, $$
where $\vec{\pi}=\vec{P} + e \vec{A}$ and the operators $V=-e \phi$ and $O=\vec{\alpha}\cdot \vec{\pi}$ are even and odd with respect to the "big" and "small" parts of the electron wavefunction. The general idea of the Foldy-Wouthuysen transformation is to choose a unitary operator $e^{i {\mathcal S}}$ such that when we apply a unitary transformation
$$ {\mathcal H}' = e^{i {\mathcal S}} {\mathcal H} e^{-i {\mathcal S}} + i \hbar (\partial_t \, e^{i {\mathcal S}} ) e^{-i {\mathcal S}},$$
the odd parts of the transformed Hamiltonian are suppressed by powers of $c$. It turns out that this is possible, and the right choice of ${\mathcal S}$ is
$$ {\mathcal S} = - \frac{i \beta O}{2 m c^2}. $$
It so happens that you actually have to apply the FW transformation twice in order to suppress $O$ by two powers of $c$. After doing this, one of the resulting terms in the twice-transformed Hamiltonian is
$$ \frac{1}{8 m^2 c^2} \left[ \vec{\alpha} \cdot \vec{\pi} , \left[ \vec{\alpha} \cdot \vec{\pi}, (-e \phi) \right] + i \hbar \partial_t (\vec{\alpha} \cdot \vec{\pi}) \right] $$
$$ = -\frac{e \hbar^2}{8 m^2 c^2} \nabla \cdot \vec{\mathcal E} - \frac{g_S \mu_B}{2 m c^2} \vec{S} \cdot (\vec{\mathcal E} \times \vec{\pi}), $$
where $\vec{\mathcal E}=-\nabla \phi$ is the Coulomb electric field. The electric field $\vec{\mathcal E}$ can be replaced by the explicit Coulomb expression $\sim Z e^2 \vec{r}/r^3$, and the numerator $\vec{r} \times \vec{p}$ is precisely the orbital angular momentum $\vec{L}$. Thus we have a term in the Hamiltonian proportional to $\vec{S} \cdot \vec{L}$, which is spin-orbit coupling. Oh and by the way, the first term above I believe is the Darwin term proportional to $\delta^3(\vec{r})$.
The Foldy-Wouthuysen transformation can also be used to derive the Pauli term $\sim \vec{\sigma} \cdot \vec{\mathcal B}$ and the familiar leading order relativistic correction to the kinetic energy, proportional to $p^4$. For a free particle, it can also be shown to account for the phenomenon of zitterbewegung.
A: Spin orbit coupling is a relativistic effect, so cannot be derived from the non-relativistic Schroedinger equation plus the magnetic moment term. This  is what the Pauli equation is once one does some manipulation.  The ${\bf L}\cdot {\bf S }$ term can  be added by hand however. 
A: The best way to derive spin orbit and Pauli interaction from the Dirac equation in my opinion is to follow the approach in Itzykson & Zuber. They solve the Dirac equation for the hydrogen atom by passing via the squared Dirac equation.
