Propagator of a real scalar field does not give an unambiguous result The propagator of a real scalar field is basically a Green function. Its evaluation requires specifying the contour which seems arbitrary. The integral depends on the choice of the contour. For four different choices for the contour, the integral converges to four different values. 


*

*Does physics tell us what should be the correct choice of the contour?

*Why does the integral converges to different values for different choices of the contour even though we are evaluating the same integral? I apologize if this is a mathematics question.
 A: When we integrate the propagator with respect to $k^0$ (i.e. the energy), we encounter two poles: one at $\omega_{\mathbf{k}}=\sqrt{\mathbf{k}^2+m^2}$ and one at $-\omega_{\mathbf{k}}=-\sqrt{\mathbf{k}^2+m^2}$, where $\mathbf{k}$ is the 3-momentum. In order to regularize this integral, we move these poles slightly off of the real line by adding or subtracting $i\epsilon$ and taking the limit as $\epsilon\to 0$. For each pole, there are two ways to do this (adding vs. subtracting $i\epsilon$), which leads to four total choices of regularized propagator:


*

*The causal advanced propagator, given by choosing $\omega_{\mathbf{k}}+i\epsilon$ and $-\omega_{\mathbf{k}}+i\epsilon$,

*The causal retarded propagator, given by choosing $\omega_{\mathbf{k}}-i\epsilon$ and $-\omega_{\mathbf{k}}-i\epsilon$,

*The time-ordered propagator, given by choosing $\omega_{\mathbf{k}}-i\epsilon$ and $-\omega_{\mathbf{k}}+i\epsilon$, and

*The anti-time-ordered propagator, given by choosing $\omega_{\mathbf{k}}+i\epsilon$ and $-\omega_{\mathbf{k}}-i\epsilon$.
It is only the time-ordered propagator (also known as the Feynman propagator) that gives you precisely the time-ordered two-point correlation function of two free scalar fields. In a similar vein, the anti-time-ordered propagator gives you the anti-time-ordered two-point correlation function of two free scalar fields.
The other two choices for regularization (that is, the causal propagators), as their names suggest, give propagators that vanish unless the one point is in the future light cone (for the retarded propagator) or past light cone (for the advanced propagator) of the other point. In other words, these correlation functions vanish at spacelike separation.
The reason the Feynman propagator is the one most commonly applied is because we are usually interested in its product, namely the time-ordered two-point correlation function. As to why we're interested in this quantity in particular, that's the subject of the first few chapters of most quantum field theory textbooks. In short, the time-ordered correlation function is the one that is most amenable to perturbative expansion.

