Second order relativistic corrections to Pauli equation from Dirac equation

Question

I'm trying to derive the full and correct Hamiltonian for spin$\frac{1}{2}$ particles from Dirac equation up to second order in $v/c$. For a potential and magnetic field constant in time.

In particular, there should be Zeeman term, relativistic kinetic term, spin-orbit term and Darwin term:

$$\hat{H}=\frac{\mathbf{P}^2}{2m}+U-e \hbar \mathbf{ \sigma} \mathbf{B}-\frac{\mathbf{P}^4}{8m^3 c^2}- \frac{e \hbar}{4 m^2 c^2} \mathbf{ \sigma} (\mathbf{F} \times \mathbf{P})- \frac{e \hbar^2}{8 m^2 c^2} (\nabla \mathbf{F}) \tag{1}$$

Here:

$$\mathbf{P}=- i \hbar \nabla-e \mathbf{A} \\ \mathbf{B}= \nabla \times \mathbf{A} \\ \mathbf{F}=- \frac{1}{e} \nabla U$$

This is taken from various sources, for example, Berestetsky, Lifshitz, Pitaevsky "Quantum electrodynamics", (33.12) though I have added the Zeeman term and changed some notation.

Now, assuming the above is correct, I want to derive this expression from Dirac equation directly.

I use the following form of the equation as a starting point:

$$\begin{cases} \displaystyle \Psi_+ = - \frac{\sigma \mathbf{P}}{mc} \Psi_-+\frac{E}{mc^2} \Psi_+ \\ \displaystyle \Psi_- = \frac{\sigma \mathbf{P}}{mc} \Psi_+ -\frac{E}{mc^2} \Psi_- \end{cases} \tag{2}$$

Where $\Psi_{\pm}$ are two component spinors and $\sigma$ is the vector of Pauli matrices.

Introducing a new energy:

$$E=mc^2+\mathcal{E} \\ \mathcal{E}=i \hbar \frac{\partial}{\partial t} -U $$

We obtain:

$$\begin{cases} \displaystyle \frac{\mathcal{E}}{mc^2} \Psi_+= \frac{\sigma \mathbf{P}}{mc} \Psi_- \\ \displaystyle \left(2+ \frac{\mathcal{E}}{mc^2} \right)\Psi_- = \frac{\sigma \mathbf{P}}{mc} \Psi_+ \end{cases} \tag{3}$$

Assuming $$\Vert \mathcal{E} \Vert \ll mc^2$$

We obtain, up to first order in $\Vert \mathcal{E} \Vert / mc^2$:

$$\mathcal{E} \Psi_+ = \frac{\sigma \mathbf{P}}{2m} \left(1- \frac{\mathcal{E}}{2mc^2} \right) \sigma \mathbf{P} \Psi_+ \tag{4}$$

I assume that this should directly give us (1).

However, after some simplifications, while I get most of the terms above right (including the most important SO term), I get an incorrect Darwin term (twice as large) and also a bunch of extra terms.

My full derivation is lengthy, so I wanted to ask first if I missed something? Is (4) correct (up to second order in $v/c$) and should it lead to (1) without any extra terms?

I have extensively used the identity:

$$(\sigma \mathbf{a}) (\sigma \mathbf{b})= \mathbf{a} \mathbf{b}+i \sigma (\mathbf{a} \times \mathbf{b})$$

though it's not always clear how to apply it to operators. In any case, it gives a lot of extra terms not present in (1) when applied to (4).

---

There's a related question, but only about Pauli equation, which is simple enough to derive.

Answer 1

Let me provide some steps of the derivation here to see where the problem had been.

$$\mathcal{E} \Psi_+ = \frac{\sigma \mathbf{P}}{2m} \left(1- \frac{\mathcal{E}}{2mc^2} \right) \sigma \mathbf{P} \Psi_+$$

Let's change the notation a little:

$$\mathcal{E}=E_c-U$$

$$E_c \Psi_+ = \frac{(\sigma \mathbf{P})^2}{2m} \Psi_+ +U \Psi_+ -\frac{E_c}{4m^2 c^2} (\sigma \mathbf{P})^2 \Psi_+ + \frac{1}{4m^2 c^2} ( \sigma \mathbf{P} U \sigma \mathbf{P} )\Psi_+ $$

$$E_c \Psi_+ = \left(1- \frac{(\sigma \mathbf{P})^2}{4m^2 c^2}\right) \left( \frac{(\sigma \mathbf{P})^2}{2m} +U + \frac{1}{4m^2 c^2} ( \sigma \mathbf{P} U \sigma \mathbf{P} ) \right)\Psi_+ $$

$$E_c \Psi_+ = \left( \frac{(\sigma \mathbf{P})^2}{2m} +U - \frac{(\sigma \mathbf{P})^4}{8m^3 c^2}- \frac{(\sigma \mathbf{P})^2 U}{4m^2 c^2} + \frac{1}{4m^2 c^2} ( \sigma \mathbf{P} U \sigma \mathbf{P} ) \right)\Psi_+ $$

It's easy to identify several correct terms from (1) immediately.

One important note: Berestetsky et al derive the expression (1) sans the Zeeman term for the case when no vector potential is present. Which is probably why I initially got all the extra terms.

In this simple case we can write (defining $\mathbf{P}=\mathbf{p}-e \mathbf{A}$):

$$(\sigma \mathbf{P})^2=\mathbf{p}^2$$

Another important note: the probability density up to $v^2/c^2$ has the following form:

$$\rho=|\Psi_+|^2+|\Psi_-|^2=|\Psi_+|^2+ \frac{1}{4m^2 c^2} |\sigma \mathbf{p} \Psi_+|^2$$

To preserve the usual definition, according to Berestetsky et al, we have to definte another wavefunction:

$$\Psi_+= \left(1- \frac{\mathbf{p}^2}{8m^2 c^2}\right) \Psi$$

Which leaves us with:

$$E_c \Psi =\left(1+ \frac{\mathbf{p}^2}{8m^2 c^2}\right) \left( \frac{\mathbf{p}^2}{2m} +U - \frac{\mathbf{p}^4}{8m^3 c^2}- \frac{\mathbf{p}^2 U}{4m^2 c^2} + \frac{1}{4m^2 c^2} ( \sigma \mathbf{p} U \sigma \mathbf{p} ) \right) \left(1- \frac{\mathbf{p}^2}{8m^2 c^2}\right)\Psi $$

Up to second order we have:

$$E_c \Psi = \left( \frac{\mathbf{p}^2}{2m} +U - \frac{\mathbf{p}^4}{8m^3 c^2}- \frac{\mathbf{p}^2 U}{8m^2 c^2}- \frac{ U \mathbf{p}^2}{8m^2 c^2} + \frac{1}{4m^2 c^2} ( \sigma \mathbf{p} U \sigma \mathbf{p} ) \right) \Psi $$

Which coincides with the expression provided in Berestetsky et al before they derive (1) (without the Zeeman term I added).

So I guess my question is answered.

If someone adds a more in depth and/or intuitive explanation for all this stuff, I will gladly award them the bounty.