Why is the Wick contraction in HFB or BCS equal to a single-particle density?

I'm trying to understand how in Hartree-Fock-Bogoliubov (HFB) or BCS theory we can write a product of creation/annihilation operators as single-particle densities under the guise of "Wick's theorem".

In this set of slides, on slide no. 20 (and I've seen this in many many, many papers) the author makes the claim $c_a^\dagger c_b = \rho_{ba} + :c_a^\dagger c_b:$, invoking Wick's theorem, where $\rho_{ba} \equiv \left.\langle\Phi\right| c_a^\dagger c_b\left|\Phi\rangle\right.$. Essentially, the Wick contraction of $c_a^\dagger c_b$ is equal to $\rho_{ba}$. I'd like to put this on firm footing. My two questions are:

1. How do we know the contraction $\langle c_a^\dagger c_b \rangle = \rho_{ba}$?
2. Given #1, the normal ordered product $:c_a^\dagger c_b:$ must be 0 in the state $\left|\Phi\rangle\right.$ by definition. Are we free to choose "normal-ordering" with respect to any state we like?

That is, can we say $c_a^\dagger c_b$ is normal-ordered with respect to the vacuum $\left|-\rangle\right.$, so $\langle c_a^\dagger c_b \rangle$ should be 0 there but isn't normal-ordered with respect to the state $\left|\Phi\rangle\right.$ where $c^\dagger_a c_b \ne :c^\dagger_a c_b:$ (whatever $:c^\dagger_a c_b:$ actually is)? Is there a way to prove that some combination $c_a^\dagger c_b$, $c_a c_b$, etc. always has a vanishing normal-ordered product (I think I can show this with Boguliubov quasiparticles in this case)?

If we have the freedom of #2, then #1 is trivial without having to calculate anything, since $$\rho_{ba} = \left.\langle\Phi\right| c_a^\dagger c_b\left|\Phi\rangle\right. = \left.\langle\Phi\right| \langle c_a^\dagger c_b \rangle + :c_a^\dagger c_b:\left|\Phi\rangle\right. = \langle c_a^\dagger c_b \rangle.$$

From there we can do the more general decomposition $$\left.\langle\Phi\right| c_a^\dagger c_b^\dagger c_d c_c \left|\Phi\rangle\right. = \rho_{ca} \rho_{db} - \rho_{da}\rho_{cb} + \kappa_{ba}^* \kappa_{cd},$$

where $\kappa_{ba} = \left.\langle\Phi\right| c_a c_b \left|\Phi\rangle\right.$, using Wick's theorem.

Normal ordering is an ordering of a product of field operators, in which all the annihilation operators are placed in to the right of all creation operators. This mean that the expectation value of a normal ordering in relation to the vacuum state is zero.

If $c_a^\dagger c_b = \rho_{ba} + :c_a^\dagger c_b:$ and $|\Phi\rangle$ is the vacuum state, then $\rho_{ba} = \left.\langle\Phi\right| c_a^\dagger c_b\left|\Phi\rangle\right.$ because $\langle \Phi|:c_a^\dagger c_b:|\Phi\rangle=0$ .

Is important to note that the definition of normal odering is dependent of the definition of vacuum state. Each set of creation and annihilation operator are defined in relation to some vacuum state, which imply that the normal ordering is implicitly related to vacuum state definition. In Fermi sea our definition of vacuum state is the state of fermi energy filled by electrons for bellow, and for this we have a suitable definition of Normal ordering compatible to the vacuum state (The particle-hole picture).

Often the consequences of Wicks theorem are thought of in the context of expectation values, however the notation quickly becomes unclear when you do this, because e.g. contractions depend implicitly on the state used for the expectation values.

At its roots, Wicks theorem should be used for operators alone, then you take the expectation value with respect to whatever state you choose. Granted, the resulting Wick expansion is only useful if you consider the expectation value with respect to the vacuum since all normal ordered terms go to zero (hence why this step is usually skipped).

Wick's theorem can let us re-write a string of operators in terms of normal ordered pieces and contractions. We use Wick's theorem for a single pair of operators to define the contraction as the difference between the time ordered set and the normal ordered set. For time independent operators, the given order is considered time ordered. e.g:

$\langle c^\dagger_a c_b \rangle \equiv \ T(c^\dagger_a c_b) - :c^\dagger_a c_b:\ =\ c^\dagger_a c_b - :c^\dagger_a c_b:$

The problem with giving the contraction a value at this point is that it depends on the basis of operators we are considering (and the corresponding vacuum). If we consider non-relativistic fermion operators, the best we can do is say that it is a c-number. With a little foresight we denote this c-number $\rho_{ba}$ (where the flipped indices are by convention).

$\langle c^\dagger_a c_b \rangle \equiv \ \rho_{ba} \qquad \text{therefore} \qquad c^\dagger_a c_b = \rho_{ba}\ + :c^\dagger_a c_b:$

Now we take a few examples of specific bases.

Single Particle Basis (bare vacuum $\lvert 0 \rangle$):

1. If the operators $c_a$ are single particle operators $a_a, a^\dagger_a$ in the single particle basis, then we can take the equation defining the contraction at face value.

$\langle a^\dagger_a a_b \rangle = a^\dagger_a a_b - :a^\dagger_a a_b:\ = 0$.

Therefore the vacuum density is $\rho_{ba} = \langle 0 \rvert a^\dagger_a a_b \lvert 0 \rangle = 0$

2. If the operators $c_a$ are quasiparticle operators $b_a, b^\dagger_a$ then we need the bogoliubov transformation coefficients to express them in the single particle basis.

$\langle b^\dagger_a b_b \rangle \equiv b^\dagger_a b_b - :b^\dagger_a b_b:\ = (V^T_{ai} c_i + U^T_{ai} c^\dagger_i) (U^\dagger_{bj} c_j + V^\dagger_{bj} c^\dagger_j)\ -\ :(V^T_{ai} c_i + U^T_{ai} c^\dagger_i) (U^\dagger_{bj} c_j + V^\dagger_{bj} c^\dagger_j)\:$

Where Einstein summations are in effect. Since all the differences are zero except the anti-normal ordered terms, we get

$\langle b^\dagger_a b_b \rangle = V_{ai}^T V^\dagger_{bj}(c_i c^\dagger_j - :c_i c^\dagger_j: ) = V_{ai}^T V^\dagger_{bj} \{c_i,c_j^\dagger \} = V_{ai}^T V^\dagger_{bj} \delta_{ij} = V^\dagger_{bi} V_{ia}$

Therefore the vacuum density is $\rho_{ba} = \langle 0 \rvert b^\dagger_a b_b \lvert 0 \rangle = V^\dagger_{bi} V_{ia}$

Quasi-Particle Basis (HFB vacuum $\lvert \phi \rangle = \prod_i b_i \lvert 0 \rangle$:

1. If the operators $c_a$ are single particle operators $a_a, a^\dagger_a$ in the quasi-particle basis, then we need the inverse bogoliubov transformation.

$\langle a^\dagger_a a_b \rangle = (V_{ai} b_i + U^*_{ai} b^\dagger_i) (U_{bj} b_j + V^*_{bj} b^\dagger_j) - :(V_{ai} b_i + U^*_{ai} b^\dagger_i) (U_{bj} b_j + V^*_{bj} b^\dagger_j):\ = V^*_{bi}V^T_{ia}$

Therefore the vacuum density is, $\rho_{ba} = \langle \phi \rvert a^\dagger_a a_b \lvert \phi \rangle = V^*_{bi}V^T_{ia}$

2. If the operators $c_a$ are quasiparticle operators in the quasiparticle basis, then again the calculation is simple.

$\langle b^\dagger_a b_b \rangle = b^\dagger_a b_b - :b^\dagger_a b_b:\ = 0$.

Therefore the vacuum density is, $\rho_{ba} = \langle \phi \rvert b^\dagger_a b_b \lvert \phi \rangle = 0$

You can play a similar game with the other 3 types of contractions, so you know how to express all the terms in your Wick expansion in a given basis. With this we can evaluate expectation values of any product state as long as we can express that state in terms of creation and annihilation operators and their corresponding vacuum. If you just write $\rho$ and $\kappa$, your Wick expansion will be valid in any basis.