For your second question, the integral changes value depending on the regularization procedure because of how contour integrals are typically done in complex analysis. In almost every case, extensive use is made of Cauchy's theorem and the related residue theorem, which state, in essence: integrating a holomorphic function around a closed curve in the complex plane will give you either zero or a specific quantity related to how many and what kind of singularities the curve encloses. When we want to integrate a function over the real line, we do this by integrating in a giant semicircle and extending the diameter of the semicircle to infinity*. Thus, we enclose either all poles in the upper half-plane or all poles in the lower half-plane (and yes, the choice of semicircle direction is part of the regularization procedure), and where we move the poles during regularization determines whether they're enclosed by our semicircle.
*The reason this procedure gives us the integral on the real line is due to Jordan's Lemma, which states, in essence, that the integral along the "circle" part of the semicircular contour vanishes as the semicircle's diameter goes to infinity, leaving only the "diameter" part along the real line.
A: The different choices of contour yield Green's functions with different boundary conditions. (Remember that Green's functions aren't unique, because different functions can satisfy the same differential equation. The Green's function is only defined up to the kernel of the differential operator.)
Strictly speaking, you're not actually evaluating the same integral for different choices of contour. If the poles lie exactly on the real axis then the integral doesn't converge, so you need to move them slightly off the real axis in order to get a convergent integral. Different choices of directions to move the poles yield slightly different integrands. The physical situation tells you which one is physically correct, but the answer depends on what problem you're trying to solve. The Feynman propagator is sometimes, but not always, the physically correct one.
A: This is not as complete of an answer as I could give but the detais aren't that hard to fill in. 
The thing is that the propagator is a very specific Green's function of the Klein-Gordon equation. This one
$$D(x-y)=\langle0|T\Phi(x)\Phi(y)|0\rangle$$
Where $|0\rangle$ is the vacuum and the $T$ means time ordering. You should commute the inner fields such that the coordinate with later time appears first. 
Instead of manipulating the contour, we can make the right poles appear inside the contour by adding a small positive number $\epsilon$ to the integrand 
$$D(x-y)=\int\frac{d^4k}{(2\pi)^4}\frac{i}{k^2-m^2+i\epsilon}e^{-ik(x-y)}$$
And integrate over $k^0$ in a semicircle for which the integral converges (this depends on the sign of $(x-y)^0$, taking care of the time ordering). 
The contour matters because the integrand has divergences in the complex plane. There are two of them, and depending on which one you keep inside the contour the complex integral takes on a different value. This is called Cauchy's integral theorem. 
A: May be  I misunderstand your question but I try to answer it. Consider Feynman propagator, which solves EOM with delta function, $-i\delta^{(d)}(x-y)$ in right hand side. This propagator has integral representation,
$$G_F=\int_p\frac{e^{-ip\cdot(x-y)}}{p^2-m^2-i\epsilon},$$
where $\epsilon$-prescription says how one should consider poles of denominator (and $\epsilon\rightarrow 0$) and we work in $d=(D+1)$-space time and
$$\int_p=\int\frac{d\omega_p}{2\pi}\int\frac{d^Dp}{(2\pi)^D}.$$
First, $p^2-m^2-i\epsilon=\omega_p^2-(m^2+{\bf p}^2)-i\epsilon$ and
$$\omega_p^2-(m^2+{\bf p}^2)-i\epsilon=(\omega_p-\sqrt{{\bf p}^2+m^2}+i\epsilon)(\omega_p+\sqrt{{\bf p}^2+m^2}-i\epsilon),$$
where we have used $\epsilon\rightarrow 0$. Then, integral over $\omega_p$ becomes
$$I=\oint\frac{d\omega_p}{2\pi}\frac{ie^{-i\omega_p(x^0-y^0)}}{(\omega_p-\sqrt{{\bf p}^2+m^2}+i\epsilon)(\omega_p+\sqrt{{\bf p}^2+m^2}-i\epsilon)}$$
In case of $t=x^0-y^0>0$, we close contour in upper-half plane (I denote thi contour as $\mathcal{K}$). For $t=x^0-y^0<0$ we close contour in lower-half plane (this contour is $\mathcal{K}'$). For $\mathcal{K}$-contour, only one pole gives contribution. You can easily calculate residue at this pole and find that
$$I=\frac{e^{-it\sqrt{{\bf p}^2+m^2}}}{2\sqrt{{\bf p}^2+m^2}}.$$
For the $\mathcal{K}'$-contour the answer is the same up to differnt sign in exp ($-$ instead of $+$).
Let me denote obtained integral (for contour $\mathcal{K}$) as
$$G_F=\int\frac{d^D p}{(2\pi)^D}e^{i{\bf p}\cdot({\bf x}-{\bf y})}I\equiv D({\bf x}-{\bf y},t).$$
In case of $\mathcal{K}'$-contour, you obtain $D({\bf x}-{\bf y},-t)$. In addition, we know that due to Lorentz invariance $D({\bf r},t)=D(-{\bf r},-t)$. Finally, you'll find that $$G_F=\theta(x^0-y^0)D(x-y)+\theta(y^0-x^0)D(x-y)=\langle 0| T\phi(x)\phi(y)|0\rangle$$
What physics says about it? In my opinion (may be I am wrong), physics says us about casual properties of this propagator. It is quiet physical when you think about contours in terms of casuality. Also, you explicitly have four different functions, but it is not unambiguous because this functions solve different equations or have different properties:


*

*Wightman propagator, $D=\langle 0|\phi(x)\phi(y)|0\rangle$, which solves EOM with zero in RHS

*Feynman propagator, $G_F=\langle 0| T\phi(x)\phi(y)|0\rangle$, which solves EOM with $-i\delta^{(d)}(x-y)$ in RHS (and anti-time ordering version)

*Advanced propagator, $G_A$, which exists only for $x^0>y^0$

*Retarded propagator, $G_R$
