Consider the case of a free scalar field, governed by the usual Lagrangian,
$$\mathcal{L} = \frac{1}{2}\partial_\mu \phi \partial^\mu \phi - \frac{1}{2}m^2 \phi^2$$
The propagator, or equivalently Green's function for the theory is a function which can be though of as a response when we use a delta function as an input in the equations of motion, i.e.
$$\square \Delta_F(x-y) \sim \delta^{(4)}(x-y)$$
Explicitly, it is given by a Fourier integral over four-momentum,
$$\Delta_F (x-y)=\int \frac{d^4 p}{(2\pi)^4} \frac{i}{p^2-m^2} e^{-ip \cdot (x-y)}$$
We encounter a singularity at $p^0 = \pm \sqrt{p^2+m^2} = \pm E_{p}$. We can choose a contour which avoids these by dipping below the first, and then above the other. However, to apply the residue theorem, it must be closed. For $x^0 > y^0$, we close it in an anti-clockwise direction in the upper-half plane, and the opposite if $y^0 > x^0$. Alternatively, we can define the Feynman propagator,
$$\Delta_F = \int \frac{d^4 p}{(2\pi)^4} \, \frac{i}{p^2-m^2 + i\epsilon} e^{-ip\cdot (x-y)}$$
The $i\epsilon$ prescription due to Feynman shifts the poles by an infinitesimal amount away from the real axis; as a result the integral going straight through the real line is equivalent to the aforementioned contour.
Depending on your purpose, it may be useful to pick a particular contour, in which case we can define the retarded propagator $\Delta_R$ as the one which chooses to go over each pole on the real line, and the advanced contour going under. See the depiction below:
To understand what they mean physically, consider in the context of response theory, the response function $\chi$ which determines how a system changes under the addition of a source $\phi$, i.e.
$$\delta \langle \mathcal{O}_i(t) \rangle = \int dt' \chi(t-t')\phi(t')$$
It's clear the aforementioned is in fact a convolution, $\chi \ast \phi$, and $\chi$ also has the interpretation of a Green's function. But we can't affect the past, so clearly,
$$\chi(t) = 0, \quad t < 0$$
For the Fourier transform, this is equivalent to stating,
$$\chi(\omega) \quad \mathrm{analytic}, \quad \mathrm{Im} \, \omega > 0$$
In other words, $\chi(\omega)$ has no poles in the upper-half plane. This object $\chi(t)$ is in fact our retarded Green's function, and is also called the causal Green's function, precisely because the above requirement is imposed.
the integral
$$ \int_0^{\infty} e^{i\epsilon t}dt$$
is not-defined for $\epsilon \in R$. For your integral, Mathematica probably also gave you the condition that ${\rm Im}\{\epsilon_\lambda - \epsilon\} < 0$. For real frequencies, we have to add a small infinitesimal in order to assure convergence. Thus
$$ \int_0^{\infty} dt e^{i\epsilon t} \to \int_0^{\infty} dt e^{i(\epsilon+i\gamma) t} = -\frac{i}{\epsilon+i\gamma}$$
the physical meaning of this $\gamma$ is to give a finite life-time to correlations in the system. The limit $\gamma \to 0^{+}$ assures us that particles and holes (in non-interacting systems) have the limit of infinite life-times as they do not decay. Sometimes you will see that this $\gamma$ is replaced with some finite inverse life-time in a phenomenological style $\gamma \to \tau^{-1}$.
Also: the imaginary part of the Green function is associated with the density of states. Adding this infinitesimal and taking the limit of $\gamma \to 0^{+}$ give us a delta function, as desired. Upon replacing $\gamma \to \tau^{-1}$, in this context, the density of states attains a finite width.
Best Answer
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.