Path integral calculations $e^{i\omega 0^+}$ When computing correlation functions using the path integral formulation, I often need to compute integrals such as
$$
\int_{-\infty}^\infty \frac{d\omega}{2\pi} \frac{1}{i\omega -\epsilon}
$$
Obviously, this integral doesn't converge, so I often see textbooks use the convention
$$
\int_{-\infty}^\infty \frac{d\omega}{2\pi} \frac{e^{i\omega 0^+}}{i\omega -\epsilon}
$$
I understand choosing $0^+$ or $0^-$ will affect how we choose the contour integral in the upper half or lower half plane. My question is before that, i.e., why do we choose $0^+$ or $0^-$ given the first integral?
So let me provide a concrete example, the toy Hubbard model. Consider the fermionic path integral with free action
$$
S_0 = \int_{-\infty}^\infty \frac{d\omega}{2\pi} \sum_{i=1,2}\bar{\psi}_i(\omega)(i\omega -\epsilon)\psi_i(\omega)
$$
And quartic-perturbation (schematically)
$$
U\int \bar{\psi}_1\bar{\psi}_2\psi_1\psi_2
$$
Then
$$
\langle \bar{\psi}_1(\omega)\psi_1(\omega)\rangle = \frac{\langle \bar{\psi}_1(\omega) \psi_1(\omega) e^{U\int \bar{\psi}_1\bar{\psi}_2\psi_1\psi_2} \rangle_0}{\langle e^{U\int \bar{\psi}_1\bar{\psi}_2\psi_1\psi_2} \rangle_0}
$$
where the subscript $0$ denotes average with respect to the free action $S_0$. Then to the first order of $U$, we have something like
$$
\langle \bar{\psi}_1(\omega)\psi_1(\omega)\rangle = \frac{2\pi \delta(0)}{i\omega -\epsilon} +\frac{2\pi \delta(0)}{(i\omega -\epsilon)^2} \int \frac{d\omega}{2\pi} \frac{1}{i\omega -\epsilon}
$$
I see authors writing the integral in the second term as
$$
\int_{-\infty}^\infty \frac{d\omega}{2\pi} \frac{e^{i\omega 0^+}}{i\omega -\epsilon}
$$
But I'm not sure about the justification.
 A: Why the integral is ‘easy’
