Which limit for Matsubara frequency summation? In the context of a simple toy problem for Feynman path integrals, I consider a two-site Hubbard model for spinless fermions. I expand the path integral to first order in the interaction $V$, which means that I have to compute an "average" of the type
$$\left\langle \bar \psi_1(\omega_n) \bar \psi_2(\omega_m) \psi_2(\omega_p) \psi_1(\omega_q) \right\rangle$$
with respect to the non-interacting action.
Since the Grassmann numbers for sites $1$ and $2$ are different numbers, this average immediately factors, and I'm left with two-variable averages, which give me just the non-interacting Matsubara Green's functions.
But now I have to calculate a Matsubara sum of the type
$$\frac{1}{\beta} \sum_n \frac{1}{i\omega_n - \epsilon}$$
where $\epsilon$ is just a constant. 
From searching in the literature, I know that contour integration is the way to go, but unfortunately this sum is ambiguous and I have to introduce a factor $e^{-i\omega_n \tau}$ with $\tau \rightarrow 0^+$ or $\tau \rightarrow 0^-$, and the result of the contour integration will be quite different depending on which limit I take. 
My question now is, how do I know what limit to take? And is there an easier way of carrying out the sum over frequency?
 A: Yes, for this particular Matsubara sum you mentioned, taking different limits will lead to two results differed by 1. This is because the summation in consideration does not converge, which can be seen from the following integral (the continuous limit of the sum) considering the ultraviolet divergence (the large $\omega$ behavior) 
$$\frac{1}{2\pi}\int_{-\infty}^{\infty}\mathrm{d}\omega\frac{1}{i\omega-\epsilon}\sim\frac{1}{2\pi i}\int_{-\infty}^{\infty}\mathrm{d}\omega\frac{1}{\omega}\rightarrow\infty.$$
Asking for the result of an essentially divergent sum will not even lead to a definite answer. Introducing the factor $e^{i\omega_n\tau}$ is to control the convergence of the sum, but the result will then depend on how one chooses to control the convergence, i.e. $\tau\rightarrow0^+$ (controlling the left-half complex plane) or $0^-$ (controlling the right-half complex plane).
The two results will be different by a shift of 1. Even this difference "1" has a physical meaning, reflecting the ambiguity in defining Fermion particle number. Because the summation
$$\frac{1}{\beta}\sum_n\frac{1}{i\omega_n-\epsilon}=G(\tau=0)=-\mathcal{T}_{\tau}\langle\psi(\tau=0)\bar{\psi}\rangle$$
corresponds to the Green's function at imaginary time 0, which physically means to count the number of Fermions on the energy level $\epsilon$. Here $\mathcal{T}_{\tau}$ stands for the time ordering operator. However $\tau=0$ can be either understood as $\tau\rightarrow0^+$ or $\tau\rightarrow0^-$ (corresponding to the two ways to control the convergence), which becomes very critical here because the time ordering operator will then order the two cases differently:
$$-\mathcal{T}_{\tau}\langle\psi(\tau=0^+)\bar{\psi}\rangle = -\langle\psi\bar{\psi}\rangle,$$
$$-\mathcal{T}_{\tau}\langle\psi(\tau=0^-)\bar{\psi}\rangle = \langle\bar{\psi}\psi\rangle.$$
According to the anticommutation relation of the Fermion operators, we have $\bar{\psi}\psi=1-\psi\bar{\psi}$, that is why it is expected that the two results will differ by a constant 1. In fact $\bar{\psi}\psi$ counts the number of particles, while $\psi\bar{\psi}$ counts the number of holes. It becomes a problem whether to define the Fermion number (with respect to the vacuum state) as the particle number or as the negative of the hole number. This depends on our definition of the vacuum state: whether we treat the vacuum as no particle state or as a state filled with particles (like the concept of Fermi sea). We know both choices are acceptable, just as we can either define electron as a particle in the empty space or as a hole in the anti-electron Fermi sea. Because of this ambiguity, one needs to specify the choice when the $\tau\rightarrow0$ limit is taken. In conclusion, both ways of taking the limit are legitimate, and the difference in the results is just a matter of the definition of Fermion vacuum state.
Similar discrepancy also exists in the Bosonic case, which is again related to the definition of Boson vacuum state. In general, such discrepancy roots in the divergence of the Matsubara summation in consideration. If the summation itself converges to a definite result, then the result will be unique no matter which limit is chosen to take. For example, if one tries to calculate
$$\frac{1}{\beta}\sum_n\frac{1}{(i\omega_n-\epsilon)^2},$$
the result will always be $\beta/4 \mathrm{sech}^2 (\beta\epsilon/2)$ (for Fermionic case) no matter $\tau\rightarrow0^+$ or $0^-$.
A: Here I will list the recipe to compute the Matsubara summation, namely how to sum the following series 
$$\dfrac{1}{\beta}\sum_{n}g(i\omega_n),$$
where $\omega_n$ is Matsubara frequencies and n will be integers.


*

*Analytic continuation:
$$g(i\omega_n) \rightarrow g(z)$$

*Find the poles $z_m$ of $g(z)$ and compute 
$$\sum_{m}f(z_m)Res(g,z_m)$$


where $f(z)$ can be Fermi or Bose distribution function depending on the question of interest and note that the $z_m$ will not stay on the imaginary axis.
For example:
$$\dfrac{1}{\beta}\sum_n\dfrac{1}{i\omega_n + E}=f(-E)Res(g(z),-E)=\dfrac{1}{e^{-\beta E}-1}$$
where $$\omega_n=\dfrac{2n\pi}{\beta}$$
are bosonic frequencies.
