Multi-channel mean field theory I have always been confused about the theoretical foundation of the mean field approximation. Below I follow the book Many-body Quantum Theory in Condensed Matter Physics by Bruus and Flensberg, Chapter 4.
For an interaction term $AB$ in a Hamiltonian, the mean field decoupling is done by ignoring higher order fluctuations of $A, B$ around some mean field average:
$$
\begin{align*}
    AB &= [\underbrace{
        (A - \langle A \rangle)
    }_{\text{fluctuation}} + \langle A \rangle]
    [\underbrace{
        (B - \langle B \rangle)
    }_{\text{fluctuation}} + \langle B \rangle]
    \\
    &\approx (A - \langle A \rangle) \langle B \rangle
    + \langle A \rangle (B - \langle B \rangle)
    + \langle A \rangle \langle B \rangle
    \\
    &= A \langle B \rangle + \langle A \rangle B
    - \langle A \rangle \langle B \rangle
    \tag{4.7}
\end{align*}
$$
However, in practice, we are usually decoupling 4 operators; thus we have multiple ways (called channels by physicists) to group them into two pairs, and perform the above decoupling. For example, consider the Hartree-Fock approximation for a system of fermions
$$
\begin{align*}
    H &= H_0 + V_{\text{int}}
    \tag{4.18a}
    \\
    H_0 &= \sum_\nu \xi_\nu c^\dagger_\nu c_\nu
    \tag{4.18b}
    \\
    V_{\text{int}} &= \frac{1}{2} 
    \sum_{\nu,\nu'} \sum_{\mu,\mu'} V_{\nu\mu\mu'\nu'}
    c^\dagger_\nu c^\dagger_\mu
    c_{\mu'} c_{\nu'}
    \tag{4.18c}
