Origin of potential term in non-relativistic QFT? The Hamiltonian of the free real scalar field is
$\hat H = \int d^3k (\hbar \omega_k) \hat a^\dagger_k \hat a_k$
where $\omega_k = c \sqrt{k^2 + \mu^2}$. In the limit where $|k| \ll \mu$,
$\omega_k \approx \mu c + \frac{ck^2}{2 \mu} = mc^2 + \frac{\hbar^2 k^2}{2m}$
where $m = \hbar \mu / c$ has the dimensions of mass. In this regime,
$\hat H \approx \int d^3 k( mc^2 + \frac{\hbar^2 k^2}{2m}) \hat a^\dagger_k \hat a_k = \int d^3{x} \Big( mc^2 \hat a(x)^\dagger \hat a(x) - \frac{\hbar^2}{2m} \hat a(x)^\dagger \nabla^2 \hat a(x) \Big).$
This, in a sense, explains how QM fits into QFT.
However, often when people give this example of a non-relatvistic quantum field theory (which is the result of taking the subset of the Hilbert space where large modes aren't excited) they also include a potential energy term to the Hamiltonian.
$$\hat V = \iint d^3{x} \hspace{0.1cm} d^3{y}  \hspace{0.1cm} \hat a(x)^\dagger \hat a(y)^\dagger V(|x - y|) \hat a(x) \hat a(y)$$
where we imagine that $V(r)$ is the Coulomb potential or something similar.
I understand in general how such a Coulomb potential can arise. Namely, it will arise when our field is coupled to another field. However, I can not come up with an argument as to why $\hat V$, in the form given, arises. (I assume it requires that particles are not moving at relativistic speeds for one.) Can anyone come up with an explanation for this explicit form of $\hat V$, similar to the one I gave for the origin of the free non-relativistic Hamiltonian?
 A: This answer concerns the motivation for the form a 2-body interaction term takes in the second-quantized formalism. Take the Hilbert space $L^2(\mathbb{R}^3)$ as the single particle Hilbert space and consider the bosonic Fock space over that. I didn't want to post this, because I got confused about a factor of 2 (please point out the mistake to me!), but since noone else seems to answer it.. The correct calculation should be included in most books on Condensed Matter physics, none of which I can call my own.
Given a 2-body potential as a function $V(|x-y|)$, defining an operator $\hat{V}$ on Fock space representing this potential amounts to finding an operator $\hat{V}$ which acts on two particle position eigenstates as
$$\hat{V}|x, y\rangle =V(|x-y|)|x,y\rangle.$$
That is to say, the state having exactly 2 particles, one of which is at $x$ and the other at $y$, hast potential energy $V(|x-y|)$.
By inserting an identity (for the 2-particle sector) $I^{(2)}=\int dx dy |x,y\rangle \langle x,y |$, one can see that this operator hast the form
$$\hat{V}^{(2)}= \int dx dy \ |x,y\rangle V(|x-y|)\langle x,y|$$
when restricted to the 2-particle subspace. This argument would be completely valid, if we had the fock space of distinguishable particles. However, we want bosons, so we insert symmetrization projectors and hope that this will fix it (for fermions one could insert anti-symmetrization operators here)
$$V^{(2)}=S \left( \int dx dy \ |x,y\rangle V(|x-y|)\langle x,y| \right) S$$
The zero-particle and one-particle restrictions are zero. As a generalization to more than 2 particles, we ask ourselves what the potential energy is when you have $n$ particles at positions $x_1, \dots,x_n$ subject to that potential. This is a classical question, which is answered by
$$V=\frac{1}{2}\sum_{i\neq j} V(|x_i - x_j|)$$
Thus we define
$$\hat{V}^{(n)}|x_1,\dots, x_n\rangle=\frac{1}{2}\sum_{i\neq j} V(|x_i - x_j|)|x_1,\dots, x_n\rangle$$
This leads by the same argument to
$$\hat{V}^{(2)}=\int dx_1\dots dx_n \ |x_1,\dots, x_n\rangle \frac{1}{2}\sum_{i\neq j} V(|x_i - x_j|)\langle x_1,\dots, x_n|$$
We have thus completely defined $\hat{V}$ on the Fock space. To see that this is the same as your expression, apply it to a position eigenstate:
$$\int dx dy\ \hat{a}(x)^{\dagger} \hat{a}(y)^{\dagger} V(|x-y|)\hat{a}(x)\hat{a}(y)|x_1,\dots x_n\rangle$$
To do this calculation first look at
$$\hat{a}(x)\hat{a}(y)|x_1,\dots x_n\rangle=\hat{a}(x)\sum_{i=1}^{n} \delta(y-x_i)|x_1,\dots (\text{no } x_i)\dots, x_n \rangle=$$
$$=\sum_{i\neq j}\delta(x-x_j) \delta(y-x_i)|x_1,\dots (\text{no } x_i, x_j)\dots, x_n \rangle$$
and now (being somewhat lazy with the symmetrization)
$$\int dx dy\ \hat{a}(x)^{\dagger} \hat{a}(y)^{\dagger} V(|x-y|)\hat{a}(x)\hat{a}(y)|x_1,\dots x_n\rangle= \int dx dy\ V(|x-y|)\hat{a}(x)^{\dagger} \hat{a}(y)^{\dagger} \sum_{i\neq j}\delta(x-x_j) \delta(y-x_i)|x_1,\dots (\text{no } x_i, x_j)\dots, x_n \rangle=S \int dx dy\ V(|x-y|) \sum_{i\neq j}\delta(x-x_j) \delta(y-x_i)|x,y,x_1,\dots (\text{no } x_i, x_j)\dots, x_n \rangle$$
with the presence of the deltas we can rename the positions in the state. Also, due to symmetrization $S$ it doesn't matter at which spot those stand anyway.
$$=S \int dx dy\ V(|x-y|) \sum_{i\neq j}\delta(x-x_j) \delta(y-x_i)|x_i,x_j,x_1,\dots (\text{no } x_i, x_j)\dots, x_n \rangle=S \int dx dy\ V(|x-y|) \sum_{i\neq j}\delta(x-x_j) \delta(y-x_i)|x_1\dots x_n \rangle$$
and the result is, as disired up to a factor of 2 (probably due a stupid mistake of mine)
$$S\sum_{i\neq j}V(|x_i-x_j|)|x_1\dots x_n \rangle$$
