Creo que el subíndice $\text{IN}$ se refiere a los operadores de la imagen de interacción y no a los campos "entrantes" (¿qué son esos?).
Existe una forma de transformar la segunda expresión en algunos operadores de interacción-imagen entre el vacío de la teoría libre:
$$ \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}. $$
Esto justifica el uso del teorema de Wick. Diagramáticamente es equivalente a la siguiente afirmación: sólo tenemos que considerar los diagramas de Feynman conectados (los que no contienen burbujas de vacío).
Por cierto, sólo pudimos utilizar el estado de vacío de la teoría libre en el contexto de la teoría interactiva en los cálculos de interacción-imagen. Esto es un truco matemático. Físicamente, no hay lugar para un estado de teoría libre en la QFT interactiva.
Como el OP lo pidió en los comentarios, he decidido proporcionar una derivación detallada de la relación entre el vacío interactivo y el vacío libre.
Consideremos la siguiente exponenciación del Hamiltoniano total:
$$ e^{-i\hat{H}(T+t_{0})}=\sum_{N}e^{-iE_{N}(T+t_{0})}\left|N\right>\left<N\right|. $$
Ahora viene un truco sucio: para un tamaño suficientemente grande $T$ todas las interacciones se "apagarían" y lo que queda es el modo de vacío:
$$ \lim_{T\rightarrow\infty}e^{-i\hat{H}(T+t_{0})}\sim e^{-iE_{\Omega}(T+t_{0})}\left|\Omega\right>\left<\Omega\right|. $$
Podríamos intentar hacer este argumento un poco más preciso desde el punto de vista matemático, a costa de perder la intuición física, estableciendo $T\rightarrow\infty(1-i\epsilon)$ . En este caso, sólo el modo de menor energía (el estado de vacío) del operador hamiltoniano sobrevive a la exponenciación, ya que el exponencial contiene ahora una pequeña parte real.
También podría entenderse así: en el gran $T$ límite todas las oscilaciones en el exponencial se vuelven mucho más rápidas que el modo de vacío, que da el comportamiento principal. La última explicación proporciona algunas ideas físicas, pero carece de precisión matemática.
En cualquier caso, utilizamos este truco para actuar con la mencionada exponencial sobre el estado de vacío de la teoría libre $\,\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>.$$
Por otro lado, elegimos el Hamiltoniano de la teoría libre $\hat{H}_{0}$ sea de orden normal, lo que significa que la energía de vacío de la teoría libre es cero:
$$\hat{H}_{0}\left|0\right>=0;\quad\Longrightarrow\quad e^{i\hat{H}_{0}(T+t_{0})}\,\left|0\right>=\left|0\right>.$$
Por lo tanto,
$$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>.$$
Combinando estos dos resultados, obtenemos:
$$\left|\Omega\right>\sim\frac{\hat{U}(t_{0},-T)\left|0\right>}{e^{-iE_{\Omega}(T+t_{0})}\left<\Omega|0\right>},$$
donde asumimos implícitamente que $T\rightarrow\infty$ . A menos que el producto de Hilbert de los dos vacíos sea cero (lo que significaría que $\hat{H}$ no puede ser tratada de ninguna manera como una perturbación de $\hat{H}_{0}$ ), el denominador es sólo una constante numérica.
Del mismo modo, actuando con la exponencial conjugada hermitiana (para seguir tomando la $T\rightarrow\infty$ límite al final) sobre el vacío del sujetador, obtenemos:
$$\left<\Omega\right|\sim\frac{\left<0\right|\,\hat{U}(T,t_{0})}{e^{-iE_{\Omega}(T-t_{0})}\left<0|\Omega\right>}.$$
La condición de normalización dicta:
$$\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>.$$
Ahora estamos preparados para derivar la expresión de las correlaciones en la imagen de la interacción.
$$\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>.$$
Podemos pegar los operadores de evolución entre los operadores de campo de interacción-imagen (dentro del símbolo de ordenación cronológica) utilizando la ley de composición:
$$\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.$$
El siguiente paso sería tratar con cuidado la ordenación del tiempo. Cada uno de los operadores de evolución ya está ordenado en el tiempo, porque hay un exponencial ordenado en el tiempo en la fórmula de Dyson. Hay más: podemos barajar los operadores dentro del símbolo de ordenación cronológica como queramos (ya que al fin y al cabo se ordenarán cronológicamente). La última observación sería que, dado que los dos únicos operadores de evolución que están fuera del símbolo de ordenación cronológica también están (por la fórmula de Dyson) ordenados en el tiempo, podríamos moverlos dentro sin cambiar nada. Al fin y al cabo, todos los operadores de evolución (reordenados como queremos, dentro de los corchetes de ordenación cronológica) pueden pegarse muy bien para dar $\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>,$$
donde $N=e^{-2iE_{\Omega}T}\cdot\left|\left<0|\Omega\right>\right|^{2}$ es el factor de normalización que sabemos (fue derivado anteriormente) es igual a
$$N=\left<0\right|\hat{U}(T,-T)\left|0\right>=\left<0\right|\text{T}\,\hat{U}(T,-T)\left|0\right>.$$
Por lo tanto, las correlaciones de los campos cuánticos que interactúan pueden expresarse formalmente en la imagen de la interacción:
$$\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>},$$
donde el operador de dispersión $\hat{S}$ se define como
$$\hat{S}=\lim_{T\rightarrow\infty}\,\hat{U}(T,-T)=\text{T}\,\exp\left\{ -i\intop_{-\infty}^{+\infty}dt\,\hat{H}_{I}(t)\right\} .$$
