There are too many Wick's Theorems! This is a follow-up question to QMechanic's great answer in this question. They give a formulation of Wick's theorem as a purely combinatoric statement relating two total orders $\mathcal T$ and $\colon \cdots \colon$ on an algebra.

I have come across "Wick's theorems" in many contexts. While some of them are special cases of the theorem [1], others are -- as far as I can see -- not. I am wondering if there is an even more general framework in which Wick's theorem can be presented, showing that all of these theorems are in fact the same combinatoric statement.


*

*Wick's theorem applies to a string of creation and annihilation operators, as described e.g. on Wikipedia:
$$ ABCD = \mathopen{\colon} ABCD \mathclose{\colon} + \sum_{\text{singles}} \mathopen{\colon} A^\bullet B^\bullet CD \mathclose{\colon} + \cdots \tag{*} $$
Here, the left hand side is "unordered" and it seems to me that [1] is not valid?

*The creation and annihilation operators in (*) can be either bosonic or fermionic.
This technicality is not a problem in [1] since it allows for graded algebras.

*Wick's theorem can also be applied to field operators:
$$ \mathcal T\, \phi_1 \cdots \phi_N = \mathopen{\colon} \phi_1 \cdots \phi_N \mathclose{\colon} + \sum_{\text{singles}} \mathopen{\colon} \phi_1^\bullet \phi_2^\bullet \cdots \phi_N \mathclose{\colon} + \cdots $$
Since the mode expansion of a field operator $\phi_k$ consists of annihilation and creation operators, normal ordering is actually not simply a total order on the algebra of field operators. Once again, we can not apply [1]?

*In a class I am taking right now, we applied Wick's theorem like this to field operators that didn't depend on time:
$$ \phi_1 \cdots \phi_N = \mathopen{\colon} \phi_1 \cdots \phi_N \mathclose{\colon} + \sum_{\text{singles}} \mathopen{\colon} \phi_1^\bullet \phi_2^\bullet \cdots \phi_N \mathclose{\colon} + \cdots $$
This seems to combine the issues of points 1 and 3...

*In probability theory, there is Isserlis' Theorem:
$$ \mathbb E(X_1 \cdots X_{2N}) = \sum_{\text{Wick}} \prod \mathbb E(X_i X_j) $$
This looks like it should also be a consequence from one and the same theorem, but I don't even know what the algebra would be here.

*My string theory lectures were quite a while ago, but I vaguely remember that there we had radial ordering instead of time ordering. Also there seems to be some connection to OPEs.
This seems to not be a problem with [1].

*In thermal field theory, the definition of normal ordering changes.
This seems to not be a problem with [1] either.
 A: Various comments to the post (v3):


*

*One may speculate that seemingly unordered operators are in practice always ordered wrt. some order. 

*-

*As long as the fields $\phi_i=\phi_i^{(+)}+\phi_i^{(-)}$ are linear in creation and annihilation operators, this should not be a problem. 

*-

*Isserlis' theorem is related to the path-integral formulation of Wick's theorem, cf. e.g. this Phys.SE post.

*-

*-
The single-most important generalization of the operator formulation of Wick's theorem (as compared to my Phys.SE answer) is to consider contractions that doesn't belong to the algebra center. This is often used in CFT, see e.g. Ref. 1.
References:


*

*J. Fuchs, Affine Lie Algebras and Quantum Groups, (1992); eq. (3.1.35). 