\end{align*}
$$
The Hartree mean field decoupling of $V_{\text{int}}$ is (I am being sloppy about possible $\{c, c^\dagger\} = 1$ when commuting the $c, c^\dagger$ operators)
$$
\begin{align*}
    V^H_{\text{int}}
    &= \frac{1}{2} 
    \sum_{\nu,\nu'} \sum_{\mu,\mu'}
    (c^\dagger_\nu c_{\nu'})
    (c^\dagger_\mu c_{\mu'})
    \\
    &= \frac{1}{2} 
    \sum_{\nu,\nu'} \sum_{\mu,\mu'} V_{\nu\mu\mu'\nu'}\Big[
        c^\dagger_\nu c_{\nu'}
        \langle c^\dagger_\mu c_{\mu'} \rangle
        + \langle c^\dagger_\nu c_{\nu'} \rangle
        c^\dagger_\mu c_{\mu'}
        - \langle c^\dagger_\nu c_{\nu'} \rangle
        \langle c^\dagger_\mu c_{\mu'} \rangle
    \Big]
    \tag{4.22}
\end{align*}
$$
The Fock mean field decoupling of $V_{\text{int}}$ is
$$
\begin{align*}
    V^F_{\text{int}}
    &= -\frac{1}{2} 
    \sum_{\nu,\nu'} \sum_{\mu,\mu'}
    (c^\dagger_\nu c_{\mu'})
    (c^\dagger_\mu c_{\nu'})
    \\
    &= -\frac{1}{2} 
    \sum_{\nu,\nu'} \sum_{\mu,\mu'} V_{\nu\mu\mu'\nu'}\Big[
        c^\dagger_\nu c_{\mu'}
        \langle c^\dagger_\mu c_{\nu'} \rangle
        + \langle c^\dagger_\nu c_{\mu'} \rangle
        c^\dagger_\mu c_{\nu'}
        - \langle c^\dagger_\nu c_{\mu'} \rangle
        \langle c^\dagger_\mu c_{\nu'} \rangle
    \Big]
    \tag{4.23}
\end{align*}
$$
In principle, we also have a third way of decoupling in the Cooper pairing channel, but it is often simply ignored when one is not caring about superconductivity:
$$
\begin{align*}
    V^C_{\text{int}}
    &= \frac{1}{2} 
    \sum_{\nu,\nu'} \sum_{\mu,\mu'}
    (c^\dagger_\nu c^\dagger_\mu)
    (c_{\mu'} c_{\nu'})
    \\
    &= \frac{1}{2} 
    \sum_{\nu,\nu'} \sum_{\mu,\mu'} V_{\nu\mu\mu'\nu'}\Big[
        c^\dagger_\nu c^\dagger_\mu
        \langle c_{\mu'} c_{\nu'} \rangle
        + \langle c^\dagger_\nu c^\dagger_\mu \rangle
        c_{\mu'} c_{\nu'}
        - \langle c^\dagger_\nu c^\dagger_\mu \rangle
        \langle c_{\mu'} c_{\nu'} \rangle
    \Big]
\end{align*}
$$
Then the Hartree-Fock mean field Hamiltonian, with $V^C_{\text{int}}$ omitted, is constructed as
$$
H_{\text{HF}}
= H_0 + V_{\text{int}}^H + V_{\text{int}}^F
\tag{4.24}
$$
Question 1:
Why directly adding $V_{\text{int}}^H$ and $V_{\text{int}}^F$? Why is this not a double-counting? Say, why not use $(V_{\text{int}}^H + V_{\text{int}}^F)/2$? Or even more radically, $\alpha V_{\text{int}}^H + \beta V_{\text{int}}^F$ with $\alpha + \beta = 1$?
Question 2:
By Wick's theorem, we can decompose
$$
\begin{align*}
    & c^\dagger_\nu c^\dagger_\mu
    c_{\mu'} c_{\nu'}
    \\
    &= N[c^\dagger_\nu c^\dagger_\mu c_{\mu'} c_{\nu'}]
    \\ &\quad
    + N[c^\dagger_\nu c_{\nu'}]
    c^{\dagger \bullet}_\mu c^\bullet_{\mu'}
    + c^{\dagger \bullet}_\nu c^\bullet_{\nu'}
    N[c^\dagger_\mu c_{\mu'}]
    + c^{\dagger \bullet}_\nu c^\bullet_{\nu'}
    c^{\dagger \circ}_\mu c^\circ_{\mu'}
    && \text{(Hartree)}
    \\ &\quad
    - N[c^\dagger_\nu c_{\mu'}]
    c^{\dagger \bullet}_\mu c^\bullet_{\nu'}
    - c^{\dagger \bullet}_\nu c^\bullet_{\mu'}
    N[c^\dagger_\mu c_{\nu'}]
    - c^{\dagger \bullet}_\nu c^\bullet_{\mu'}
    c^{\dagger \circ}_\mu c^\circ_{\nu'}
    && \text{(Fock)}
    \\ &\quad
    + N[c^\dagger_\nu c^\dagger_\mu]
    c^\bullet_{\mu'} c^\bullet_{\nu'}
    + c^{\dagger\bullet}_\nu c^{\dagger\bullet}_\mu
    N[c_{\mu'} c_{\nu'}]
    + c^{\dagger\bullet}_\nu c^{\dagger\bullet}_\mu
    c^\circ_{\mu'} c^\circ_{\nu'}
    && \text{(Cooper)}
\end{align*}
$$
Here $N$ is the normal-ordering symbol (all operators are at equal time, so time-ordering is omitted); contracted pairs of operators are indicated by bullets (following the notation in Wikipedia). Different ways to contract the operators correspond to the 3 channels described above. Can we rigorously relate the mean field decoupling to Wick's theorem? References that establish the relation between them are also welcome.
 A: While decoupling is physically intuitive, I find that the variational approach to mean field the most consistent mathematically. It relies on Bogoliubov's inequality and is even presented in wikipedia Mean Field.
The idea is that you want to calculate the canonical ensemble of Hamiltonian $H$ and temperature $T$, you can use its defining variational principle, namely it's the ensemble minimising the free energy:
$$
F = \langle H\rangle-TS
$$
In QM, the ensemble is typically represented by a density matrix $\rho$, so:
$$
\langle H\rangle = Tr (\rho H) \\
S = -Tr(\rho\ln\rho)
$$
and you can check that the minimum is indeed given by:
$$
\rho = \frac{1}{Z}e^{-\beta H} \\
Z = Tr(e^{-\beta H})
$$
The idea is that since you cannot do the calculation for the original $H$, you restrict yourself to a family of parametrised ensembles, and find the parameters minimising the original problem.  The condition of minimising $F$ (or practically being a stationary point) is precisely the self-consistent equation of mean field. This is the finite temperature version of the variational method for estimating the ground state of a quantum system.
In practice, these family of ensembles are canonical ensembles of quadratic Hamiltonians, which you can easily solve and use Wick's theorem (his is where it comes to play) to calculate any expected values. Note that since only $\rho$ is relevant for the variational principle, adding an overall constant to the Hamiltonian ansatz is irrelevant. Indeed, it gets cancelled out by dividing by the partition function.
Say in general that you are interested in the Hamiltonian:
$$
H = H_0+V_{int}
$$
and you replace $V_{int}$ by a simpler term $V_{mf}$ depending on various parameters $\lambda$ a parameter to be optimised using the variational principle. In the case of the form $V_{mf}=\lambda X$, the self consistent equation gives the intuitive result:
$$
\lambda = \frac{\partial \langle V_{int}\rangle_{mf} }{\partial \langle X\rangle_{mf}} \tag{1}
$$
with the average $\langle ...\rangle_{mf}$ is taken using the canonical ensemble of $H_{mf} = H_0+V_{mf}$.
Question 1
This is not double counting. The Hartree mean field is about using the ansatz:
$$
V_{mf} =     \sum_{\nu,\nu'} \sum_{\mu,\mu'} V_{\nu\mu\mu'\nu'}\left[\lambda_{\mu\mu'}^Hc^\dagger_\nu c_{\nu'} + \lambda_{\nu\nu'}^Hc^\dagger_\mu c_{\mu'}\right]
$$
with the $\lambda_{\mu\mu'}^H,\lambda_{\nu\nu'}^H$ freely varying parameters, while the Fock mean fields about using the ansatz:
$$
V_{mf} =     \sum_{\nu,\nu'} \sum_{\mu,\mu'} V_{\nu\mu\mu'\nu'}\left[\lambda_{\mu'\nu}^Fc^\dagger_\mu c_{\nu'} + \lambda_{\mu\nu'}^Fc^\dagger_\nu c_{\mu'}\right]
$$
with the $\lambda_{\mu\mu'}^F,\lambda_{\nu\nu'}^F$ freely varying parameters.
For a general $V_{\nu\mu\mu'\nu'}$, these two ansatz are different, so it is not redundant.
Question 2
Wick's theorem now comes from the fact that you want to calculate $\langle V_{int}\rangle_{mf}$. The average is obtained by doing all contractions possible, theoretically giving you all three channels. However, depending on your mean field ansatz, some of them will trivially vanish, which is why often times not all of them are accounted for.
Note that calculating the average using Wick's theorem and applying the self consistent equation $(1)$, you obtain the correct values of the parameters, namely the corresponding expected value.
Hope this helps.
Example
Including hopping and interactions, consider the following Hamiltonian:
$$
H = H_0+V_{int} \\
H_0 = \sum_x h_{xy}c_x^\dagger c_y \\
V_{int} = \frac{1}{2}\sum_{x,y} V_{xy}c^\dagger_xc^\dagger_yc_yc_x
$$
with $\xi$ real and $V_{xy}=V_{yx}$ real. Say you want to approximate by the mean field Hamiltonian:
$$
H_{mf} = H_0+V^H_{mf}+V^F_{mf} \\
V^H_{mf} = \sum_x \lambda_x c_x^\dagger c_x \\
V^F_{mf} = \sum_{x\neq y} \mu_{xy}c_x^\dagger c_y
$$
$\lambda$ real parameters and $\mu$ complex parameters satisfying $\mu_{xy} = \mu_{yx}^*$ to be determined by the variation principle. You then have:
$$
\begin{align}
\langle V_{int}\rangle_{mf} &= \sum_{x,y} V_{xy}\langle c^\dagger_xc^\dagger_y c_yc_x\rangle_{mf}  \\
&= \sum_{x,y} V_{xy}\left[\langle c^\dagger_xc_x\rangle_{mf} \langle c^\dagger_y c_y\rangle_{mf} -\langle c^\dagger_xc_y\rangle_{mf} \langle c^\dagger_y c_x\rangle_{mf}\right]  \\
\end{align}
$$
From $(1)$ I obtain:
$$
\begin{align}
\lambda_x &= \sum_{y\neq x}V_{xy}\langle c^\dagger_yc_y\rangle_{mf} \\
\mu_{xy} &= -V_{xy}\langle c^\dagger_yc_x\rangle_{mf} 
\end{align}
$$
so the mean field Hamiltonian is:
$$
H_{mf} = H_0+\sum_{y\neq x}V_{xy}\langle c^\dagger_yc_y\rangle_{mf} c^\dagger_xc_x-V_{xy}\langle c^\dagger_yc_x\rangle_{mf}  c_x^\dagger c_y
$$
which is equivalent to the decoupling method. Note that since $H_{mf}$ is defined up to an irrelevant additive constant, you don't need to add the fully contracted term.
