I am looking to gain a more rigorous and deeper understanding as to how an $i\varepsilon$ prescription actually changes the end result of a divergent integral, specifically in regards to Green's functions. I understand the use of $i\varepsilon$ allows one to impose desired boundary conditions. A very helpful Phys.SE answer here suggests that we should think of the $i\varepsilon$ as a contribution due to damping, which all real systems possess. At the end, we take the limit of $\varepsilon \to 0$. Yet I am puzzled by this fact. It seems that if we do not inlcude $i\varepsilon$ at the offset, we will not obtain the same answer (see example below). Both the principle value and the $i\varepsilon$ method both serve as valid Green's functions (with different B.C.'s I suppose), but I am under the impression that the $i\varepsilon$ is fundamentally what we are after.
The closest answer that I have come across comes from Shankar's Prinicples of Quantum Mechanics, which in regards to inverting the Helmholtz operator, he says we must consider a slightly different different operator $\nabla^2 +k^2 +i\varepsilon$, where the eigenfunctions are plane waves of complex wavenumber. He goes on to say that such functions are not part of the space we are restricting oursleves to, namely, the space of functions normalized to unity or the Dirac delta function. Thus $\nabla^2 +k^2 +i\varepsilon$ may be inverted within the physical Hilbert space. And then $i\varepsilon$ is sent to zero at the end of the calculation.
I am curious as to how, by working with functions outside our restricted space, we may thus appropriately invert an operator? I believe this sort of generalization to complex operators (and their associated space they act on) is on the right track.
Here is a basic example on the inequivalence between $i\varepsilon$ and the Cauchy principle value. I have found this result extends to the case with two poles on the real axis.
If we start with a function such as $$ f(x) = \int_{-\infty}^\infty dk\ \frac{e^{ik x}}{k}\quad\text{for}\quad x>0 $$ and seek to assign a non-divergent value, we may find the principle value by integrating over a contour of $$\text{P.V}[ f(x)] + \int_\delta + \int_{arc} = \oint. $$ Here $\text{P.V}$ is the principle value, $\delta$ represents a small semicircle contribution around the pole, and $arc$ is the contribution from the arc at $+\infty$. With Jordan's lemma, we close above and thus the $arc$ at infinity makes no contribution. We are then left with $$\text{P.V}[ f(x)] = \oint - \int_\delta = \pm i\pi a_{-1} +2\pi i\sum_{res}.$$ Here $a_{-1}$ represents the first term of a Laurent expansion (this is contribution from $\int_\delta$). This term is positive (negative) for a clockwise (counter-clockwise) semi-circle of vanishing radius $\delta$.
To see that principle value is unique, we can make a detour around the pole either clockwise (in which we do enclose the pole) or counter clockwise (no pole enclosed, no residue). Either way $$\text{P.V}[ f(x)] = i\pi.$$
Okay, now lets use the idea of pushing the pole upwards or downwards by $$k \to k\pm i\varepsilon$$ and take the limit as $\varepsilon \to 0$. We may now carry the integration across the real line unhindered and we find, when closing above as we must, $$ \lim_{\varepsilon\to 0} \int_{-\infty}^\infty dk\ \frac{e^{i(k\pm i\varepsilon) x}}{k\pm i\varepsilon} = \begin{cases} 0&\quad\text{for}\quad k \to k+ i\varepsilon \\ 2\pi i&\quad\text{for}\quad k \to k-i\varepsilon \end{cases} $$ This last part of obtaining $0$ is useful for retarded/advanced Green's functions.