A: I wanted to expand on Wick's theorem form the point of view of 5. This is from the point of view of Euclidean free path integrals. I think this perspective is very illuminating but certain points of it are not emphasized in the literature. The discussion is more straightforward when thinking of a finite-dimensional analogue. In here we take our space of fields to be the finite-dimensional space $V=\mathbb{R}^n$. We can put on this space the linear coordinates $\{\phi^i|i\in\{1,\dots,n\}\}$. To clarify the physical intuition, one should think of the index $i$ as a position in a discrete spacetime with $n$ points.
In this setting a free theory is determined by (unnormalized) correlation function obtained via a Gaussian path integral, which in this case is just a finite-dimensional integral. Observables are determined by polynomial functions $F(\phi)$ and the correlation functions are of the form
$$\langle F(\phi)\rangle=\int\text{d}^n\phi\, e^{-\frac{1}{2}\phi^iA_{ij}\phi^j}F(\phi),$$
for a symmetric and positive-definite $A_{ij}$.
Wick's theorem in version 5 can be easily proved following the discussion in https://arxiv.org/abs/1202.1554. This is obtained by noting that the integral of a total derivative vanishes since the exponential decays at the boundary due to the positive-definiteness of $A_{ij}$. Indeed, by the product rule
$$0=\int\text{d}^n\phi\,\frac{\partial}{\partial\phi^i}\left(e^{-\frac{1}{2}\phi^iA_{jk}\phi^j}\phi^{r_1}\cdots\phi{r_s}\right)=-A_{ij}\langle\phi^j\phi^{r_1}\cdots\phi^{r_s}\rangle+\sum_{t=1}^s\delta^{r_l}_i\langle \phi^{r_1}\cdots\widehat{\phi^{r_l}}\cdots\phi^{r_s}\rangle,$$
where the $\widehat{\phi^{r_l}}$ means we skip this term. Denoting by $A^{ij}$ the inverse matrix $A^{ij}A_{jk}=\delta^i_k$, we can solve this equation to
$$\langle \phi^i\phi^{r_1}\cdots\phi^{r_s}\rangle=\sum_{t=1}^sA^{ir_l}\langle \phi^{r_1}\cdots\widehat{\phi^{r_l}}\cdots\phi^{r_s}\rangle.$$
This is Wick's theorem! It says that in order to compute the correlation function we just need to consider all possible contractions of $\phi^i$ with all other terms, each contraction taking out a propagator $A^{ir_l}=\langle \phi^i\phi^{r_l}\rangle$. Then a simple induction shows that
$$\langle \phi^{r_1}\cdots\phi^{r_s}\rangle=\sum_{P\in\text{Pair}(s)}\prod_{\{a,b\}\in P}\langle\phi^a\phi^b\rangle$$
Now, the other Wick's theorems can be obtained from this one as follows. First, we need to define the notion of normal ordering in this setting. This definition is particularly physical. Let $F(\phi)$ be a monomial in $\phi$ like $F(\phi)=\phi^{i_1}\cdots\phi^{i_u}$. We define the normal ordering $:F(\phi):$ to be the polynomial such that all correlations $\langle:F(\phi):G(\phi)\rangle$ for a polynomial $G(\phi)$ is obtained by considering all Wick contractions contributing to $\langle F(\phi)G(\phi)\rangle$ except those with contractions of two fields within the monomial $F(\phi)$.
From this definition it is not clear that such a polynomial exists or, if it does, if it is unique. Uniqueness should be some consequence of a theorem saying that a polynomial is completely determined by moments. In any case, to proof existence one can give an explicit construction. Uniqueness is more or less clear from it.
For the normal ordering of a bilineal monomial the construction is clear from Wick's theorem
$$\langle \phi^i\phi^j\phi^{r_1}\cdots\phi^{r_s}\rangle=\langle\phi^i\phi^j\rangle\langle \phi^{r_1}\cdots\phi^{r_s}\rangle+\sum_{t=1}^s\langle\phi^i\phi^{r_l}\rangle\langle \phi^j\phi^{r_1}\cdots\widehat{\phi^{r_l}}\cdots\phi^{r_s}\rangle.$$
The correlation $$\langle:\phi^i\phi^j:\phi^{r_1}\cdots\phi^{r_s}\rangle$$ should only consist of the last term. Then it is clear what to do, define
$$:\phi^i\phi^j:=\phi^i\phi^j-\langle\phi^i\phi^j\rangle.$$
One can repeat this for higher order monomials but I won't do it here since the computations do get a bit complicated.
In general, we have Wick's theorem
$$:\phi^{i_1}\cdots\phi^{i_u}:=\phi^{i_1}\cdots\phi^{i_u}-\sum_{\{a,b\}}\langle\phi^{i_a}\phi^{i_b}\rangle:\phi^{i_1}\cdots\widehat{\phi^{i_a}}\cdots\widehat{\phi^{i_b}}\cdots\phi^{i_u}:-\sum_{\{a,b\},\{c,d\}}\langle\phi^{i_a}\phi^{i_b}\rangle\langle\phi^{i_c}\phi^{i_d}\rangle:\phi^{i_1}\cdots\widehat{\phi^{i_a}}\cdots\widehat{\phi^{i_b}}\cdots\widehat{\phi^{i_c}}\cdots\widehat{\phi^{i_d}}\cdots\phi^{i_u}:-\cdots$$
where the first sum is over 1-contractions, the second term is over 2-contractions, and so on. Although the combinatorics of the proof can get a little messy, the big picture is rather simple. The terms of the form $\langle\phi^{i_a}\phi^{i_b}\rangle:\phi^{i_1}\cdots\widehat{\phi^{i_a}}\cdots\widehat{\phi^{i_b}}\cdots\phi^{i_u}:$ appearing in the first sum are the ones cancelling all contributions to correlation functions containing a single Wick contraction within $\phi^{i_1}\cdots\phi^{i_u}$. Similarly, the terms of the form $\langle\phi^{i_a}\phi^{i_b}\rangle\langle\phi^{i_c}\phi^{i_d}\rangle:\phi^{i_1}\cdots\widehat{\phi^{i_a}}\cdots\widehat{\phi^{i_b}}\cdots\widehat{\phi^{i_c}}\cdots\widehat{\phi^{i_d}}\cdots\phi^{i_u}:$ cancel correlation functions containing only two Wick contractions within $\phi^{i_1}\cdots\phi^{i_u}$. This is the form of Wick's theorem appearing in version 4. It gives an explicit inductive formula for normal ordering.
Let me comment now on version 3. In our setting we have defined normal ordering through its behaviour in correlation functions. These are computed by path integrals, which automatically time order. This means that in the operator formalism these correspond to matrix elements of a time ordered operator $\mathcal{T}\hat{\phi}^{i_1}\cdots\hat{\phi}^{i_u}$. So version 4 of Wick's theorem corresponds to version 3, the former being in the path integral formalism while the second in the operator formalism.
In order to go from version 4 to version 5, one just needs to note that ⟨:():⟩=0. Indeed, in order to obtain a non-zero answer one needs to add at least a monomial of degree equal to that of (). Only then will one start having contractions that don’t pair any two elements within (). Incidentally, this also clarifies the relationship with the creation/annihilation statement, since the normal ordering there precisely annihilates the vacuum expectation values by placing annihilation operators to the right. More precisely, it can be seen that creation/annihilation normal ordering for a product of two fields (linear in creation and annihilation operators) is also given by
$$:\phi^i\phi^j:=\mathcal{T}\phi^i\phi^j-\langle\phi^i\phi^j\rangle.$$
This normal ordering also satisfies the recursion relation imposed by Wick's theorem to obtain normal ordering of higher order monomials in the fields. We conclude that both normal orderings coincide on bilinears and satisfy the same recursion relation. They must then coincide always.
The OPEs can also be understood from this point of view of the path integral formalism. The main idea of the free case though is the following. To compute the operator product expansion of a group of operators, we would like to express them as a series of well-defined operators at a single point in spacetime weighted by a coefficient functions depending on the positions of the original operators which may diverge as these positions get close to one another. Being well-defined just means that its correlation functions with other operators far away are all convergent. This is most easily done writing the product of the operators using Wick's theorem. This is because the divergent parts appear inside correlation functions and are thus numerical coefficients. All other operators appear inside normal ordering and thus, when inserted into correlation functions, are never contracted with one another. There are thus no divergences when computing correlation functions with far away operators.
The discussion above is made clear with an example. Consider the operator product expansion of $\phi(0)\phi(x)$ in a free scalar field theory. One could try to write this a series of operators at $0$ by Taylor expanding
$$\phi(0)\phi(x)=\phi(0)^2+x\phi(0)\partial\phi(0)+\frac{1}{2}x^2\phi(0)\partial^2\phi(0)+\cdots.$$ However, in this series all operators are ill-defined. For example,
$$\langle\phi(0)^2\phi(x)\phi(y)\rangle=\langle\phi(0)^2\rangle\langle\phi(x)\phi(y)\rangle+2\langle\phi(0)\phi(x)\rangle\langle\phi(0)\phi(y)\rangle$$
and this term diverges, even when $x$ and $y$ are away from each other and $0$. On the other hand, we can Taylor expand after using Wick's theorem
$$\phi(0)\phi(x)=:\phi(0)\phi(x):+\langle{\phi(0)\phi(x)}\rangle=\langle{\phi(0)\phi(x)}\rangle+:\phi(0)^2:+x:\phi(0)\partial\phi(0):+\frac{1}{2}x^2:\phi(0)\partial^2\phi(0):+\cdots.$$
This is precisely in the OPE form. The first term is a numerical function that diverges as $x\rightarrow 0$ multiplied by a well defined operator, the identity operator. The rest of the terms are well defined operators as well. For example now
$$\langle:\phi(0)^2:\phi(x)\phi(y)\rangle=2\langle\phi(0)\phi(x)\rangle\langle\phi(0)\phi(y)\rangle,$$
which is well-defined as long as $x$, $y$, and $0$ are apart from one another. In particular, we see that the divergent part of the OPE is
$$\phi(0)\phi(x)\sim\langle{\phi(0)\phi(x)}\rangle,$$
which is heavily used in free field theory.
This procedure can be expanded to the interacting case using perturbation theory. For definiteness, let me explain this using $\phi^4$ theory. In perturbation theory we have
$$\langle\phi(0)\phi(x)\cdots\rangle=\int\mathcal{D}\phi e^{-\frac{1}{2\hbar}\int\text{d}^D y\phi(-\Delta)\phi+\frac{\lambda}{4!\hbar}\int\text{d}^Dy\phi^4}\phi(0)\phi(x)\cdots=\sum_{n=0}^\infty\frac{\lambda^n}{4!^n\hbar^n}\int\text{d}^Dy_1\cdots\text{d}^Dy_n\int\mathcal{D}\phi e^{-\frac{1}{2\hbar}\int\text{d}^D y\phi(-\Delta)\phi}\phi(0)\phi(x)\phi(y_1)^4\cdots\phi(y_n)^4\cdots=\sum_{n=0}^\infty\frac{\lambda^n}{4!^n\hbar^n}\int\text{d}^Dy_1\cdots\text{d}^Dy_n\langle\phi(0)\phi(x)\phi(y_1)^4\cdots\phi(y_n)^4\cdots\rangle_G$$
The subscript $G$ indicates that the last correlation is taken in the free theory. Accordingly, we can apply Wick's theorem on each of these terms individually.
For the $n=0$ term, we have the contributions to the operator product expansion
$$\langle\phi(0)\phi(x)\cdots\rangle_G=\langle:\phi(0)\phi(x):\cdots\rangle_G+\langle\phi(0)\phi(x)\rangle\langle\cdots\rangle_G$$
For the first term, we can do a Taylor expansion around $x=0$, just like we did on the free case. This yields the $\lambda^0$ contribution to the OPE. In terms of $\hbar$, the first term contributes at order $\hbar^0$ while the second at order $\hbar$. Only the second has divergent terms as $x\rightarrow 0$. Moreover, we can use Feynman diagrams to keep track of these

