Evaluating propagator without the epsilon trick Consider the Klein–Gordon equation and its propagator:
$$G(x,y) = \frac{1}{(2\pi)^4}\int d^4 p \frac{e^{-i p.(x-y)}}{p^2 - m^2} \; .$$
I'd like to see a method of evaluating explicit form of $G$ which does not involve  avoiding singularities by the $\varepsilon$ trick. Can you provide such a method?
 A: Expanding on dmckee's comment:  
The $+i\epsilon$-trick has the blessing of OCD mathematicians because it follows directly from a deep fact about the group of spacetime translations:  the group $\{e^{-i\langle P,x\rangle/\hbar}| x \in \mathbb{R}^n\}$ of spacetime translations is the boundary of an analytic semigroup $\{e^{-i\langle P,\xi\rangle/\hbar}| x \in \mathbb{C}^n \mbox{ and } Im(\xi) \leq 0\}$.  
Many quantities in field theory are expressed in terms of these translations, and frequently these quantities can be computed more easily by analytically continuing from real "Minkowski" time to imagininary "Euclidean" time, where the delicate cancellation of phases becomes the crude suppression of exponential damping.  When you use the $+i\epsilon$-trick, what you're really doing is saying that the particular cancellation of phases you want is the one which respects this analyticity. This is precisely what's happening when you use the $+i\epsilon$-trick to evaluate the Klein-Gordon propagator.  You've got an integral which does not converge absolutely, and you're picking out a certain resummation which does.  The $+i\epsilon$ is not just a trick here; it's really the definition of the quantity you're after.
A: Before answering the question more or less directly, I'd like to point out that this is a good question that provides an object lesson and opens a foray into the topics of singular integral equations, analytic continuation and dispersion relations. Here are some references of these more advanced topics: Muskhelishvili, Singular Integral Equations; Courant & Hilbert, Methods of Mathematical Physics, Vol I, Ch 3; Dispersion Theory in High Energy Physics, Queen & Violini; Eden et.al., The Analytic S-matrix. There is also a condensed discussion of `invariant functions' in Schweber, An Intro to Relativistic QFT Ch13d.
The quick answer is that, for $m^2 \in\mathbb{R}$, there's no "shortcut." One must choose a path around the singularities in the denominator. The appropriate choice is governed by the boundary conditions of the problem at hand. The $+i\epsilon$ "trick" (it's not a "trick") simply encodes the boundary conditions relevant for causal propagation of particles and antiparticles in field theory.
We briefly study the analytic form of $G(x-y;m)$ to demonstrate some of these features.
Note, first, that for real values of $p^2$, the singularity in the denominator of the integrand signals the presence of (a) branch point(s). In fact, [Huang, Quantum Field Theory: From Operators to Path Integrals, p29] the Feynman propagator for the scalar field (your equation) may be explicitly evaluated:
\begin{align}
G(x-y;m) &= \lim_{\epsilon \to 0} \frac{1}{(2 \pi)^4} \int d^4p \, \frac{e^{-ip\cdot(x-y)}}{p^2 -  m^2 + i\epsilon}  \nonumber \\
&= \left \{ \begin{matrix}
-\frac{1}{4 \pi} \delta(s) + \frac{m}{8 \pi \sqrt{s}} H_1^{(1)}(m \sqrt{s}) & \textrm{ if }\, s \geq 0 \\
 -\frac{i m}{ 4 \pi^2 \sqrt{-s}} K_1(m \sqrt{-s}) & \textrm{if }\, s < 0.
\end{matrix} \right.
\end{align}
where $s=(x-y)^2$.
The first-order Hankel function of the first kind $H^{(1)}_1$ has a logarithmic branch point at $x=0$; so does the modified Bessel function of the second kind, $K_1$. (Look at the small $x$ behavior of these functions to see this.) 
A branch point indicates that the Cauchy-Riemann conditions have broken down at $x=0$ (or $z=x+iy=0$). And the fact that these singularities are logarithmic is an indication that we have an endpoint singularity [eg. Eden et. al., Ch 2.1]. (To see this, consider $m=0$, then the integrand, $p^{-2}$, has a zero at the lower limit of integration in $dp^2$.)
Coming back to the question of boundary conditions, there is a good discussion in Sakurai, Advanced Quantum Mechanics, Ch4.4 [NB: "East Coast" metric]. You can see that for large values of $s>0$ from the above expression that we have an outgoing wave from the asymptotic form of the Hankel function. 
Connecting it back to the original references I cited above, the $+i\epsilon$ form is a version of the Plemelj formula [Muskhelishvili]. And the expression for the propagator is a type of Cauchy integral [Musk.; Eden et.al.]. And this notions lead quickly to the topics I mentioned above -- certainly a rich landscape for research.
A: As far as my experience goes, the problem stems from writing the right solution for all reals to the problem: 
$$
(p-m)G(p)=1.
$$
which reads:
$$
G(p)=\text{P.v.} \frac {1}{p-m}+c_0\delta(p-m)
$$
where $\text{P.v.}$ stands for principal value. The $\delta(\epsilon-\omega)$ function appears as it is the Kernel of $(\omega-\epsilon)$ and $c_0$ is some constant to be fixed. If we now take the Fourier transform we obtain:
$$
\int e^{ipt }G(p)=\:\left( i\pi  \text{sign}(t)+c_0\right)e^{i m t}
$$
$c_0$ must now be fixed according to boundary conditions; for the retarded and advanced  Green functions, one has $c_0=\pm i\pi$ and the solution given by the $i \epsilon$ trick is recovered. In my opinion though, it is a rather bad method as it only works when you have poles or first order, as the $\delta$ function can be approached by square integrable functions. If you now are looking for the solution of: 
$$
(p-m)^kG(p)=1.
$$
with k an integer, you now have
$$
G(p)=\text{P.v.} \frac {1}{(p-m)^k}+\sum_{j=0}^kc_j\delta^{(j)}(p-m)
$$
with $\delta^{(k)}$ the k-th derivative of the delta function. The Fourier transform reads
$$
G(t)=\left( i\pi \frac{(i t)^{k-1}}{k-1!}\text{sign}(t)+\sum_{j=0}^kc_j (-it)^j\right)e^{i m t}
$$
And again, the $c_j$s are fixed depending on boundary conditions. Yet I don't know any way to recover this solution with the $i\epsilon$ trick.
