Recovering the position-space probability density from QFT Lately I've been trying to wrap my head around the relationship between quantum fields and the wave function of non-relativistic quantum mechanics. It is well-known that QFT, through its demotion of position to a simple label, lacks a position operator $\hat{x}$ and hence a position basis $|x\rangle$. A direct consequence of this is that you cannot construct a position-space wavefunction $\psi(x)=\langle x|\psi\rangle$. Despite this, experiments such as the electron double slit experiment show that there is indeed some kind of position-space probability density you can reconstruct, at least in the non-relativistic limit. A natural question then arises: what form does the probability density approximated by $|\psi(x)|^2$ take in QFT?
This excellent answer by Chiral Anomaly cleared up many of my doubts regarding a similar question in the context of the photon double slit experiment. The classical EM wave interference pattern is realized by a superposition of coherent states corresponding to the two slits. This pattern acts as a probability distribution for localized detections of single photons, strikingly similar to that expected by a naive Born rule interpretation. From how I understood it, the localization of the photons works because, while the creation/annihilation operators are strictly non-local, their commutators with the field are rapidly-decreasing functions of space-like separation: $\left[F^\pm(\vec x,t),\,F(\vec y,t)\right] \sim |\vec x-\vec y|^{-4}$. So for coarse enough resolution in the detector they can be considered local.
This same kind of localization concept applies for electrons, only now the commutator is an exponentially decreasing function of space-like separation $\left[\psi^\pm(\vec x,t),\,\psi(\vec y,t)\right] \sim e^{-m|\vec x-\vec y|}$. However, the interference pattern for electrons arises in a totally different way: by diffraction and interference of the electron wave function. The wave-like nature of light is explained through coherent states, but how is the position wave function of the electron explained by QFT?
Following Chiral Anomaly's work, I'd like to present here my initial attempt at trying to make sense of this mathematically. Please let me know whether this makes any sense or if I've badly misunderstood something :)
For simplicity, consider a massive scalar field $\phi$. We introduce the following particle number operator corresponding to each pixel, $i$, of a detector: $$N(\beta_i,t) := \int d^3x\ \beta_i(\vec x)\phi^+(\vec x,t)\phi^-(\vec x,t)\tag*{(1)}$$ where $\phi^\pm(\vec x,t)$ are the creation and annihilation parts of the field expansion and $\beta_i(\vec{x})$ is a real-valued smearing function that is zero outside the detector pixel $i$. This operator is non-local (so accurately counts particle number) but has a characteristic scale $1/m$ ($\hbar/mc$ in non-natural units) from the exponential decay of the commutator mentioned above. So as long as the pixel is larger than this length scale the operator can be considered local to a good approximation.
Consider an arbitrary state given by $$|\psi\rangle=\int\frac{d^3p}{(2\pi)^3}\frac{1}{2\omega_\vec{p}}\tilde{\psi}(\vec{p})|\vec{p}\rangle,\tag*{(2)}$$ where $\tilde{\psi}(\vec{p})$ is some square-integrable function which we will call the "momentum-space wave function". We then have that
\begin{align*}
\phi^-(\vec{x},t)|\psi\rangle&=\int\frac{d^3q}{(2\pi)^3}\frac{1}{\sqrt{2\omega_\vec{q}}}a_\vec{q} e^{-iq_\mu x^\mu}\int\frac{d^3p}{(2\pi)^3}\frac{1}{2\omega_\vec{p}}\tilde{\psi}(\vec{p})|\vec{p}\rangle\\
&=\int\frac{d^3q}{(2\pi)^3}\frac{d^3p}{(2\pi)^3}\frac{1}{2\omega_\vec{p}\sqrt{2\omega_\vec{q}}}e^{-iq_\mu x^\mu}\tilde{\psi}(\vec{p})a_\vec{q}|\vec{p}\rangle\tag*{(3)}
\end{align*}
using $a_\vec{q}|\vec{p}\rangle=(2\pi)^3\sqrt{2\omega_\vec{p}}\delta^{(3)}(\vec{p}-\vec{q})|0\rangle$, this reduces to
\begin{align*}
\phi^-(\vec{x},t)|\psi\rangle&=\int\frac{d^3p}{(2\pi)^3}\frac{1}{2\omega_\vec{p}}e^{-ip_\mu x^\mu}\tilde{\psi}(\vec{p})|0\rangle\tag*{(4)}
\end{align*}
In the non-relativistic limit we have that $|\vec{p}|\ll m$, so $\omega_\vec{p}\approx m$ and we arrive at
\begin{align*}
\phi^-(\vec{x},t)|\psi\rangle&\approx\frac{e^{-imt}}{2m}\left(\int\frac{d^3p}{(2\pi)^3}e^{i\vec{p}\cdot\vec{x}}\tilde{\psi}(\vec{p})\right)|0\rangle\tag*{(5)}
\end{align*}
We recognize the integral in the parentheses as the Fourier transform of the momentum-space wave function, which in QM is precisely the position space wave function $\psi(\vec{x})$. Inserting this into the expectation value for the particle number, we obtain $$\langle\psi|N(\beta_i,t)|\psi\rangle=\int d^3x\ \beta_i(\vec{x})\,\langle\psi|\phi^+(\vec x,t)\phi^-(\vec x,t)|\psi\rangle=\frac{1}{4m^2}\int d^3x\ \beta_i(\vec x)\,\psi^*(\vec{x})\psi(\vec{x}),\tag*{(6)}$$ which, for $\beta_i(\vec{x})=1$ inside the pixel, is simply the integrated probability density of the wave function across the pixel as in ordinary QM.
So this seems to work out...almost: the time-dependence has disappeared! So I'm not sure I'm going about this in the correct way (I suspect my equating of $\omega_p$ with $m$ may not be valid). My question is: does this work, and if not, what is the right way to recover the position space probability density in QFT?