As we see, in these diagrams all external legs are automatically normal ordered, so that it is understood that in full correlation functions they shouldn't be contracted with one another. In particular, we can expand in a Taylor series when these legs are close to 0.
Now, let us consider the expansion via Wick's theorem of the order $\lambda$ term
$$\frac{\lambda}{4!\hbar}\int\text{d}^Dy\langle\phi(0)\phi(x)\phi(y)^4\cdots\rangle_G$$.
The first term has no contractions
$$\frac{\lambda}{4!\hbar}\int\text{d}^Dy\langle:\phi(0)\phi(x)\phi(y)^4:\cdots\rangle_G.$$
We can represent this with the following Feynman diagram

As before, all external legs are normal ordered. We also see that the external legs coming from the vertex do not carry propagators. This vertex contributes at order $\hbar^{-1}$ and has no divergences as $x\rightarrow 0$.
There are 4 terms coming from having 1 contraction, which are represented by the Feynman diagrams

All of these contribute at order $\hbar^0$ and only the first diverges as $x\rightarrow 0$. However, this divergence is in a certain way already captured from a term at order $\lambda^0$. In fact, we can resum all of the terms with this divergence
$$\langle\phi(0)\phi(x)\rangle_G\langle:e^{\frac{\lambda}{4!\hbar}\int\text{d}^Dy\phi^4}:\cdots\rangle_G$$
The terms with two contractions are of the from

