Solution (source code)

= 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=\langle\psi_i^\dagger\psi_i\rangle$ and the <one-body density matrix> $g_{ij}=\langle\psi_i^\dagger\psi_j\rangle$. <Wick's theorem> pairs each creation operator with an annihilation operator. The direct pairing gives $n_1n_2$; the exchanged pairing gives $\epsilon g_{12}g_{21}$, with $\epsilon=+1$ for <bosons> and $-1$ for <fermions>. Hence
$$
\boxed{C_{12}=n_1n_2+\epsilon|g_{12}|^2}.
$$
<normal 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.

\b[A bosonic Fock product state does not generally obey the asserted identity.] The <bosonic number-state corrections to Gaussian density correlations> follow from $\langle a_\alpha^{\dagger2}a_\alpha^2\rangle=N_\alpha(N_\alpha-1)$ and give
$$
\boxed{C_{12}^{\rm Fock}=n_1n_2+|g_{12}|^2
-\sum_\alpha N_\alpha(N_\alpha+1)|\varphi_\alpha(\mathbf r_1)|^2
|\varphi_\alpha(\mathbf r_2)|^2}.
$$
For one occupied mode this is $N(N-1)|\varphi_1|^2|\varphi_2|^2$, whereas the uncorrected Gaussian expression gives $2N^2|\varphi_1|^2|\varphi_2|^2$. \b[This counterexample disproves the asserted identity for arbitrary bosonic number-product states.] The formulas below give both the Gaussian result, with its additional state assumption, and the exact bosonic number-state alternative. The ensemble distinction is also explicit in a later primary treatment of ideal-gas correlations.<Quantum Correlations of Ideal Bose and Fermi Gases in the Canonical Ensemble by Kazumasa Tsutsui and Takafumi Kita>[Abstract: canonical reduced density matrices versus grand-canonical Wick factorization; 2016]{ref}