$i\epsilon$ in the expression of Feynman Propagator

In Peskin, the Feynman's propagator for a real scalar field is first presented in a form without $i\epsilon$

$$D_F(x-y)=\int\frac{dp^3}{(2\pi)^3}\int\frac{dp^0}{2\pi}\frac{ie^{-i(x-y)}}{p^2-m^2}$$ with appropriate way of going around the poles at $p^0=E_p$ and $p^0=-E_p$ in the complex $p^0$ plane.

Subsequently, an expression for the same $D_F(x-y)$ with $i\epsilon$ which is more common is introduced

$$D_F(x-y)=\int\frac{dp^3}{(2\pi)^3}\int\frac{dp^0}{2\pi}\frac{ie^{-i(x-y)}}{p^2-m^2+i\epsilon}$$

and the poles of the integrand with $i\epsilon$ are located at $p^0=E_p-i\epsilon$ and $p^0=-E_p+i\epsilon$, with $\epsilon$ being a positive infinitesimal.

I have several questions regarding the usage of the second expression for the Feynman's propagator which seems to be standard in QFT literature.

1. In any actual computation, the residue of the second integrand at the pole $p^0=E_p-i\epsilon$ is equivalent to the residue of the first integrand at $E_p$ because $\epsilon$ is infinitesimal (similarly for the other pole) Am I correct?

2. The second expression can also be formally understood as $$D_F(x-y) = \lim_{\epsilon\to 0}\int\frac{dp^3}{(2\pi)^3}\int\frac{dp^0}{2\pi}\frac{ie^{-i(x-y)}}{p^2-m^2+i\epsilon}$$ Besides the propagtor, any expression with $i\epsilon$ in subsequent diagrammatic calculations are understood to be implicitly associated with a $\lim_{\epsilon\to 0}$ which is always suppressed in literature.

3. Some authors emphasize that the $i\epsilon$ trick displaces the poles infinitesimally and hence allow us to do the integral along the real axis in the complex $p^0$ plane. What is the importance/physical significance of being able to integrate along the real axis?

4. In my opinion, if our purpose is simply the construction of the Feynman's propagator for a wave equation/field theory, which expression to choose is only a matter of convention/taste. However, the $i\epsilon$ form seems to be the standard one in the literature of QFT. I observe that the the association of an $i\epsilon$ with the time or Hamiltonian is an essential ingredient of the arguments leading to the final expressions for correlations function or scattering amplitude in many QFT books. Is this the reason why the $i\epsilon$ expression is the preferred expression for $D_F$ in QFT literature?

The $$i\epsilon$$ is used as a prescription to tell you how to integrate in the complex $$p^{0}$$ plane. In a sense, it enforces the boundary conditions of your propagator, and uniquely fixes the contour over which you have to integrate to obtain the Green's function.

To make this more precise, note the following identity,

$$\lim_{\epsilon\rightarrow 0^+}\frac{1}{x + i\epsilon} = P\frac{1}{x} - i\pi \delta(x).$$

Now, by employing this in the expression for the Green's function you have given in your question, it is clear that the Green's function decomposes into two parts, both of which obey the formal equation $$(\partial^2 + m^2)G = 0$$. In this sense, the solution to the Green's function equation is nonunique, and the introduction of the $$i\epsilon$$ tells you the boundary conditions which fix the solution uniquely, viz. for a source supported in the past, positive energies contribute and for sources supported in the future, negative energies contribute.

This is a general overview as to why this prescription is used. The answers to the specific parts of your question are as follows.

1. Yes, in practice, once you perform the integral, the residue is indeed the one obtained for the corresponsing real value of energy, except the conditions under which you evaluate the residue are different.

2. This is a matter of convention. You could very well have taken $$-i\epsilon$$ where the infinitesimal aproaches zero from the negative axis. (I am not sure if this answers your question though.)

3. Again, the prescription gives you the boundary conditions allowing you to obtain a unique solution. Later, while evaluating loops, the $$i\epsilon$$ makes it possible to do a Wick rotation and perform the loop integrals in Euclidean space.

4. This has to do with a conceptual point. Without going into the detailed calculations, you can schematically see that the path integral has an integrand of the form $$\exp\left(i\frac{1}{2}\phi(-\partial_{\mu}\partial^{\mu} - m^2)\phi\right)$$. The $$i\epsilon$$ factor ensures the convergence of the above integrand when the path integral is performed.