So it is important to understand that the non-convergence in this case is “not so bad” in the sense that if you just do the integral you would first multiply both top and bottom by complex conjugate  as $$I=\int_{-\infty}^\infty\frac{d\omega}{2\pi}~\frac{-i\omega - \epsilon}{\omega^2 + \epsilon^2},$$ and then you can argue that the term $\omega/(\omega^2 + \epsilon^2)$ is odd in $\omega$ and should be regarded as having integral zero, so that all that is left is $$
\begin{align}
I&=-\int_{-\infty}^\infty\frac{d\omega}{2\pi}~\frac{\epsilon}{\omega^2 + \epsilon^2}=-\int_{-\infty}^\infty\frac{\epsilon~du}{2\pi}~\frac{\epsilon}{\epsilon^2u^2 + \epsilon^2}\\
&=-\frac{1}{2\pi}\int_{-\infty}^{\infty}\frac{du}{u^2 + 1}=-\frac{1}{2\pi}\left(\tan^{-1}(\infty) - \tan^{-1}(-\infty)\right)\\
&=-\frac{\pi}{2\pi} =-\frac12.
\end{align}
$$
The point is that this is in some sense a very tractable problem. Formally to denote that result you have to preface the integral with the letters $\operatorname{PV}$ denoting Cauchy’s principal value, and then we do indeed have the above result; the idea that “well, $\omega/(\omega^2 + \epsilon^2)$ is odd and therefore zero” hushes over creative ways of expanding one of these integral bounds exponentially faster than the other one and thereby extracting some nonzero terms from an integral that looks like $\int dx/x.$
Brief residue theorem refresher
But what we want to do is to use some tools of complex analysis so that we can be lazy.
One tool of complex analysis is called the residue theorem, and it is a special case of a more broad thing called contour deformation. The idea is that if you have an integration contour which is a circle about some point, you can then arbitrarily deform that contour through any regions of the complex plane on which your function is analytic (generally meaning both nonsingular and non-branch-cut, though you can also deform into branch cuts if you are careful to go onto a new branch sheet). In particular, you can shrink it into a tiny circle around that point. Now around that point $z_0$ the function is usually described by some Laurent series $$f(z) = \sum_{k=-\infty}^\infty c_k~(z-z_0)^k,$$ but letting $z = z_0 + r e^{i\theta}$ for $0< \theta<2\pi$ for some smallish nonzero $r$, let's call that the contour $C_r$; this integral becomes $$\oint_{C_r} dz~f(z) =\sum_{k=-\infty}^\infty \int_0^{2\pi}i~r~e^{i\theta}d\theta~c_k r^k e^{ik\theta}.$$Irrespective of the exact $r$, this is sinusoidal in $\theta$ and thus integrates to zero unless $k=-1$: and that term has no $r$- or $\theta$-dependence and simply gives $$\oint_{C_r} dz~f(z) = 2\pi i~c_{-1}.$$This term $c_{-1}$ is known as the residue of an isolated singularity; if you have multiple isolated singularities the argument just gives $2\pi i$ times a sum over their individual residues. Note that due to our choices about $\theta$ this is a counterclockwise contour; for a clockwise contour you would get $-2\pi i.$
But we do not have a closed contour
The residue theorem is extremely easy to apply here: the thing has a pole at $\omega_0=-i\epsilon$ and looks like $(2\pi i)^{-1} (\omega + i \epsilon)^{-1}$ and so $1/(2\pi i)$ is the residue, for any closed contour that wraps around $-i\epsilon.$ (As you have more experience you literally do not even do any of this reasoning, you just look at it and say, “oh, $1/(2\pi i \omega)$, cancel the $\omega.$”) 
Since the $\omega$ axis passes over this singularity we are looking at a clockwise contour. So we will get $(-2\pi i)/(2\pi i) = -1,$ easy peasy, right?
Well, the first clue is that we are not getting the same answer both ways. And it’s not due to the subtleties of the Cauchy principal value—that would have contributed something imaginary, not real. The problem is ultimately traced back to the fact that we proved the residue theorem for closed contours, this is an open contour.
So we introduce a trick, we add and then subtract a semicircle of the contour off at infinity, $R e^{i\theta}$ for very large $R$ and for $\theta$ ranging from $0$ to $-\pi.$ Then the part that we have added, combines with the integral over the real axis to form a closed contour. The part that we have subtracted is all we need to really look at. Alternatively you can think that we deform the original contour from the real line to the semicircle $R e^{i\theta}$ for large $R$ and $\theta$ ranging from $-\pi$ to $0$, but it leaves behind a little clockwise contour about the point at $-i\epsilon.$ Both are the same thing, just depends on whether you like adding-and-subtracting contours or pushing-them-around.
The remaining contour is then $$\frac1{2\pi i}\int_{-\pi}^0 iRe^{i\theta}d\theta~\frac{1}{Re^{i\theta} + i \epsilon}.$$For large $R$ we can then expand this as $1/(1 + (\epsilon/R)e^{-i\theta})$ and do a geometric series but essentially we will find that the limit $R\to\infty$ is the same as the limit $\epsilon\to 0$ and either way this evaluates directly as $i\pi/(2\pi i) = \frac12.$ So that’s the rest of it, $-1 + \frac12 = -\frac12.$ It is a residue that we can just see by inspection plus a contour that is not terribly hard to do.
But I am even lazier than that
But you might be wondering if you can be even lazier. Surely there is a general principle such that I don’t need to reason about contours at infinity, right? And indeed there is, and that is the source of this term you are looking at.
We multiply $f(z)$ (which does not have a zero contour off at infinity), by $e^{q z}$ for some nonzero $q$ that will make the value on the infinity-contour zero. 
What does $q$ have to be? Well, on this contour, $z$ will be $x + i y$ for two extremes: either $x$ is going to be both very very positive or very very negative (and $y$ will be near zero) or else $y$ will be very very negative (and $x$ may be near zero). The big-$x$ side of things we want to solve by making this term $e^{qx}$ oscillate wildly and thus not integrate to anything useful: so $q$ is pure-imaginary, say $q=is.$ The $y$ side of things we want to solve by making this term $e^{qz} = e^{is~iy} = e^{-sy}$ just become an $e^{-\infty} = 0$ expression. This way, the whole contour is zero. For very very negative $y$, this requires that $s$ also be negative. So we multiply by $e^{-i0^+\omega}$ because it gives us a space out at $\omega \approx -i\infty$ where contours are “free,” they are all exponentially damped out to zero.
We can also multiply by $e^{i0^+\omega}$ in which case contours at $\omega \to +i\infty$ are “free.” This is one of many implications of Jordan’s lemma. The mnemonic is really simple, $e^{+i\eta z}$ to close in the positive imaginary side, $e^{-i\eta z}$ to close in the negative imaginary side.
You may have noticed that you have two options: one is to close up top, in which case the contour contains no poles and the answer is zero, or to close below, in which case the contour contains one pole at $-i\epsilon$ and the answer is $-1$ as discussed above. So actually we are seeing a discontinuity in this function as a function of $s$, and that means we cannot meaningfully take the limit as $s \to 0.$ This comes from the $e^{i\eta \omega}$ terms having an odd imaginary component which is interacting with the existing odd imaginary component in $f(\omega)$ to produce an even real component. 
The solution is to actually use $\cos(0^+ \omega) \to 1$, which is purely real, which effectively means that you close both top and bottom and take the average of the two. So I see $0$ and $-1$ and I know that the average is $-1/2$, which should be the integral. Very lazy, indeed!
A: I revisited this problem after some time, and even though I think @CR Drost provided a very informative answer, I don't quite think it's "exactly" the solution to my problem.
Indeed, after doing the calculation rigorously and following each step carefully, I think the answer is that the fermionic action (at zero temperature) is actually supposed to be
$$
S_0 = \int_{-\infty}^\infty \frac{d\omega}{2\pi} \bar{\psi}(\omega)(i\omega+0^+ -(\epsilon -\mu) e^{i\omega 0^+})\psi(\omega)
$$
where the $0^+$ and $e^{i\omega 0^+}$ are factors that occurs when applying the Lie-trotter formula, i.e., $e^{-\beta H} = (e^{-\beta H/M})^M$ and we take $M \to \infty$.
In this case, if I were to say calculate the average particle number at zero temperature, I would have
$$
\langle N \rangle =\frac{1}{\beta Z} \frac{\partial}{\partial \mu}Z, \quad \beta \to \infty
$$
And by taking the derivative, the extra factor $e^{i\omega 0^+}$ occurs in the integral.
One small detail, I haven't checked whether we can rigorously interchange limits $M\ \to \infty$ with partial derivatives, integrals and $\beta \to \infty$ (maybe I'll try in the future and come back to this question). However, unless someone else points out a significant flaw in the logic, I believe that this is actually where the $e^{i\omega 0^+}$ comes from, and (at least most) other answers are just "possible justifications" of understanding an ill-defined integral.
