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
The first question we have to ask is: what is a one particle state in an interacting theory? It is reasonable to require that they are states that are both momentum eigenstates and energy eigenstates. (In fact, as the Hamiltonian and the momentum operator commute, these are not two different conditions.) Weinberg, in his famous textbook, says that particle states are those which transform under an irreducible representation of the Poincare group, but we need not fuss around with the Poincare group here.
All we will say is that, in the interacting theory, there are some single particle states, labelled by
$$|\lambda k \rangle$$
where $k$ is the four-momentum, and $\lambda$ is whatever other labels we need for our particles. (In this answer I will be working with just a real scalar field, but even in the spin-0 case there can still be extra data that distinguishes our particles in an interacting theory.)
Now, we know know we have a set of momentum and energy eigenstates $|\lambda k\rangle$ that represent the stable particles of our theory. We can now "smudge out" these definite momentum states into wave packets, using a Gaussian window function $f_W$ that has some momentum uncertainty $\kappa$. We will denote these smudged out, approximate energy and momentum eigenstates with a subscript $W$ for "window."
$$|\lambda k\rangle_W \equiv \int d^3{\mathbf{k}'} f_W(\mathbf{k} - \mathbf{k}') |\lambda k\rangle$$
We will come back to these.
Now, the free vacuum $|0\rangle$ of $\hat H_0$ and the true vacuum $|\Omega\rangle$ of $\hat H = \hat H_0 + \hat H_{\rm int}$ are very different states. Particles in the interacting theory must indeed be defined to be formed from the action of the "creation operator" on the true vacuum, as long as we properly define what we mean by the "creation operator" in the interacting theory.
To create an annihilate particles, we will use the Klein Gordon inner product. (We suppress $\hbar$ and $c$.)
$$(\psi_1, \psi_2)_{KG} \equiv i \int d^3 x (\psi^*_1 \partial_t \psi_2 - \partial_t \psi^*_1 \psi_2)$$
The motivation for defining this is that in the FREE theory, the Klein Gordon inner product gives us an inner product between single particle states. If we have two single particle states (in the free theory) $|\Psi_1\rangle$ and $|\Psi_2 \rangle$, we have
$$\langle \Psi_1 | \Psi_2 \rangle = ( \psi_1, \psi_2 )_{KG}$$
where we used the "single particle wave functions" of the states defined by
$$\psi_i(x) \equiv \langle 0| \hat \phi(x) |\Psi_i\rangle$$
The niceties of free field theory come from the simple algebra of the creation and annihilation operators, combined with the fact that the annihilation operator annihilates the vacuum. We will try to recreate those relationships using the Klein Gordon inner product. However, to do this, we will need to use widely separated wave packets.
From here on out, everything will be in the interacting theory.
For a given function $\psi$, we define the creation and annihilation operators that "create" the state corresponding to that wave function as follows. $$ \hat a^\dagger_i (t) \equiv -\big( \psi^*_i(t, \cdot), \hat \phi(t, \cdot) \big)_{KG} $$ $$\hat a_i(t) = \big( \psi_i(t, \cdot), \hat \phi(t, \cdot) \big)_{KG}$$
(In the free theory, this creation operator literally would create the single particle state with the single particle wave function $\psi_1$.)
(Something I must mention about these operators is their time evolution. It is a point of notational confusion that $\hat a^\dagger_{1}(t)$ depends explicitly on a time $t$, given that we usually have defined time dependence such that $e^{i \hat H t} \hat{O}(t') e^{-i \hat H t} = \hat {\mathcal{O}}(t'+t)$. This is not the case here.)
Now, sadly, in the interacting theory, the annihilation operator defined above will not annihilate the vacuum. However, we can recover something close:
$$\langle \Omega| \hat a_1(t) |\Omega\rangle = i \int d^3{x}\langle \Omega| \big( \psi_1^*(t, \vec x) \partial_t \hat \phi (t, \vec x) - \partial_t \psi_1^*(t, \vec x) \hat \phi(t, \vec x) \big) |\Omega\rangle$$ $$= i \int d^3{x} \big( \psi_1^*(t, \vec x) \partial_t \langle \Omega| \hat \phi(t, \vec x) |\Omega\rangle - \partial_t \psi^*_1(t, \vec x) \langle \Omega| \hat \phi(t, \vec x) |\Omega\rangle \big)$$ $$= i \langle \Omega| \hat \phi(t, \vec x) |\Omega\rangle \int d^3{x} (- \partial_t \psi_1^*(t, \vec x))$$
The fact that $\partial_t \langle \Omega| \hat \phi(t, \vec x) |\Omega\rangle = 0$ follows directly from the fact that the vacuum state has zero energy, so $e^{- i \hat H t} |\Omega\rangle = |\Omega\rangle$. Now as we want $\langle \Omega| \hat a_1(t) |\Omega\rangle = 0$ for any $\psi_1$, we can see that this is achieved if and only if $\langle \Omega| \hat \phi(x) |\Omega\rangle = \langle \Omega| \hat \phi(0) |\Omega\rangle = 0$. We will assume this is the case.
In the free theory, $\langle 0| \hat a_1(t) \hat a_2^\dagger(t) |0\rangle = \langle \Psi_1 | \Psi_2\rangle = (\psi_1, \psi_2)_{KG}$. In an interacting theory, for any $\hat a_1$ and state $|\Psi_2\rangle$ (not just a single particle state) we have
$$\langle \Omega| \hat a_1(t) |\Psi_2\rangle = \langle \Omega| \big( \psi_1(t, \cdot) , \hat \phi(t, \cdot) \big)_{KG}|\Psi_2\rangle$$ $$= \big( \psi_1(t, \cdot), \langle \Omega| \hat \phi(t, \cdot) |\Psi_2\rangle\big)_{KG} $$ $$= \big( \psi_1(t, \cdot), \psi_2(t, \cdot) \big)_{KG}$$
$$\langle\Psi_2| \hat a_1(t) |\Omega\rangle = \big( \psi_1(t, \cdot) , \psi_2^*(t, \cdot) \big)_{KG}$$
Remember our single particle states? We're now going consider the "single particle wave function" of those states. Namely, they have to be plane waves.
$$\langle \Omega| \hat \phi(x) |\lambda k\rangle = C_\lambda e^{-ikx}$$
where $C_\lambda$ is a constant that depends on $\lambda$.
We now want to see what our states $\hat a^\dagger_1|\Omega\rangle$ have to do with these true particle wave packets $| \lambda k \rangle_W$. To do this, we will see what the inner product of these two states are. Just from our simple algebra above, for an annihilation operator $\hat a_{\lambda_1 k_1} = (\psi_{k_1}, \hat \phi)_{KG}$ where $k_1^2 = m_{\lambda_1}^2$, we have
\begin{equation*} \begin{split} \langle \Omega | \hat a_{\lambda_1 k_1} (t) |\lambda_2 k_2\rangle_W = \big( \psi_{ k_1}(t, \cdot), \langle \Omega |\hat \phi(t, \cdot)|\lambda_2 k_2\rangle_W \big)_{KG} = C_{\lambda_2}\big(\psi_{ k_1}(t, \cdot), \psi_{ k_2 }(t, \cdot) \big)_{KG}\\ {}_W \langle \lambda_2 k_2 | \hat a_{\lambda_1 k_1}(t) |\Omega\rangle = \big( \psi_{ k_1}(t, \cdot), {}_W \langle \lambda_2 k_2 |\hat \phi(t, \cdot) |\Omega\rangle \big)_{KG} = C_{\lambda_2}\big(\psi_{ k_1}(t, \cdot), \psi^*_{ k_2 }(t, \cdot) \big)_{KG}. \end{split} \end{equation*} We desire for the top expression to be $\propto \delta_{\lambda_1 \lambda_2} \delta^3(\mathbf{k}_1 - \mathbf{k}_2)$ and for the bottom expression to be 0. If this were the case, then the only single particle state $\hat a^\dagger_{ k_1}(t) |\Omega\rangle$ would overlap with would be $|\lambda_1 k_1\rangle$, and $\hat a_{k_1}(t)$ could still functionally "annihilate" the vacuum, even though we need to keep ${}_W \langle \lambda k |$ on the left. Defining $\omega_{\lambda k} \equiv (m^2_{\lambda} + \mathbf{k}^2)^\frac{1}{2}$, we have
\begin{equation*} \begin{split} \big(\psi_{ k_1}(t, \cdot), \psi_{ k_2 }(t, \cdot) \big)_{KG} = (2 \pi)^3 \int d^3{\mathbf{k}} f_W(\mathbf{k}_1 - \mathbf{k}) f_W(\mathbf{k}_2 - \mathbf{k}) (\omega_{\lambda_1 k} + \omega_{\lambda_2 k})e^{it(\omega_{\lambda_1 k} - \omega_{\lambda_2 k})} \\ \big(\psi_{ k_1}(t, \cdot), \psi_{ k_2 }^*(t, \cdot) \big)_{KG} = (2 \pi)^3 \int d^3{\mathbf{k}} f_W(\mathbf{k}_1 - \mathbf{k}) f_W(\mathbf{k}_2 + \mathbf{k}) (\omega_{\lambda_1 k} - \omega_{\lambda_2 k})e^{it(\omega_{\lambda_1 k} + \omega_{\lambda_2 k})}. \\ \end{split} \end{equation*}
The top expression is not $\propto \delta_{\lambda_1 \lambda_2} \delta^3_W(\mathbf{k}_1 - \mathbf{k}_2)$ and the bottom expression is not $0$. However, if we take $\kappa \ll |\mathbf{k}_1 - \mathbf{k}_2|$ and also take $t \to \pm \infty$, they are! This hinges on our assumption that $m_{\lambda_1} \neq m_{\lambda_2}$ if $\lambda_1 \neq \lambda_2$. The $e^{it (\ldots)}$ term will oscillate wildly in both integrals if $\lambda_1 \neq \lambda_2$, causing them to be 0. In the top integral, this oscillation does not occur when $\lambda_1 = \lambda_2$. Furthermore, the top integral will be negligible unless $\mathbf{k}_1 = \mathbf{k}_2$. Taking the $f_W(\mathbf{k}) \to \delta^3(\mathbf{k})$ and $t \to \pm \infty$ limit, we can now write
\begin{equation*} \begin{split} \langle \lambda_2 k_2 | \hat a_{\lambda_1 k_1}^\dagger (\pm \infty) |\Omega\rangle = C_{\lambda_2} (2 \pi)^3 2 \omega_{\lambda_2 k_2} \delta_{\lambda_1 \lambda_2} \delta^3(\mathbf{k}_1 - \mathbf{k}_2) \\ \langle \lambda_2 k_2 | \hat a_{\lambda_1 k_1} (\pm \infty) |\Omega\rangle = 0. \end{split} \end{equation*} These properties are even more important than I let on. This is because the states $| \lambda k \rangle$ are so generally defined: they are just momentum eigenstates with all the extra necessary data stuffed into $\lambda$. As they diagonalize the momentum operator, they form a basis of our entire state space! Therefore, we can immediately see from the first equation that
$$\hat a^\dagger_{\lambda k}(\pm \infty) |\Omega\rangle = -C_\lambda \big( e^{ikx}, \hat \phi( x)\big)_{KG} |\Omega\rangle \vert_{t = \pm \infty} = | \lambda k \rangle$$ where we have chosen the normalization $\langle\lambda k | \lambda' k' \rangle = C_{\lambda}^* C_{\lambda'} (2 \pi)^3 (2 \omega_{\lambda k}) \delta_{\lambda \lambda'} \delta^3(\mathbf{k} - \mathbf{k}')$. From the second equation, we can immediately see that
$$\langle\Psi | \hat a_{\lambda k}(\pm \infty) |\Omega\rangle = 0 \hspace{0.15 cm} \text{ for all } \langle\Psi | \hspace{0.5 cm} \Longrightarrow \hspace{0.5 cm} \hat a_{\lambda k}(\pm \infty) |\Omega\rangle = 0.$$ Apparently our asymptotic creation and annihilation operators behave almost exactly like our good old creation and annihilation operators from the free theory!
There's another important property I must mention, which is that two creation/annihilation operators that have different $\lambda k$ data will commute. This is a direct consequence of the fact that our creation/annihilation operators are spacial integrals weighted by wave packets that are spatially separated at large times. (For operators with the same $k$ but different $\lambda$, as $m_\lambda$ is different the wave packets will propagate at different speeds and will still succeed to separate.) Note that spatial separation is a property of wave packets but not of plane waves. This is another place where it is necessary to view plane waves as a limit of wave packets in order to properly understand your theory. In fact, the operators will not commute unless they are defined with this limiting procedure.
We are finally ready to define our incoming and outgoing multi-particle states. As our asymptotic creation operators only change the ground state in localized spatial regions and each spatial excitation is justifiably called a "particle state" we can say that acting with a few of them on the ground state will create a perfectly good multi-particle state. We will now define our incoming (created at $t = -\infty$) and outgoing (created at $t = + \infty$) multi-particle asymptotic states.
$$ |\lambda_1 k_1, \ldots, \lambda_n k_n\rangle_{\rm in} \equiv \hat a^\dagger_{\lambda_1 k_1}(-\infty) \ldots \hat a^\dagger_{\lambda_n k_n}(-\infty) |\Omega\rangle \\ |\lambda_1 k_1, \ldots, \lambda_n k_n\rangle_{\rm out} \equiv \hat a^\dagger_{\lambda_1 k_1}(+\infty) \ldots \hat a^\dagger_{\lambda_n k_n}(+\infty) |\Omega\rangle$$
The four-momenta $k_i$ will have masses $k_i^2 = m_{\lambda_i}^2$ and no $|\lambda_i k_i\rangle$ is allowed to equal another. Some people prefer to rescale $\hat \phi$ in order to hide those $C_\lambda$ prefactors but I will not. The nature of these prefactors will be explored much later. It is important to note that the total momenta of these states are approximately the sum of all $\mathbf{k}_i$, and the energy is approximately the sum of all $\omega_{\lambda_i k_i}$. This lends more credence to the notion that these are "multi-particle" states.
Now that we have successfully defined our incoming and outgoing asymptotic multi particle states and derived some important properties of our newly constructed asymptotic creation and annihilation operators, we have completed the framework necessary to derive the LSZ reduction formula. Using the properties defined here, you should be able to justifiably go through the steps as outlined in Srednicki.
To answer your doubt 2: In order to get our states to have the right properties, we needed these to be wave-packets that are widely separated in the distant past and future. Therefore, these states are only approximately momentum and energy eigenstates (although you can get as close as you want). As they're not perfect energy eigenstates, some time evolution will occur. You particles will start far apart, come together, interact, then (different) particles will leave.
TLDR: If you define creation and annihilation operators properly, using the Klein Gordon inner product with widely separated wave packets in the far past/future, you will get your actual particle states when acting with these operators on the true vacuum $|\Omega\rangle$.