Zero-Temperature Limit of Matsubara Sums Consider the fermionic Matsubara sum for
$$\frac{1}{\beta}\sum_n G^2(k,\omega_n) = \frac{1}{\beta}\sum_n \frac{1}{(i\omega_n - \epsilon_k)^2}$$
where $G(k,\omega_n)$ is the free fermion Green's function. To evaluate this sum, one can do the standard procedure by introducing the fermionic distribution function $n_F$, providing poles at $(2n+1)\pi/\beta$ with residue $-1/\beta$. This gives the result
$$\frac{1}{\beta}\sum_n G^2(k,\omega_n) = n'_F(\epsilon_k)$$
In the $T\rightarrow 0$ limit, this evaluates to $\delta(\epsilon_k)$.
However, now let us first take the $T\rightarrow 0$ limit first. We can then change the sum into an integral, giving us
$$\int_{-\infty}^\infty \frac{d\omega}{2\pi} \frac{1}{(i\omega - \epsilon_k)^2}$$
Evaluating this integral by contour integration, however, gives $0$ due to the fact that it is a double pole.
Is there something wrong with my manipulations here, or is there something inherently tricky about taking the $T\rightarrow 0$ limit?
 A: One needs to be careful when taking the $T\rightarrow 0$ limit. Consider $\epsilon_k = 0$. In this case, the integral form diverges. Therefore, the statement that this integral is always $0$ is clearly incorrect. 
One takes a clue of the proper procedure by realizing that $\delta(\epsilon_k)$ is a distribution, and hence must be "properly" understood within an integral. Therefore, we really want to be considering the quantity:
$$\lim_{\beta\rightarrow\infty} \int d\epsilon_k f(\epsilon_k) \frac{1}{\beta}\sum_{n}\frac{1}{(i\omega_n - \epsilon_k)^2}$$
where we will take $f(\epsilon_k)$ to be as nice as possible (smooth, vanishes sufficiently quickly at $\infty$, etc). Now, we need to be painstakingly careful about the order of operations. 
First, let's evaluate this quantity in the order shown.
\begin{align}
&\lim_{\beta\rightarrow\infty} \int d\epsilon_k f(\epsilon_k) \frac{1}{\beta}\sum_{n}\frac{1}{(i\omega_n - \epsilon_k)^2} \\
=& \lim_{\beta\rightarrow\infty} \int d\epsilon_k f(\epsilon_k) n_F'(\epsilon_k) \\
=& f(0)
\end{align}
Instead, we could have "moved the limit inside the integral" by treating $n_F(\epsilon_k)$ as a distribution, so that $\lim_{\beta\rightarrow \infty} n_F'(\epsilon_k) = \delta(\epsilon_k)$ in the sense of distributions, and arrived at the same result.
Now, how should we understand the limiting procedure converting the sum into an integral in the sense of distributions? Using our "sufficiently nice" properties of $f(\epsilon_k)$, we can interchange the sum and the integral:
\begin{align}
&\lim_{\beta\rightarrow\infty} \frac{1}{\beta}\sum_{n} \int d\epsilon_k f(\epsilon_k) \frac{1}{(i\omega_n - \epsilon_k)^2} \\
=& \int \frac{d\omega}{2\pi} \int d\epsilon_k f(\epsilon_k) \frac{1}{(i\omega - \epsilon_k)^2}
\end{align}
where in the second line we used the definition of the (Riemann) integral. We wish to interchange these two integrals so that we evaluate $\int d\omega$ first; classically this is not allowed since the integral is not absolutely convergent. However, let us do the following trick:
\begin{align}
&\int \frac{d\omega}{2\pi} \int d\epsilon_k f(\epsilon_k) \frac{1}{(i\omega - \epsilon_k)^2} \\
=& \lim_{\epsilon \rightarrow 0} \int_{\mathbb{R}-(-\epsilon,\epsilon)} \frac{d\omega}{2\pi} \int d\epsilon_k f(\epsilon_k) \frac{1}{(i\omega - \epsilon_k)^2} \\
=& \lim_{\epsilon \rightarrow 0} \int d\epsilon_k f(\epsilon_k) \int_{\mathbb{R}-(-\epsilon,\epsilon)} \frac{d\omega}{2\pi} \frac{1}{(i\omega - \epsilon_k)^2} \\
=& \lim_{\epsilon \rightarrow 0} \int d\epsilon_k f(\epsilon_k) \frac{1}{\pi}\frac{\epsilon}{\epsilon^2 + \epsilon_k^2} \\
=& f(0)
\end{align}
In the second line, we remove the problematic region so that we can interchange the integrals; we take the principal value of the integral. Note that in the second to last line, a nascent delta-function $\frac{1}{\pi}\frac{\epsilon}{\epsilon^2 + \epsilon_k^2}$ appears!
In other words, in order to take the $\beta\rightarrow \infty$ limit first to convert the sum into an integral, we need to be calculating the principal value of the integral, with the $\epsilon \rightarrow 0$ limit being taken in the sense of distributions (ie after the $\int d\epsilon_k$ integral). Concretely, this means
\begin{align}
&\frac{1}{\beta}\sum_n\frac{1}{(i\omega_n -\epsilon_k)^2}\\
=&\int_{\mathbb{R}-(-\epsilon,\epsilon)} \frac{d\omega}{2\pi} \frac{1}{(i\omega - \epsilon_k)^2} \\
=& \int_{\mathbb{R}-(-\epsilon,\epsilon)} \frac{d\omega}{2\pi} \frac{1}{(i\omega - \epsilon_k)^2} \\
=& \frac{1}{\pi}\frac{\epsilon}{\epsilon^2 + \epsilon_k^2} \\
=& \delta(\epsilon_k)
\end{align}
