Second order relativistic corrections to Pauli equation from Dirac equation 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.
 A: 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. 

