Suppose that $\hat{\Psi}^\dagger(x)$, $\hat{\Psi}(x)$ are the usual field operators in second quantisation for some identical particle, and that $\hat{c}^\dagger_n$, $\hat{c}_n$ are the creation and annihilation operators in some discrete single-particle basis, with spatial wavefunctions $\phi_n(x)$ for each mode.
Suppose I know the one body density matrix in that discrete basis - that is, I know $\langle \hat{c}_n^\dagger \hat{c}_m\rangle$ for all $n$ and $m$. I want to use this to calculate the spatial density $\rho(x)=\langle\hat{\Psi}^\dagger(x)\hat{\Psi}(x)\rangle$. If I were to guess this by intuition, I would have rather confidently said that it would just be the sum of the spatial densities in each mode - that is: $$ \rho(x)=\sum_n |\phi_n(x)|^2\langle\hat{c}^\dagger_n\hat{c}_n\rangle $$ However, this seems to be incorrect. If I take the definition of $\rho(x)$ and expand out the field operators as $\hat{\psi}(x)=\sum_n \phi_n(x)\hat{c}_n$, $\hat{\psi}^\dagger(x)=\sum_n \phi_n^*(x)\hat{c}^\dagger_n$, then I find: $$\rho(x)=\langle\hat{\Psi}^\dagger(x)\hat{\Psi}(x)\rangle=\sum_{nm}\phi_n^*(x)\phi_m(x)\langle \hat{c}_n^\dagger \hat{c}_m\rangle$$ where we have off-diagonal terms for which $n\neq m$ contributing! This surprises me. Is there a good physical intuition for why these inter-mode terms contribute to the spatial density (or have I made a mistake, or do they cancel out somehow)?