$\newcommand{\ket}[1]{\left|#1\right\rangle}$
$\newcommand{\bra}[1]{\left\langle#1\right|}$
$\newcommand{\braket}[1]{\left\langle#1\right\rangle}$
After reading from Yeh's notes on Advanced Condensed Matter Theory, Section II.9, I've decided that Kleinert was right about the numerator, but the denominator actually comes from normalization.
So, you start from the definition of the Green's function.
$$G^{(n)}(x_1,\cdots, x_n)=\bra{\Omega}\mathcal{T}\phi_{1H}(x_1)\cdots\phi_{nH}(x_n)\ket{\Omega}$$
Tangent...regarding my question on the use of $\mathcal{T}$ to change pictures:
Note that $\mathcal{T}$ in the expression above is implicitly defined to time-order all of the visible symbols inside of its grasp, that is, it operates symbolically on the $\phi_{iH}(t)$. You can't use it to "cheat" like this sketchy misproof:
$$\phi_I(t)=\mathcal{T}\left[\phi_I(t)\right]\\=\mathcal{T}\left[e^{i\int_0^t dt'H_0(t')}\phi_Se^{-i\int_0^t dt'H_0(t')}\right]\\\neq \mathcal{T}\left[e^{i\int_0^t dt'H_0(t')}e^{-i\int_0^t dt'H_0(t')}\phi_S\right]\\=\mathcal{T}[\phi_S]=\phi_S$$
The middle "equality" actually fails because the mistaken proofwriter accidentally changed definitions of $\mathcal{T}$ without noticing. The first line works assuming $\mathcal{T}$ is defined to rearrange the set of explicitly visible operators (in this case, just trivially $\phi_I(t)$). But the middle equality would only be possible if the $\mathcal{T}$ had been intended to order all the embedded $H_0(t)$ as well. Since the $H_0(t)$ in line 2 do not happen to be all ordered already, one cannot make this switch of definitions for $\mathcal{T}$ without changing the value of the expression. One should thus be careful to be consistent on how $\mathcal{T}$ is defined to avoid contradictions. It doesn't operate on the value inside the brackets, it operates in a specific manner on the symbols.
Back to the flow of the proof
We want to substitute in our Gell-Mann and Low expression into the definition of the Green's function, but notice that $\ket{\Omega}$ is not just $\ket{\Psi^{\pm}}$... because $\ket{\Psi^{\pm}}$ is not normalized. So really
$$G^{(n)}(x_1,\cdots, x_n)=\frac{\bra{\Psi^+}\mathcal{T}\phi_{1H}(x_1)\cdots\phi_{nH}(x_n)\ket{\Psi^-}}{\braket{\Psi^+|\Psi^-}}$$
(where I'm using the fact that actually $\ket{\Psi^+}=\ket{\Psi^-}$ and are just different expressions for the same ket.) Expanding, we find
$$G^{(n)}(x_1,\cdots, x_n)=\frac{\bra{0}U_I(\infty,0)\mathcal{T}\left[\phi_{1H}(x_1)\cdots\phi_{nH}(x_n)\right]U_I(0,-\infty)\ket{0}}{\braket{0|U_I(\infty,0)|0}\braket{0|U_I(0,-\infty)|0}}\\\quad \times 1/\frac{\braket{0|U_I(\infty,0)U_I(0,-\infty)|0}}{\braket{0|U_I(\infty,0)|0}\braket{0|U_I(-0,\infty)|0}}$$
The unpaired Moller operators in the denominator cancel:
$$G^{(n)}(x_1,\cdots, x_n)=\frac{\bra{0}U_I(\infty,0)\mathcal{T}\left[\phi_{1H}(x_1)\cdots\phi_{nH}(x_n)\right]U_I(0,-\infty)\ket{0}}{\braket{0|U_I(\infty,-\infty)|0}}$$
What we CAN do with $\mathcal{T}$
Now we just have to clean up the numerator. Assume WLOG that, for a given $\{x_i\}$, those Heisenberg operators are time ordered when written from 1 to $n$ as above. Then write them out in terms of the interaction operators using
$$O_H(t)=U^\dagger(t,0)U_0(t,0)O_I(t)U_0^\dagger(t,0)U(t,0)$$
and
$$U_I(t, t')=U_0^{\dagger}(t, 0)U(t, t')U_0(t', 0)$$
(where subscripts indicate full, free, and interaction evolution operators) to reach
$$G^{(n)}(x_1,\cdots, x_n)=\frac{\bra{0}\mathcal{T}\left[U(\infty,t_1)\phi_{1H}(x_1)U_I(t_1,t_2)\phi_{2H}(x_2)\cdots U_I(t_{n-1},t_n)\phi_{nH}(x_n)U(t_n,-\infty)\right]\ket{0}}{\braket{0|U(\infty,-\infty)|0}}$$
I played some tricks here; let me explain. Formally, the $\mathcal{T}$ was only assumed at the start to order the $\phi_{iH}(t)$ around as unbreakable chunks (as discussed above). But since each of the individual evolution operators contains only the already time-ordered $H_I(t)$ from $t_i$ to $t_{i+1}$, then, if the $\phi_{iH}(t)$ are in order, that means all of the various implicit $H_I(t)$ in the entire expression are actually in time-order for us. So now we can actually treat the $\mathcal{T}$ as ordering the $\phi_{iH}(t)$ and the $H_I(t)$ without changing the value of our expression. Then we can also bring the outermost evolution operators inside the $\mathcal{T}$ since they are already at earlier/later times than all other operators.
Now, having played all those games with the definition of the time-ordering operator to show that we can have it apply to all of the various symbols inside the brakets for this expression...we are free to rearrange all those symbols ourselves. Let's put all the evolution operators together to the right, and voila
$$G^{(n)}(x_1,\cdots, x_n)=\frac{\bra{0}\mathcal{T}\left[\phi_{1H}(x_1)\cdots \phi_{nH}(x_n)S\right]\ket{0}}{\braket{0|S|0}}$$
where the $\mathcal{T}$ is understood to order together all of the $\phi_I(t)$ and the $H_I(t)$ inside the Dyson series for the $S$-matrix (even though the $H_I(t)$ inside the $S$-matrix are, of course, already time-ordered within themselves by the Dyson series!) Phew!
Best Answer
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\} .$$