All of these contribute at order $\hbar$ and only the first (possibly) diverges as $x\rightarrow 0$ (well, the second one also diverges but we have already discussed these type of terms above). This term is really interesting and it is thoroughly explored in https://pirsa.org/18030064. There is is shown that it does diverge in $D=4$, and in fact, its divergence is of the form
$$\frac{\lambda}{2\hbar}\int\text{d}^Dy\langle\phi(0)\phi(y)\rangle_G\langle\phi(x)\phi(y)\rangle_G:\phi(0)^2:,$$
when expanding the external legs around $0$.
Finally, we have the terms with 3 contractions

these all contribute at order $\hbar^2$ but only the second has a new divergence. This divergence multiplies the identity operator.
In summary, for OPEs in the interacting case we sum up over diagrams of the type above. Disconnected diagrams either don't have any divergences as $x\rightarrow 0$ (if there is no path connecting the $\phi(0)$ and $\phi(x)$ vertices), or their divergences already appear in a connected diagram of a lower order in perturbation theory. As a final comment, all of these diagrams also suffer from loop divergences that have to be renormalized as usual in perturbative quantum field theory.
A: I will give an answer to explain why there are too many Wick's theorems in condensed matter physics or many-body physics.
Actually, the importance of Wick's theorem is closely related to the calculation of Green's function. Green's function techniques in condensed matter physics or in many-body physics usually rely on expansion of the Green's function in question (generally contains quartic terms in Hamiltonian) in an infinite series of higher Green's functions for a noninteracting solvable system and a subsequent contraction into products of one-particle Green's function. This decomposition is greatly simplified by the use of suggestive diagrammatic representations. The rigorous foundation of this procedure is known as Wick's theorem.


