How does QFT predict the probability density to find a particle at x? In quantum mechanics, the probability density of a particle's position is $$\rho(x)=|\langle x|\psi\rangle|^2$$
What is the corresponding expression in QFT to predict this distribution? Since $\rho(x)$ can be measured (at least to some accuracy) in experiment, this seems like a fair question to ask.
I have asked this question to a few people at my university. To anticipate a few answers which in my opinion did not resolve this question, let me point out that


*

*In QFT, because particle number is not conserved, it has been suggested that the question I am asking might not make sense. A slightly reformulated version of the question does, though: "what is the probability density of finding (for example) exactly one electron at x?".

*I have read some answers to similar questions like this one, which mentions a position operator defined by Newton & Wigner, so that its eigenstates give the probability of a position outcome as in QM. According to the same post, this attempt at an X operator failed because the probability was not Lorentz invariant.

*U(1) symmetry in $\mathscr{L}$ suggests something like a probability density, but this quantity can be negative, which is why it is interpreted as charge density instead (as far as I have read).

*I am aware that QFT is not usually used to predict position probability densities. However such things are commonly measured, for example in the double slit experiment. Therefore if the theory is at least as general as QM, it is desirable that it has something to say here.


Of course if one of these bullet points somehow does resolve the question I would be very interested to better understand how.
 A: The required amplitude is an overlap between a localized particle state and a field configuration state.  
For the case of a free scalar field, Roman Jackiw in his work on Functional representations for quantized fields has performed the computation of the exactly required amplitude (equation 2.14A ). 
I'll repeat here some details of his derivation and recast the result into a coherent state basis, which I think it is more amenable to experimental measurement.
In the Schrodinger picture, we associate a vector $|\phi\rangle$ in Hilbert space to each field configuration. Please notice that all the following expressions are infinite dimensional generalizations of the harmonic oscillator. (The following computations are not rigorous, the rigor can be improved by weighting and smearing, but this will not be done here).
The action of the field operator at time $t=0$ on the configuration vectors is given by:
$$\hat{\Phi}(x)|\phi\rangle = \phi(x)|\phi\rangle$$
($x$ is a vector in three dimensions )
The action of the canonical momenta (at $t=0$ )is given by:
$$\hat{\Pi}(x)|\phi\rangle = \frac{\delta}{\delta\phi(x)}|\phi\rangle$$
For a quadratic Hamiltonian:
$$\hat{H} = \int \mathrm{d}x \, \hat{\Pi}(x)^2 + \int \mathrm{d}x \, \mathrm{d}y\,  \hat{\Phi}(x)\Omega(x,y) \hat{\Phi}(y) $$
($\Omega$ can be viewed as an infinite dimensional mass matrix)
The vacuum functional has the form:
$$\Psi_{\Omega}(\phi) = \langle \phi | \Omega \rangle = \mathrm{det}^{\frac{1}{4}}\left({\frac{\Omega}{\pi}}\right) e^{-\frac{1}{2}  \int \mathrm{d}x \, \mathrm{d}y\,  \Phi(x)\Omega(x,y) \Phi(y)} $$
This is an infinite dimensional generalization of the Harmonic oscillator vacuum wave function.
In infinite dimensions, a wave functional with a different mass matrix will be orthogonal to the vacuum wave functional and to all of its excited state, thus will belong to another superselection sector. There are advantages in keeping a general mass matrix because in infinite dimensions they allow approximate treatment of weakly interacting fields through a Bogoliubov transformation. However, for simplicity, I'll continue with a diagonal mass matrix:
$$\Omega(x,y) \propto \delta(x-y)$$
(The corresponding vacuum and its wave functional will be denoted by: $|0\rangle$ and $\Psi_{0}(\phi) $ respectively). Also, I'll not bother with the infinite normalizations.
In analogy, the particle creation and annihilation operators are given by:
$$A(x) =  \frac{1}{\sqrt{2}}(\hat{\Phi} (x)+ i  \hat{\Pi} (x))$$
$$A^{\dagger}(x) =  \frac{1}{\sqrt{2}}(\hat{\Phi} (x) - i  \hat{\Pi} (x))$$
Thus, the required amplitude is given by:
$$\langle \phi | x\rangle = \langle \phi|A^{\dagger}(x)| 0\rangle> = \phi(x) \langle \phi|0\rangle =  \phi(x) e^{-\frac{1}{2}  \int \mathrm{d}x \,  \phi(x)^2} $$ 
The Schrodinger representation of a coherent state is given by analogy to the harmonic oscillator:
$$\Psi_{\alpha}(\phi) =  N e^{-\frac{1}{2}  \int \mathrm{d}x \,  (\Phi(x) - \sqrt{2}\alpha(x))^2 }$$
and the coherent state is given by:
$$|\alpha\rangle= \int \mathcal{D}\phi \Psi_{\alpha}(\phi) |\phi\rangle$$
(The integration is over all the field configurations). It is easy to see that $|\alpha\rangle$ is an eigenstate of the annihilation operator.
Thus we obtain:
$$\langle \alpha| x \rangle = \int \mathcal{D}\phi \Psi_{\alpha}(\phi) \langle \phi | x\rangle = \alpha(x)$$
In quantum systems described by coherent states of field theories , the function $\alpha (x)$ is usually referred to as the macroscopic wave function. The required amplitude has a very simple expression in terms of the macroscopic wave function.
Interpretation of the results
The requested amplitude $\langle \phi | x\rangle$ in the question is an overlap between a field configuration $\phi$ and a single particle state. The computed overlap is proportional to $\phi(x) e^{-\frac{1}{2}  \int \mathrm{d}x \,  \phi(x)^2}$. This result, also computed by Jackiw, makes sense, since it is proportional to the field strength at $x$. However, this amplitude will be in general be small unless the configuration $\phi$ is peaked at $x$. Thus, in practice this amplitude should be useful only when the system is in a lumped state around $x$, for example in a solitonic state. (If I had taken proper care of the normalizations, the space integral  the squared modulus of the result would have been automatically normalized to $1$).
The last computation was intended to find the amplitude $\langle \alpha | x\rangle$, in a coherent state rather in a Schrödinger state. The second last equation is just the change of basis between the coherent state basis and the Schrödinger basis. The complex field $\alpha$ is called the macroscopic wave function for the following reason: 
In a nonrelativistic system described by a complex Schrödinger field, the canonical momentum is equal to the complex conjugate of the field: $\hat{\Pi}(x)= \hat{\Phi}^{*}(x)$. Since we must have $[\hat{\Phi}(x), \hat{\Pi}(x')] = \delta(x-x')$, we must interpret the field as the annihilation operator and its momentum as the creation operator (not as in the relativistic case, where they are combinations of creation and annihilation operators). In this case, the expectation value of the field in the coherent state basis will be:
$$\langle \alpha| \hat{\Phi}(x) | \alpha \rangle = \langle \alpha| A(x) | \alpha \rangle  =  \alpha(x)$$
(The identity $A(x) |\alpha \rangle =\alpha(x)|\alpha \rangle $, valid in the coherent state basis, was used).
In relativistic systems such as radiation or relativistic plasmas we can still call the function $\alpha(x)$ a macroscopic wave function although the above reasoning is not valid. 
The reason that I gave you also the result in a coherent state basis is that whenever a distributed system is in a coherent state (such as a Bose-Einstein condensate); the macroscopic wave function is a measurable quantity. 
