When can I use Wick's theorem? Wick's theorem means that for fermions, a four point correlation function (for example) can be written in terms of two point correlation functions:
\begin{equation}
\langle b_l^\dagger b_l b_m^\dagger b_m \rangle = \langle b_l^\dagger b_l \rangle \langle b_m^\dagger b_m \rangle - \langle b_l^\dagger b_m^\dagger \rangle \langle b_l b_m \rangle+  \langle b_l^\dagger b_m \rangle \langle b_l b_m^\dagger \rangle
\end{equation}
My question is when can I use this? In particular, I'm interested in finite temperature many-body perturbation theory and calculating the correlation functions from something like
\begin{equation}
Z[\overline{f},f] = \int \mathcal{D} (\overline{\phi},\phi) e^{-S_0 +S_{int}+\int_0^\beta d\tau \sum_l (\overline{f}_l b_l + b_l^\dagger f_l)},
\end{equation} 
where $S_0$ is the unperturbed part of the system, $S_{int}$ is the perturbation, and the final part of the exponential allows us to calculate the correlation functions via functional derivatives. 
Are there any circumstances under which I would need to calculate the four point correlation function itself? Or can I always use Wick's theorem?
 A: You can always use Wick's theorem when you're describing expectation values with respect to free field states (that is, non-interacting, like a single Slater determinant state). For example, in Hubbard-Stratonovich or (generally) Variational Mean-Field Theory (like Bogliubov-deGennes) you take the interaction terms and decouple them into a model that is quadratic in fermion operators, but in a background of classical fields.
In a path-integral sense then, what I'm saying is that Wick's theorem is a statement about expectation values with respect to Gaussian distributions. Since we generally can't calculate in anything other than Gaussian distributions this is handy. Consider the action
$S=ax^2+bx^4$
and we want to integrate $e^{-S}$ over all $x$. Well, this is clearly the expectation value of $e^{-bx^4}=1-bxxxx+\mathcal{O}(b^2)$ with respect to a Gaussian distribution, and Wick's trick tells me that (since $x$ is just a scalar there's no antisymmetry...): $$\langle xxxx\rangle=\langle xx\rangle+\langle xx\rangle + \langle xx\rangle$$
so to any order in $b$ we can calculate the integral and we only have to know $\langle xx\rangle$ and do a little bit of combinatorics.
The generalization of this to operators follows basically the same logic and is surprisingly straightforward (... for path integrals! Standard books always present Wick's theorem in an operator formulation that presumably doesn't even make sense at finite temperatures and I hate it.)
A: In the operator$^1$ formulation, the vital assumption, which makes the standard Wick's theorem hold, is the assumption that the contractions are in the center of the pertinent operator algebra. This is often stated casually as the contractions should be $c$-numbers, meaning that the contractions should (super)commute with all pertinent operators.

$^1$ A standard lore in physics states that the operator formalism is equivalent to the path integral formalism, although the actual map between the two formalisms may be quite subtle.