Edit: Including the $p^2$ term in non-relativistic limit.
For $|\vec{p}|/m\ll 1$ we have, up to second order in $p$,
\begin{align*}
&\omega_p\approx m+\frac{\vec{p}^2}{2m}\tag*{(7)}\\\frac{1}{\omega_p}&\approx \frac{1}{m^2}\left(m-\frac{\vec{p}^2}{2m}\right)\tag*{(8)}
\end{align*}
So instead of equation 5 we have
\begin{align*}
\phi^-(\vec{x},t)|\psi\rangle&\approx\frac{e^{-imt}}{2m^2}\left(m\int\frac{d^3p}{(2\pi)^3}e^{-i\frac{p^2}{2m}t}e^{i\vec{p}\cdot\vec{x}}\tilde{\psi}(\vec{p})-\int\frac{d^3p}{(2\pi)^3}\frac{p^2}{2m}e^{-i\frac{p^2}{2m}t}e^{i\vec{p}\cdot\vec{x}}\tilde{\psi}(\vec{p})\right)|0\rangle\\
&=\frac{e^{-imt}}{2m^2}\left(m\,\psi(\vec{x},t)+\frac{1}{2m}\nabla^2\psi(\vec{x},t)\right)|0\rangle\tag*{(9)}
\end{align*}
where $\psi(\vec{x},t)$ is the (now) time-dependent position wave function and satisfies the free non-relativistic Schrödinger equation.
Inserting this into the integral for the number expectation value then gives, to second order in $|\vec{p}|,$
\begin{align*}
\langle\psi|N(\beta_i,t)|\psi\rangle&=\frac{1}{4m^2}\int d^3x\ \beta_i(\vec x)\left(\psi^*(\vec{x})\psi(\vec{x})+\frac{1}{2m^2}(\psi\nabla^2\psi^*+\psi^*\nabla^2\psi)\right)\tag*{(10)}
\end{align*}
Which is the ordinary probability density plus an extra "current" term $\sim (\psi\nabla^2\psi^*+\psi^*\nabla^2\psi)$$\sim (\psi\partial_t\psi^*-\psi^*\partial_t\psi)$ that doesn't seem to vanish when integrated as far as I can tell.
 A: When expanding the relationship $\omega=\sqrt{\vec p^2 + m^2}$ in powers of $\vec p^2/m^2$, we need to keep the $\vec p^2$ term:
$$
 \omega=m+\frac{\vec p^2}{2m} + O(\vec p^4).
\tag{1}
$$
The time-dependence survives in the $\vec p^2$ term because that term can't be factored out of the integral over $p$.
In any approximation, we should keep at least the lowest-order term that retains the feature we care about. If we care about a particle's ability to propagate through space (so that the expectation value of $N(\beta,t)$ varies with time), then we should at least keep the $\vec p^2$ term. To see this another way, recall that the non-relativistic Hamiltonian for a single particle is derived using (1). The $\vec p^2/2m$ term in (1) leads to the familiar term $-\nabla^2/2m$ in the Hamiltonian when single-particle states are represented by wavefunctions $\psi(\vec x)\sim\int d^3p\ e^{i\vec p\cdot\vec x}\psi(\vec p)$. The kinetic term $\nabla^2/2m$ is why the particle can propagate in the non-relativistic single-particle model, and that's why we need to retain the $\vec p^2$ term to get a time-dependent expectation value of $N(\beta,t)$.
