Free Vacuum vs Interacting Vacuum and Wick's theorem I'm studying perturbation theory in QFT and I stumbled on a conceptual problem.
My understanding of the interplay between LSZ reduction formula and the Gell-Mann & Low perturbation series is that:
The LSZ formula (disregarding for a moment the normalization issue) allows to write the matrix elements of the S-matrix in terms of correlation functions, which are Vacuum expectation values of time ordered product of fields.
Namely:
$$\langle \Omega|T[\phi(x_1)\phi(x_2)...\phi(x_n)]|\Omega\rangle$$
This vacuum, $\Omega$, to my understanding is the vacuum of the interacting theory.
Such correlation functions can be found by using Gell-Mann & Low formula which relates them to vacuum expectation values of incoming fields, namely:
$$\langle \Omega|T[\phi_{IN}(x_1)\phi_{IN}(x_2)...\phi_{IN}(x_n) \exp(-i\int H_1^{IN}(z)d^4z)]|\Omega\rangle$$
Again, for how I saw this formula derived on my teacher's notes and on Bjorken Drell's book this looks to be the same vacuum as before, the vacuum of interacting theory.
Am I right? Does this matter at all?
Now comes my issue: why do we use Wick theorem to express products of $\phi_{IN}$ fields in terms of normal ordered products?
After all these are free fields and there should be no hope for the VEV to be zero on normal ordered products.
My guess, which looks to me nothing more than a desperate hand waving, is that those incoming fields are made to be equal to the real interacting fields at $-\infty$ and so it doesn't make sense to distinguish the two vacuum as eventually our interacting fields will coincide with the incoming ones at infinite times, which is where we set up the experiment and where we measure the outcome.
It seems to me that mine is a really poor understanding of the matter and I would be grateful if someone would help me to make it clear.
EDIT: I just read for the first time about Haag's theorem and I suspect my question may be tangentially related to it.
 A: I think the subscript $\text{IN}$ refers to operators in interaction picture rather than "incoming" fields (what are those anyway?).
There exists a way to morph the second expression into some interaction-picture operators sandwiched between the free theory vacuum:
$$ \langle \Omega|T[\phi(x_1)\phi(x_2)...\phi(x_n)]|\Omega\rangle = \langle \Omega|T[\phi_{IN}(x_1)\phi_{IN}(x_2)...\phi_{IN}(x_n) \exp(-i\int H_1^{IN}(z)d^4z)]|\Omega\rangle = $$
$$ = \frac{\langle 0|T[\phi_{IN}(x_1)\phi_{IN}(x_2)...\phi_{IN}(x_n) \exp(-i\int H_1^{IN}(z)d^4z)]|0\rangle}{\langle 0|T[\exp(-i\int H_1^{IN}(z)d^4z)]|0\rangle}. $$
This justifies the use of Wick's theorem. Diagrammatically it is equivalent to the following statement: we only have to consider connected Feynman diagrams (the ones which don't contain vacuum bubbles).
By the way, we could only use the free theory vacuum state in the context of interacting theory in the interaction-picture calculations. This is a mathematical trick. Physically, there is no place for a free theory state in the interacting QFT.

Since OP asked for it in the comments, I've decided to provide a detailed derivation of the relation between the interacting vacuum and the free vacuum.
Lets consider the following exponentiation of the total Hamiltonian:
$$ e^{-i\hat{H}(T+t_{0})}=\sum_{N}e^{-iE_{N}(T+t_{0})}\left|N\right>\left<N\right|. $$
Now comes a dirty trick: for large enough $T$
  all the interactions would “die out” and what remains left is the vacuum mode:
$$ \lim_{T\rightarrow\infty}e^{-i\hat{H}(T+t_{0})}\sim e^{-iE_{\Omega}(T+t_{0})}\left|\Omega\right>\left<\Omega\right|. $$
We could try to make this argument a little more mathematically precise at the cost of losing physical intuition by setting $T\rightarrow\infty(1-i\epsilon)$.
  In this case only the lowest-energy mode (the vacuum state) of the Hamiltonian operator survives exponentiation since the exponential now contains a small real part.
It could also be understood as follows: in the large $T$
  limit all the oscillations in the exponential become much more rapid than the vacuum mode, which gives the leading behavior. The last explanation provides some physical insights, but lacks mathematical precision.
Anyways, we use this trick to act with the mentioned above exponential on the free theory vacuum state $\,\left|0\right>$:
$$ \lim_{T\rightarrow\infty}e^{-i\hat{H}(T+t_{0})}\,\left|0\right>\sim e^{-iE_{\Omega}(T+t_{0})}\left<\Omega|0\right>\,\left|\Omega\right>.$$
On the other hand, we choose the free theory Hamiltonian $\hat{H}_{0}$
  to be normal-ordered, which means that the free theory vacuum energy is zero:
$$\hat{H}_{0}\left|0\right>=0;\quad\Longrightarrow\quad e^{i\hat{H}_{0}(T+t_{0})}\,\left|0\right>=\left|0\right>.$$
Therefore,
$$e^{-i\hat{H}(T+t_{0})}\left|0\right>=e^{-i\hat{H}(T+t_{0})}e^{i\hat{H}_{0}(T+t_{0})}\left|0\right>=e^{i\hat{H}((-T)-t_{0})}e^{-i\hat{H}_{0}((-T)-t_{0})}\left|0\right>=\hat{U}^{\dagger}(-T,t_{0})\,\left|0\right>=\hat{U}(t_{0},-T)\,\left|0\right>.$$
Combining these two results, we get:
$$\left|\Omega\right>\sim\frac{\hat{U}(t_{0},-T)\left|0\right>}{e^{-iE_{\Omega}(T+t_{0})}\left<\Omega|0\right>},$$
where we implicitly assume that $T\rightarrow\infty$
 . Unless the Hilbert product of the two vacuums is zero (which would mean that $\hat{H}$
  can in no way be treated as a perturbation of $\hat{H}_{0}$
 ), the denominator is just a numerical constant.
Similarly, by acting with the Hermitian conjugate exponential (in order to still take the $T\rightarrow\infty$
  limit in the end) on the bra vacuum, we obtain:
$$\left<\Omega\right|\sim\frac{\left<0\right|\,\hat{U}(T,t_{0})}{e^{-iE_{\Omega}(T-t_{0})}\left<0|\Omega\right>}.$$
The normalization condition dictates:
$$\left<\Omega|\Omega\right>=\frac{\left<0\right|\hat{U}(T,t_{0})\,\hat{U}(t_{0},-T)\left|0\right>}{e^{-iE_{\Omega}(T-t_{0})}e^{-iE_{\Omega}(T+t_{0})}\left<\Omega|0\right>\left<0|\Omega\right>}=\frac{\left<0\right|\hat{U}(T,-T)\left|0\right>}{e^{-2iE_{\Omega}T}\cdot\left|\left<0|\Omega\right>\right|^{2}}=1;$$
$$e^{-2iE_{\Omega}T}\cdot\left|\left<0|\Omega\right>\right|^{2}=\left<0\right|\hat{U}(T,-T)\left|0\right>.$$
Now we are ready to derive the expression for the correlations in the interaction picture.
$$\left<\phi_{1}(x_{1})\dots\phi_{n}(x_{n})\right>=\left<\Omega\right|\text{T}\,\hat{\phi}_{1}(x_{1})\dots\hat{\phi}_{n}(x_{n})\left|\Omega\right>=$$
$$=\frac{1}{e^{-2iE_{\Omega}T}\cdot\left|\left<0|\Omega\right>\right|^{2}}\cdot\left<0\right|\hat{U}(T,t_{0})\,\text{T}\left\{ \hat{U}(t_{0},x_{1}^{0})\,\hat{\phi}_{I\,1}(x_{1})\,\hat{U}(x_{1}^{0},t_{0})\dots\right\} \,\hat{U}(t_{0},-T)\left|0\right>.$$
We can glue together the evolution operators between interaction-picture field operators (inside the chronological ordering symbol) by using the composition law:
$$\dots\hat{\phi}_{I\,k}(x_{k})\,\hat{U}(x_{k}^{0},t_{0})\hat{U}(t_{0},x_{k+1}^{0})\,\hat{\phi}_{I\,{k+1}}(x_{k+1})\dots=\dots\hat{\phi}_{I\,k}(x_{k})\,\hat{U}(x_{k}^{0},x_{k+1}^{0})\,\hat{\phi}_{I\,{k+1}}(x_{k+1})\dots.$$
The next step would be to deal carefully with time-ordering. Each of the evolution operators is already time-ordered, because there is a time-ordered exponential in the Dyson formula. There is more: we can shuffle the operators inside the chronological ordering symbol any way we want (since they will get ordered chronologically after all). The last observation would be that since the only two evolution operators outside the chronological ordering symbol are (because of Dyson formula) also time-ordered, we could move them inside without changing anything. After all, all the evolution operators (reshuffled the way we want, inside the chronological ordering brackets) can be nicely glued together to give $\hat{U}(T,-T)$:
$$\left<\phi_{1}(x_{1})\dots\phi_{n}(x_{n})\right>=\frac{1}{N}\cdot\left<0\right|\text{T}\left\{ \hat{U}(T,-T)\,\hat{\phi}_{I1}(x_{1})\dots\hat{\phi}_{In}(x_{n})\right\} \left|0\right>,$$
where $N=e^{-2iE_{\Omega}T}\cdot\left|\left<0|\Omega\right>\right|^{2}$ is the normalization factor which we know (it was derived above) is equal to
$$N=\left<0\right|\hat{U}(T,-T)\left|0\right>=\left<0\right|\text{T}\,\hat{U}(T,-T)\left|0\right>.$$
Therefore, correlations of interacting quantum fields can be expressed formally in the interaction picture:
$$\left<\phi_{1}(x_{1})\dots\phi_{n}(x_{n})\right>=\frac{\left<0\right|\text{T}\left\{ \hat{S}\,\cdot\,\hat{\phi}_{I1}(x_{1})\dots\hat{\phi}_{In}(x_{n})\right\} \left|0\right>}{\left<0\right|\text{T}\,\hat{S}\left|0\right>},$$
where the scattering operator $\hat{S}$ is defined to be
$$\hat{S}=\lim_{T\rightarrow\infty}\,\hat{U}(T,-T)=\text{T}\,\exp\left\{ -i\intop_{-\infty}^{+\infty}dt\,\hat{H}_{I}(t)\right\} .$$
