 Solution

There is a state-assumption issue. For a number-conserving Gaussian quantum state with no anomalous contractions, put ${\psi }_i=\psi ({\mathbf{r}}_i)$, $n_i=⟨{\psi }_i^{†}{\psi }_i⟩$ and the one-body density matrix $g_{ij}=⟨{\psi }_i^{†}{\psi }_j⟩$. Wick's theorem pairs each creation operator with an annihilation operator. The direct pairing gives $n_1{n}_2$; the exchanged pairing gives $ϵg_{12}{g}_{21}$, with $ϵ=+1$ for bosons and $-1$ for fermions. Hencenormal ordering removes the equal-position delta-function self-counting term. Fermionic occupation-number product states are Slater determinants and obey this Wick theorem. Independent thermal bosonic modes also obey it, with $N_{\alpha }$ understood as mean occupations.

A bosonic Fock product state does not generally obey the asserted identity. The bosonic number-state corrections to Gaussian density correlations follow from $⟨a_{\alpha }^{†2}{a}_{\alpha }^2⟩=N_{\alpha }(N_{\alpha }-1)$ and giveFor one occupied mode this is $N(N-1)|{\phi }_1{|}^2|{\phi }_2{|}^2$, whereas the uncorrected Gaussian expression gives $2N^2|{\phi }_1{|}^2|{\phi }_2{|}^2$. TODO: the bosonic state specification needs the Gaussian/thermal assumption, or the number-state correction must be retained. The Gaussian formulas below give the intended result under that stated qualification. The ensemble distinction is also explicit in a later primary treatment of ideal-gas correlations.^[ref]