Long range order in the BCS ground state

I am trying to prove the following equivalent form for long-range-order in superconductivity (Annett's book states something like this) :

\begin{equation}\lim_{|R-R'|\to \infty}\langle\psi^{\dagger}(R')\psi(R)\rangle =\langle\psi^{\dagger}(R)\psi(R)\rangle\end{equation}

where the brackets are evaluated on the BCS ground state : $|\psi_{BCS}\rangle= ∏_{k}(u_k^{*}+v_k^* P_k^{\dagger})|0\rangle$ ($|0\rangle$ being the zero-temperature Fermi sea), $P_k^{\dagger}=c_{k\uparrow}^{\dagger}c_{-k\downarrow}^{\dagger}$ when :

$\psi^{\dagger}(R) = \int d^3r~\psi(r) \phi_{\uparrow}^{\dagger} (R+r/2)\phi_{\downarrow}^{\dagger}(R-r/2)$

($\psi(r)$ being some wavefunction) and :

$\phi_{\sigma}(r) = \sum_{k}e^{ikr} c_{k\sigma}$ ($c$ being the commun fermion annihilator, and the $P$s commuting like bosons for this particular state).(Sorry for the lengthy definitions)

The gap function being defined as: $\Delta \propto \sum_k <P_k>$, we can see relatively easily that the right side of the first equation will be proportional to $|\Delta|^2$. But what I've found quite puzzling is what to do with the inhomogeneous terms (especially in the left term) where the $c$ and $c^{\dagger}$ do not form a pairing operator and the Fermi sea not being ideal to deal with them...

Could somebody help me with this calculation? Or am I taking the problem on the wrong side?