Where is the pole prescription in quantum mechanics?
The pole prescription comes about when you solve time-dependent wave equations via Fourier transforms. In quantum mechanics that would take the time-dependent Schrödinger equation to the time-independent Schrödinger equation. However the corresponding time-dependent Green's functions do not converge when you Fourier transform them, which is how the $i\epsilon$ prescription comes about. For more detail see my answer to this question, which was also linked by the OP as related.
Are there four independent Green's functions in quantum mechanics or just two?
In quantum mechanics, a time-dependent Green function is a solution of
$$(i \frac{\partial}{\partial t} - H ) G(t) = \mathbb{I}\delta(t) .$$
The solution is not unique, therefore there are as many Green's functions as boundary conditions (BCs) you can think up. A particularly useful one is the BC
$$ G^+(t) = 0 \quad \textrm{for } t<0,$$
since it gives a solution to the initial value problem via
$$\psi(t) = iG^+(t-t_0)\psi(t_0) \quad \textrm{for } t>t_0.$$
Note that this is exactly the same in relativistic quantum mechanics/quantum field theory, only that the corresponding wave equation differs and you have to Fourier transform a few more dimensions. E.g. for the Klein-Gordon equation, which was used as an example in the question, we have
$$ (\partial^2 + m^2) \psi(t) = 0,$$
such that
$$(\partial^2 + m^2) G(t) = \mathbb{I}\delta(t).$$
Imposing the same boundary condition as above and manipulating a little bit gives the Green function in Fourier space as
$$ G(p) \propto \frac{1}{p^2 - m^2 + i\epsilon}.$$
Note that this is unique now (details), the sign of $i\epsilon$ has implicitly been chosen by the boundary condition of the time-dependent Green function. I can never remember if that is called the advanced/retarded/Feynman Green's function and I think the terms also differ in the literature (e.g. in scattering theory you would call this the retarded Green function, but in QFT it seems to be the Feynman Green function). Either way it is the one that is useful for solving initial value problems.
Is the quantum mechanical path integral secretly related to the Feynman Green's function, or is it indifferent?
Other people will know more about this than me. Hope the answer helps anyway.
Any solution to the differential equation for the Green's function can
be used. If you solve for the Green's function using contour integration
over the frequencies, the solutions for the portion of the contour
near the pole are solutions of the homogeneous wave equation. Choosing
different contours corresponds to adding different solutions to the
homogeneous wave equation and therefore give different boundary conditions.
You could also find these solutions using other methods to
solve the differential equation without contour integration.
The reason to choose various boundary conditions is that when you use
Green's theorem to write the desired solution, you have terms where
you integrate over the source, and terms which come from the
space-time surface. If you can choose the boundary conditions for your
Green's function so that the surface terms are zero, your calculation
becomes simpler.
Every Green's function gives one particular solution.
Choosing the
retarded Green's function will make the time surface terms zero if the
fields are all zero before the source turns on. Often this is the solution
you want. Other boundary conditions solve other cases directly.
Best Answer
The $\epsilon$ in (3.59) does not shift the poles but rather regulates the integral over 3-momentum.
More concretely, we start with (note the extra $i$ in fig. 3) $$ iD^+(t,\vec x) = \frac{1}{(2\pi)^4}\int_{\gamma^+} dp_0\int \ d^3\vec p \,\frac{e^{-ip_0t+i\vec p\cdot \vec x}}{p_0^2-|\vec p|^2}, $$ where $\gamma^+$ is a contour that encircle only the pole in the right half-plane. By the residue theorem, we find $$ iD^+(t,\vec x) =\frac{i}{(2\pi)^3} \int \ d^3\vec p\, \frac{e^{-i|\vec p|t+i\vec p\cdot \vec x}}{2|\vec p|}. $$ The $\epsilon$ now enters the game to regulate the above integral as $t \rightarrow t-i\epsilon$ such that it is convergent. The rest of the details are relatively straightforward $$ \begin{align*} iD^+(t,\vec x) &=\frac{i}{(2\pi)^2} \int_{-1}^1d(\cos\theta)\int _0^\infty d|\vec p|\, |\vec p|^2 \frac{e^{-|\vec p|(\epsilon+it)+i|\vec p| |\vec x| \cos \theta}}{2|\vec p|}\\ &=\frac{i}{8\pi^2}\frac{1}{i|\vec x|}\int _0^\infty d|\vec p|\,e^{-|\vec p|(\epsilon+it)}(e^{i|\vec p||\vec x|}-e^{-i|\vec p||\vec x|})\\ &=\frac{i}{8\pi^2}\frac{1}{i|\vec x|}\left(\frac{1}{\epsilon+it-i|\vec x|}-\frac{1}{\epsilon+it+i|\vec x|}\right)\\ &=\frac{i}{8\pi^2}\frac{1}{i|\vec x|}\frac{2i|\vec x|}{(\epsilon+it)^2+|\vec x|^2} = -\frac{i}{4\pi^2}\frac{1}{(t-i\epsilon)^2-|\vec x|^2} \end{align*} $$ The regulator can be "removed" using the Sokhotski–Plemelj theorem as explained in (3.60).