How to write the single electron spin-orbit coupling under an external magnetic field? As we know, without the external magnetic field,  the single electron spin-orbit coupling(SOC) has the form $\boldsymbol{\sigma}\cdot(\boldsymbol{\nabla} V\times \mathbf{p})$ up to a coefficient, where $V$ is the static electric potential.
Now if we add an static external magnetic field $\mathbf{B}=\boldsymbol{\nabla}\times \mathbf{A}$, does the SOC preserve the above form $\boldsymbol{\sigma}\cdot(\boldsymbol{\nabla} V\times \mathbf{p})$ or change into $\boldsymbol{\sigma}\cdot(\boldsymbol{\nabla} V\times (\mathbf{p}-e\mathbf{A}/c))$? The notation $\mathbf{p}$ always represents $-i\hbar(\partial_x,\partial_y,\partial_z)$ in all of our formulas.
This reminds me another question, may be worthless: Does this expression $\mathbf{r}\times(\mathbf{p}-e\mathbf{A}/c)$ make sense? How to write the orbital angular momentum under an external magnetic field?
Thank you very much.
 A: Yes, the substitution $\mathbf{p}\rightarrow\mathbf{p}-e\mathbf{A}/c$ you are suggesting should be mathematically valid. You’ll need to pick specific systems and look at the appropriate parameters to know if this substitution is physically relevant or not.  I say that because this is a higher-order effect. You may ask: higher order with respect to what? It is the next higher order with respect to the standard Zeeman term. Using the properties of the Pauli matrices ($\boldsymbol{\sigma}$) it is easy to see that the following equality is always true $$\frac{\mathbf{p}^{2}}{2m}=\frac{\left(\boldsymbol{\sigma}\cdot\mathbf{p}\right)^{2}}{2m}\, .$$ In the presence of an external gauge field $\mathbf{A}$ we get
\begin{align}
\frac{1}{2m}\left[\boldsymbol{\sigma}\cdot\left(\mathbf{p}-\frac{e\mathbf{A}}{c}\right)\right]^{2}&=\frac{1}{2m}\left(\mathbf{p}-\frac{e\mathbf{A}}{c}\right)^{2}+\frac{i}{2m}\boldsymbol{\sigma}\cdot\left[\left(\mathbf{p}-\frac{e\mathbf{A}}{c}\right)\times\left(\mathbf{p}-\frac{e\mathbf{A}}{c}\right)\right]\\&=\frac{1}{2m}\left(\mathbf{p}-\frac{e\mathbf{A}}{c}\right)^{2}-\frac{e\hbar}{2m}\boldsymbol{\sigma}\cdot\mathbf{B}
\end{align}
The second term $\propto \boldsymbol{\sigma}\cdot\mathbf{B}$ gives rise to the Zeeman effect; this is, in fact, the leading-order effect. This comes out naturally from the Dirac equation in the nonrelativistic limit. I borrowed the above expression from equation (3.19) of 

J. J. Sakurai. “Advanced quantum mechanics.” (1967).

I will make a couple of more references to Sakurai in this post. If you go through the derivations in section 3.3 carefully you’ll see how you can start with the Dirac equation, perform some series expansions in orders of $(v/c)^2$, and finally obtain a new Schrödinger equation, given by (3.83) in Sakurai (shown below), which is first order in $(v/c)^2$: $$\left[\frac{\mathbf{p}^{2}}{2m}+eA_{0}-\frac{\mathbf{p}^{4}}{8m^{3}c^{2}}-\frac{e\hbar}{4m^{2}c^{2}}\boldsymbol{\sigma}\cdot(\mathbf{E}\times\mathbf{p})-\frac{e\hbar^{2}}{8m^{2}c^{2}}\boldsymbol{\nabla}\cdot\mathbf{E}\right]\Psi=E^{({\rm NR})}\Psi\; .$$ A couple of comments about the above expression:


*

*The Zeeman term is missing because, in the derivation in the subsection “Approximate Hamiltonian for an electrostatic problem,” Sakurai assumes $\mathbf{A}=0$ to avoid messy expressions.

*The spin-orbit coupling term $\propto\boldsymbol{\sigma}\cdot(\mathbf{E}\times\mathbf{p})$ has appeared.
If you follow all the steps of the derivation carefully you’ll note that we can replace $\mathbf{p}\rightarrow\mathbf{p}-e\mathbf{A}/c$ everywhere in the above expression (including the spin-orbit coupling term). You will also need to add back the Zeeman term since now $\mathbf{A} \ne 0$. However, as I mentioned before, being a first order (in $(v/c)^2$) effect, the new $\boldsymbol{\sigma}\cdot(\mathbf{E}\times(\mathbf{p}-e\mathbf{A}/c))$ term will be very weak compared to the zeroth order Zeeman term.
Another thing to note: in the derivation, in equation (3.82), where Sakurai evaluates the commutator $$\left[\boldsymbol{\sigma}\cdot\mathbf{p},-ie\hbar\boldsymbol{\sigma}\cdot\mathbf{E}\right]=-e\hbar^{2}\boldsymbol{\nabla}\cdot\mathbf{E}-2e\hbar\boldsymbol{\sigma}\cdot(\mathbf{E}\times\mathbf{p})$$ he sets $\boldsymbol{\nabla}\times\mathbf{E}=0$ which is consistent with his earlier assumption of $\mathbf{A}=0$. Since I could not find an equivalent expression elsewhere (for $\mathbf{A} \ne 0$), I worked out the following: $$\left[\boldsymbol{\sigma}\cdot\mathbf{p},-ie\hbar\boldsymbol{\sigma}\cdot\mathbf{E}\right]=-e\hbar^{2}\boldsymbol{\nabla}\cdot\mathbf{E}-ie\hbar^{2}\boldsymbol{\sigma}\cdot(\boldsymbol{\nabla}\times\mathbf{E})-2e\hbar\boldsymbol{\sigma}\cdot(\mathbf{E}\times\mathbf{p})$$ Using Maxwell’s equations we can easily replace that with $$ \boldsymbol{\nabla}\times\mathbf{E}=-\frac{\partial\mathbf{B}}{\partial t}$$ and accommodate non-static magnetic fields as well. Note: this new term is kind of like the time derivative of the (zeroth order) Zeeman term.
The bottom line is that all effects related to the behavior of an electron, under any conditions, eventually come from the Dirac equation. Although not very convenient, it’s often safest to start at the source (i.e. the Dirac equation which is exact) and make approximations on a case-by-case basis. As for your last question: yes, one can meaningfully define an angular momentum in an electromagnetic field. In addition to Zeeman effect, you also might want to look at the Paschen-Back effect.
