Consider a scalar field $\phi$ described by the Klein-Gordon Lagrangian density
$L = \frac{1}{2}\partial_\mu \phi^\ast\partial^\mu \phi – \frac{1}{2} m^2 \phi^\ast\phi$.
As written in every graduate QM textbook, the corresponding conserved 4-current $j^\mu = \phi^\ast i \overset{\leftrightarrow}{\partial^\mu} \phi$ gives non-positive-definite $\rho=j^0$. If we are to interpret $\phi$ as a wave function of a relativistic particle, this is a big problem because we would want to interpret $\rho$ as a probability density to find the particle.
The standard argument to save KG equation is that KG equation describes both particle and its antiparticle: $j^\mu$ is actually the charge current rather than the particle current, and negative value of $\rho$ just expresses the presence of antiparticle.
However, it seems that this negative probability density problem appears in QFT as well. After quantization, we get a (free) quantum field theory describing charged spin 0 particles. We normalize one particle states $\left|k\right>=a_k^\dagger\left|0\right>$ relativistically:
$$ \langle k\left|p\right>=(2\pi)^3 2E_k \delta^3(\vec{p}-\vec{k}), E_k=\sqrt{m^2+\vec{k}^2} $$
Antiparticle states $\left|\bar{k}\right>=b_k^\dagger \left|0\right>$ are similarly normalized.
Consider a localized wave packet of one particle $\left| \psi \right>=\int{\frac{d^3 k}{(2\pi)^3 2E_k} f(k) \left| k \right>}$, which is assumed to be normalized. The associated wave function is given by
$$ \psi(x) = \langle 0|\phi(x)\left|\psi\right> = \int{\frac{d^3 k}{(2\pi)^3 2E_k} f(k) e^{-ik\cdot x}} $$
$$ 1 = \langle\psi\left|\psi\right> = \int{\frac{d^3 k}{(2\pi)^3 2E_k} |f(k)|^2 } = \int{d^3x \psi^\ast (x) i \overset{\leftrightarrow}{\partial^0} \psi (x)}$$.
I want to get the probability distribution over space. The two possible choices are:
1) $\rho(x) = |\psi(x)|^2$ : this does not have desired Lorentz-covariant properties and is not compatible with the normalization condition above either.
2) $\rho(x) = \psi^\ast (x) i \overset{\leftrightarrow}{\partial^0} \psi(x)$ : In non-relativistic limit, This reduces to 1) apart from the normalization factor. However, in general, this might be negative at some point x, even if we have only a particle from the outset, excluding antiparticles.
How should I interpret this result? Is it related to the fact that we cannot localize a particle with the length scale smaller than Compton wavelength ~ $1/m$ ? (Even so, I believe that, to reduce QFT into QM in some suitable limit, there should be something that reduces to the probability distribution over space when we average it over the length $1/m$ … )
Best Answer
The following does not completely answer OP's question, rather, this is going to be a clarification of subtleties and difficulties on the issue. Notice that sometimes I am going to use the same notations as OP used, but not necessarily the exact same meanings and I will make them clear in the context.
The usage of "wavefunction"
It has two possible meanings when people refer to something as wavefunctions :
(1). The collection of some functions $g(x)$ furnishes a positive-energy, unitary representation of the underlying symmetry group(in our case just the Poincare group).
(2). In addition to (1) being satisfied, we should be able to interpret $|g(x)|^2$ as the probability distribution of finding the particle in $(\mathbf{x},\mathbf{x}+d\mathbf{x})$. It really has to be of the form $|g(x)|^2$ according to the standard axioms of quantum mechanics, provided you interpret $g(x)$ as an inner product $\langle x|g\rangle$ where $\langle x|$ is an eigen bra of the position operator. I will discuss what position operator means later.
The 2nd meaning is of course much stronger and OP is searching for a wavefunction in this sense. However, the $\psi(x)$ written by OP is only a wavefunction in the 1st sense, because clearly $\langle 0|\phi(x)$ cannot be eigen bras of any Hermitian position operator, easily seen from the fact that they are not even mutually orthogonal, i.e. $\langle 0|\phi(x)\phi^\dagger(y)|0\rangle\neq0$ even when $x$ and $y$ are spacelike separated. As a consequence, $\int d^3\mathbf{x}|\psi(x)|^2\neq1$ as OP has already noted. Moreover, $|\psi(x)|^2$ is invariant under Lorentz transformation, but a density distribution should transform like the 0th component of a 4-vector in relativistic space-time, as OP has also noticed.
The localized states and position operator(Newton-Wigner)
In this section I will mostly rephrase(for conceptual clarity in sacrifice of technical clarity) what is written in this paper.
What is a sensible definition of position operator of single-particle states? First we need to think about what the most spatially-localized states are, and then it would be natural to call these states $|\mathbf{x} \rangle$, then it is also natural to call the operator having these states as the eigenstates the position operator. It seems pretty reasonable and not too much to ask to require localized states to have the following properties:
(a). The superposition of two localized states localized at the same position in space should again be a state localized at the same position.
(b). Localized states transform correctly under spatial rotation, that is, $|\mathbf{x} \rangle \to |R\mathbf{x} \rangle$ under a rotation $R$.
(c). Any spatial translation on a localized state generates another localized state that is orthogonal to the original, that is, $\langle\mathbf{x}+\mathbf{y}|\mathbf{x} \rangle=0$ if $\mathbf{y}\neq 0$.
(d). Some technical regularity condition.
It turns out these conditions are restrictive enough to uniquely define localized states $|\mathbf{x}\rangle$. It can be worked out that, borrowing OP's notations and normalization convention, if $|\psi\rangle=\int\frac{d^3 k}{(2\pi)^3 2E_k} f(k) |k\rangle$, then(including time dependence) $$\langle x|\psi\rangle=\int\frac{d^3 k}{(2\pi)^3 \sqrt{2E_k}} f(k) e^{-ik\cdot x},$$
and this (unsurprisingly) gives $\int d^3\mathbf{x}|\langle x|\psi\rangle|^2=1$. However, this is not the full solution to OP's problem, because we can show that, although not as bad as transforming invariantly, $|\langle x|\psi\rangle|^2$ is not as good as transforming like a 0th component, either. The underlying reason is, as already realized by Newton and Wigner, that a boost on a localized state will generate a delocalized state, so the interpretation is really frame dependent.
As I disclaimed, I do not know if there is a complete satisfactory solution, or if it is even possible, but I hope it helps to clarify the issue.
Appendix: Some interesting properties of Newton-Wigner(NW) states and operator
I decide to make it an appendix since I think this is not directly relevant yet very interesting(all the following are for scalar field, and NW also discusses spinor field in their paper):
(1). A state localized at the origin, projected to the bras $\langle 0|\phi(x)$, has the form
$$\langle 0|\phi(x)|\mathbf{x}=0\rangle=\left(\frac{m}{r}\right)^{\frac{5}{4}}H_{\frac{5}{4}}^{(1)}(imr),$$ where $r=(x_1^2+x_2^2+x_3^2)^{\frac{1}{2}}$ and $H_{5/4}^{(1)}$ is a Hankel function of the first kind. So it is not a delta function under such bras.
(2)The NW position operator $q_i (i=1,2,3)$ acting on momentum space wavefunction(defined as $f(k)$ in OP's notation) is $$q_if(k)=-i\left(\frac{\partial}{\partial k_i}+\frac{k_i}{2E_k}\right)f(k),$$ and in nonrelativistic limit the second term approaches 0, giving the familiar expression of a position operator. We can also get, if $\Psi(x)=\langle 0|\phi(x)|\Psi\rangle$, then $$q_i\Psi(x)=x_i\Psi(x)+\frac{1}{8\pi}\int\frac{\exp(-m|\mathbf{x-y}|)}{|\mathbf{x-y}|}\frac{\partial \Psi(y)}{\partial y_i}d^3\mathbf{y},$$ and here the nonrelativistic limit is hidden in the unit of $m$, when converting back to SI units, the $m$ on the exponent is really the inverse of Compton wavelength, which can be taken as $\infty$ for low energy physics, so again the 2nd term vanishes.
(3)$[q_i,p_j]=i\delta_{ij}$.