Ultrarelativistic limit for neutrinos: Why is this approximation working? In this section of the Wikipedia article on neutrino oscillations, a neutrino mass eigenstate $\left|\nu_i\right>$ is written as
$$\left|\nu_i(t)\right> = e^{-i(E_it-\vec p_i\cdot\vec{x})} \left|\nu_i(0)\right>$$
where $E_i$ is the energy of the mass eigenstate $i$ (rest energy + kinetic energy, as I understand). Then the first approximation is made, using $|\vec p_i| = p_i \gg m_i$ to obtain
$$E_i = \sqrt{p_i^2 + m_i^2} \approx p_i + \frac{m_i^2}{2p_i}.$$
I'm fine with that step since that is just a Taylor expansion. But the next approximation confuses me:
$$p_i + \frac{m_i^2}{2p_i} \approx E + \frac{m_i^2}{2E}$$
where $E$ is "the total energy of the particle". How is this energy $E$ different than the energy $E_i$ from above? Even if $E$ is the "total energy", this sounds like we would expect that $E_i$ is smaller than the "total" energy $E$. But our approximation gives $E_i \approx E + \frac{m_i^2}{2E}$, so $E_i$ should be larger than $E$.
 A: You have a point that there is a bit of doublethink in the WP article and many texts confusing subtlety mavens as reviewed, and resolved, in  Lipkin 2006, published.  (NB This paper, however, prefers the equivalent fixed energy option, unlike the equivalent old-fashioned fixed momentum option used here.)
Much of the confusion arises from unstated assumptions and "dancing about QM" thinking, instead of pursuing full QM calculations which lose the reader. I would just state the PDG 2019 summary/mnemonic p 9, (14.38). I'm calling it "mnemonic", as it is a memorable picture of several simplistic summaries, starting from different assumptions, yielding the same answer (and meant to discourage you from different summaries yielding wrong mnemonics).
The conceit is that your wave packet consisting of, say, two mass eigenstates, with $m_1<m_2$ are produced in entanglement at the meson decay "source" at a  sloppy point, so with common momentum p. Since their masses differ, slightly, so do their energies and hence their velocities,
$$
E_i=\sqrt{p^2+m_i^2}\approx p+ \frac {m_i^2}{2p}.
$$
So the wavepacket simulacrum ~ $e^{-itE_1} +ae^{-itE_2}$  evolves  like
$$
e^{-it(p+m_1^2/2p)} \left (1+ ae^{-it(m_2^2-m_1^2)/2p}\right ),
$$
so with "total" energy $E\approx p$ unobservable in the pure collective phase of the wavepacket, and differing from the momentum by $O(m^2)$, so we might as well put it in the denominator of the interesting  oscillation phase following a in lieu of the momentum, since  this only contributes slop of $O(m^4)$!
Again, $t\approx L $ for the wave packet, so the interesting oscillating probability at the detection point  amounts to
$$
\left |1+ae^{-iL(m_2^2-m_1^2)/2E}\right |^2=1+a^2 + 2a \sin \left ( {L(m_1^2-m_2^2)\over 2E}\right ) ~.
$$