*

*The first meet


We first meet Wick's theorem is to formulate the many-body perturbation expansion of zero-temperature Green's function in which the problem can be described by Hamiltonian:
$$H=H_0+H_i$$ 
where $H_i$ is the complex many-body interaction.


*

*The second meet


We will meet Wick's theorem again when we perform the many-body expansion of the finite temeprature Green's function in which the problem can also be described by Hamiltonian $H=H_0+H_i$. The big difference compared to zero temperature Green's function is that the system is no longer in a ground state instead of a mixed state by the density matrix
$$\rho = \dfrac{e^{-\beta H}}{Tr[e^{-\beta H}]}.$$
One can see the equilibrium many-body density matrix also contains many-body interactions. To formulate the simultaneous expansion on both density matrix and the time evolution operator:
$$U(t)=e^{-i H t/\hbar}$$ 
Matsubara's strategy: replace $\tau=it$ and treat $\tau$ as a real number. As a result of this replacement, many-body perturbation expansion becomes possible.


*

*The third meet


Keldysh formalism: which is suitable for the investigation of nonequilibrium many-body problem. (Here the Wick's theorem is much like zero temperature one.)


*

*$\cdots \cdots$


The following links are the recommended literature to prove Wick's theorem and discuss the interrelations between different versions of Wick's theorem.
1.Wick theorem for general initial states;
2.Equilibrium and nonequilibrium many-body perturbation theory;
