Solution (source code)

= Solution

Write $b_k=\langle b_k\rangle+\delta b_k$. The <mean-field approximation> drops the quadratic product of the fluctuations:
$$
b_k^\dagger b_{k'}\simeq b_k^\dagger\langle b_{k'}\rangle+\langle b_k^\dagger\rangle b_{k'}-\langle b_k^\dagger\rangle\langle b_{k'}\rangle.
$$
For a Hermitian pairing interaction, the two linear terms are Hermitian conjugates. Defining the <BCS anomalous average> and pairing field by $\Delta_k=-\sum_{k'}V_{k,k'}\langle b_{k'}\rangle$, and writing $\xi_k=\epsilon_k-\mu$, gives
$$
\boxed{K\equiv H-\mu N=\sum_{k,\sigma}\xi_kc_{k\sigma}^\dagger c_{k\sigma}-\sum_k(\Delta_kb_k^\dagger+\Delta_k^*b_k)-\sum_{k,k'}V_{k,k'}\langle b_k^\dagger\rangle\langle b_{k'}\rangle}.
$$
The final constant corrects the double counting in the <reduced BCS pairing Hamiltonian>.