Could you say exactly what the density matrix $\rho(x,y)$ describes when $x \neq y$ in QM? I understand that if $\psi$ is a many-body electronic state function and
$$
g(x) = \int \psi^*(x,x_2,\cdots,x_n) \psi(x,x_2,\cdots,x_n) \, dx_2\cdots dx_n,
$$
then $g$ describes the probability to measure a particle in a volume element at $x$. And if I add the prefactor $n$, it becomes the density, simply because the density is the probability of measuring one electron in a volume element scaled by the number of electrons in the system (as the electrons are identical). But if instead I focus on this function
$$
f(x,y) = \int \psi^*(x,x_2,\cdots,x_n) \psi(y,x_2,\cdots,x_n) \, dx_2\cdots dx_n,
$$
what exactly is the physical interpretation of $f$? Alternatively, what is the physical interpretation of $n\cdot f$ (if it's easier to understand)? Please explain in the simplest possible terms, I am trying to wrap my head around it, but struggling.
 A: This is known as the one-body correlation function, sometimes noted $g^{(1)}(x,y)$. In some way it describes the link between the wavefunction at two different points in space, when tracing over all other particles (in terms of second quantization, it describes the correlation between the field $\psi(x)$ and its conjugate $\psi^{\dagger}(y)$).
For $y=x$ it is equal to the local density up to some prefactor (depending on your normalization of the wavefunction), but as $y$ differs more and more from $x$ it can be thought as a description of how the field at position $x$ affects the field at position $y$.
For systems with long-range order, such as a Bose-Einstein condensate, $g^{(1)}(x,y)$ stays finite when $|x-y| \to \infty$, otherwise it goes to $0$.
Since it can be hard to visualize the meaning of this "field", you can also consider the two-body correlation function, $$g^{(2)}(x,y) = \int \psi^*(x,y,x_3,...,x_N) \psi(x,y,x_3,...,x_N) dx_3 ... dx_N$$
It can be shown that $g^{(2)}$ is linked to the density-density correlations of the system, that is how the presence of a particle at position $x$ affects the probability to find a particle at position $y$. Furthermore, when considering non-interacting Bosonic (resp. Fermionic) systems, $g^{(2)}$ can be reduced in a combination of $g^{(1)}$ terms.
$$\langle n(x) n(y) \rangle \sim g^{(2)}(x,y) = \langle n(x) \rangle \langle n(y) \rangle + \left| g^{(1)}(x,y) \right|^2 \mathrm{(resp. } -\left| g^{(1)}(x,y) \right|^2)$$
Hope that helps.
