I am reading through the derivation of the Gross-Pitaevskii equation in the Heisenberg picture, and I am having some trouble interpreting the following identificiation. In the derivation the field operator $\psi(r)$ is defined as the position space equivalent of the creation/annihilation operator $\hat{a}\_k$. $$\hat{\psi}(r,t) = \sum\_k \phi\_k(r, t) \hat{a}\_k.$$ Here $\phi\_k(r,t)$ is a single particle state of momentum mode $k$. The GP equation is now defined in terms of the following c-number. $$\psi\_0(r,t) = \langle \hat{\psi}(r,t)\rangle.$$ The derivation I understand, but what confused me is that this complex function $\psi\_0(r,t)$ is identified as the "wavefunction of the condensate". In other words the absolute value squared of this function gives the local particle density, and can thus be used to examine the spatial profile of the BEC.

I cannot seem to find a satisfactory explanation for this identification. Surely the density is actually related to the expectation value of the number operator, e.g. $\langle\psi^{\dagger}\psi\rangle$? Also I don't understand how this is even a wavefunction at all, what does this function have to do with probability of position?

Any help would be greatly appreciated!