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.
