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\} .$$
Best Answer
As Lubos Motl mentions in a comment, for all practical purposes, OP's sought-for eq. (1) is proved via Wick's Theorem.
It is interesting to try to generalize Wick's Theorem and to try to minimize the number of assumptions that goes into it. Here we will outline one possible approach.
I) Assume that a family $(\hat{A}_i)_{i\in I}$ of operators $\hat{A}_i\in{\cal A}$ lives in a (super) operator algebra ${\cal A}$
with (super) commutator $[\cdot,\cdot]$, and
with center $Z({\cal A})$.
Here
the index $i\in I$ runs over an index set $I$ (it could be continuous), and
the index $i$ contains information, such as, e.g., position $x$, time instant $t$, annihilation/creation label, type of field, etc., of the operator $\hat{A}_i$.
II) Assume that $$ \forall i,j\in I~: \qquad [\hat{A}_i,\hat{A}_j] ~\in~Z({\cal A}). $$
III) Assume that there are given two ordering prescriptions, say $T$ and $::$. Here $T$ and $::$ could in principle denote any two ordering prescriptions, e.g. time order, normal order, radial order, Weyl order$^1$, etc. This means that the index set $I$ is endowed with two strict total orders, say, $<$ and $\prec$, respectively, such that
The $T$ symbol is (graded) multilinear wrt. supernumbers.
$ T(\hat{A}_{\pi(i_1)}\ldots\hat{A}_{\pi(i_n)})~=~(-1)^{\sigma_{\pi}} T(\hat{A}_{i_1}\ldots\hat{A}_{i_n} )$ is (graded) symmetric, where $\pi\in S_n$ is a permutation of $n$ elements, and $(-1)^{\sigma_{\pi}}$ is a Koszul sign factor.$^2$
$ T(\hat{A}_{i_1}\ldots\hat{A}_{i_n} )~=~\hat{A}_{i_1}\ldots\hat{A}_{i_n}$ if $i_1 > \ldots > i_n$.
In the special case where some of the $ i_1 , \ldots , i_n$ are equal$^3$ (wrt. the order <), then one should symmetrize in appropriate (graded) sense over the corresponding subsets. For instance, $$ T(\hat{A}_{i_1}\ldots\hat{A}_{i_n} )~=~\hat{A}_{i_1}\ldots\hat{A}_{i_{k-1}}\frac{\hat{A}_{i_k}\hat{A}_{i_{k+1}}+(-1)^{|\hat{A}{i_k}||\hat{A}{i_{k+1}}|}\hat{A}_{i_{k+1}}\hat{A}_{i_k}}{2}\hat{A}_{i_{k+2}}\ldots\hat{A}_{i_n}$$ if $i_1 > \ldots > i_k=i_{k+1}> \ldots > i_n$.
[Similar conditions 1-4 should hold for the second ordering $(::,\prec)$.]
IV) It then follows from assumptions I-III that the (generalized) contractions $$ \hat{C}_{ij}~=~T(\hat{A}_i\hat{A}_j)~-~:\hat{A}_i\hat{A}_j:~\in~Z({\cal A}) $$ belong to the center $Z({\cal A})$. The contractions are graded symmetric $$ \hat{C}_{ij}~=~(-1)^{|\hat{A}_i||\hat{A}_j|} \hat{C}_{ji}. $$
V) Assume furthermore that the contractions $\hat{C}_{ij}$ do not depend on the operators $\hat{A}_k$, i.e. $$ \frac{\partial \hat{C}_{ij}}{\partial \hat{A}_k}~=~0 $$ in order to simplify combinatoric arguments below.
VI) It is now a straightforward exercise to establish the corresponding Wick's Theorem $$ T(f(\hat{A})) ~=~ \exp\left(\frac{1}{2}\sum_{i,j\in I}\hat{C}_{ij}\frac{\partial}{\partial\hat{A}_j}\frac{\partial}{\partial\hat{A}_i} \right):f(\hat{A}):, $$ meaning a rule for how to re-express one ordering prescription $T(f(\hat{A}))$ [where $f$ is a sufficiently nice function of the $(\hat{A}_i)_{i\in I}$ family] in terms of the other ordering prescription $::$ and (multiple) contractions $\hat{C}_{ij}$. And vice-versa with the roles of the two orderings $T$ and $::$ interchanged: $$ :f(\hat{A}): ~=~ \exp\left(-\frac{1}{2}\sum_{i,j\in I}\hat{C}_{ij}\frac{\partial}{\partial\hat{A}_j}\frac{\partial}{\partial\hat{A}_i} \right)T(f(\hat{A})). $$ Such Wick's Theorems can now be applied successively to establish nested Wick's Theorems, such as, e.g.,$^4$ $$ T(:f(\hat{A})::g(\hat{A}):) ~=~ \left. \exp\left(\sum_{i,j\in I}\hat{C}_{ij}\frac{\partial}{\partial\hat{A}_j}\frac{\partial}{\partial\hat{B}_i} \right) :f(\hat{A}) g(\hat{B}): \right|_{\hat{B}=\hat{A}}. $$ These Wick's Theorems may be extended to a larger class of operators than just the $(\hat{A}_i)_{i\in I}$ family through (graded) multilinearity.
VII) Let us now assume that the operators $\hat{A}_i$ are Bosonic for simplicity. A particular consequence of a nested Wick's Theorem is the following version
$$T(:\hat{A}^2_i::\hat{A}^2_j:) ~=~ 2\hat{C}_{ij}^2 + 4 \hat{C}_{ij}:\hat{A}_i\hat{A}_j: + :\hat{A}^2_i\hat{A}^2_j:$$
of OP's sought-for eq. (1). Finally, let us mention that Wick's Theorem, radial order, OPE, etc., are also discussed in this and this Phys.SE posts.
--
Footnotes:
$^1$ Example: The Weyl/symmetric ordering satisfies $$W(f(\hat{A})) ~=~\left. \exp\left(\sum_{i\in I}\hat{A}_i \frac{\partial}{\partial a_i} \right) f(a) \right|_{a=0}. $$ For more details, see e.g. my Phys.SE answer here.
$^2$ The Koszul sign convention produces a minus sign every time two Grassmann-odd objects are permuted. In this answer $|\hat{A}_i|=0,1 \pmod 2$ denotes the Grassmann-parity of $\hat{A}_i$.
$^3$ Being equal wrt. an order is in general an equivalence relation, and it is often a weaker condition than being equal as elements of $I$.
$^4$ A nested Wick's Theorem (between radial order and normal order) is briefly stated in eq. (2.2.10) on p. 39 in J. Polchinski, String Theory, Vol. 1. Beware that radial order is often only implicitly written in CFT texts. By the way, a side-effect/peculiarity of nested ordering symbols are discussed in this Phys.SE post.