Matsubara Sums and Multiple Poles In Mahan's book, equation (4.127), he claims that
\begin{align}
&\frac{1}{\beta}\sum_{ik_n} \frac{1}{ik_n-\xi_1}\frac{1}{ik_n-\xi_2}\frac{1}{ik_n-\xi_3} \\
=& \frac{n_F(\xi_1)}{(\xi_1-\xi_2)(\xi_1-\xi_3)} + \frac{n_F(\xi_2)}{(\xi_2-\xi_1)(\xi_2-\xi_3)} + \frac{n_F(\xi_3)}{(\xi_3-\xi_1)(\xi_3-\xi_2)} - \pi^2 n_F(\xi_1)\delta(\xi_1-\xi_2)\delta(\xi_1-\xi_3)
\end{align}
where $ik_n$ are fermionic (half-integer) frequencies and $n_F$ is the Fermi distribution function. Because one eventually integrates over $\xi_1,\xi_2,\xi_3$, one needs to be careful about the choice of contour when computing the sum via residues; namely, one should avoid crossing the real axis.
The contour one should take is two infinite-radius semicircles, one $i\epsilon$ above and one $i\epsilon$ below the real axis, both with counter-clockwise orientation. Hence, the Matsubara sum evaluates to the difference between the integrals $\int_{-\infty+i\epsilon}^{\infty+i\epsilon} - \int_{-\infty-i\epsilon}^{\infty-i\epsilon}$.
The first three terms of (4.127) are simple to obtain, and indeed what you would expect if one just fixes $\xi_1,\xi_2,$ and $\xi_3$ to be some constant value. However, I am having trouble demonstrating that the last term is actually correct. 
Naively, one would want to use the Sokhotski–Plemelj theorem 
$$\frac{1}{x+i\epsilon} = P\frac{1}{x} - i\pi\delta(x)$$
term-by-term, but such a relation only holds when the function against which you are integrating is analytic. The delta-functions in the 4th term of (4.127) clearly put all the poles at the same spot, disallowing this naive replacement term-by-term.
How should one go about demonstrating that this 4th term is really the correct $\epsilon\rightarrow 0$ limit? Or is Mahan even correct?
Edit: Also, this formula implies that interchanging the Matsubara sum and the momentum integrals changes the answer; but I would've expected this not to be the case. It is curious to me that you are free to interchange with two factors, but not with three or more. Why is this the case, and which one should I consider to be "correct" order of operations and why? 
 A: The end result looks  one form of the  Poincare-Bertrand theorem.This is  a much deeper result that  the usual Sokhotski–Plemelj theorem. There is a discussion in my lecture notes: https://courses.physics.illinois.edu/phys509/sp2019/bmaster.pdf
on 381. When I wrote these notes  I did know that Hardy's theorem is now known by the name of  Poincare-Bertrand despite the fact G H  Hardy had it first. It's general principle that a theorem  is seldom  known by the name of the first person to discover it  but instead by the name of the last person to discover it.
I've not worked out how Matsubara sums work here, but I can explain some of the wierd math of Poincare-Bertrand. 
One way to understand P-B is via the Hilbert transform and the convolution theorem. 
The Hilbert transform of a test function $f(x)$ is 
$$
{\mathcal H}f(x) = \frac{P}{\pi } \int_{-\infty}^\infty \frac{f(y)}{x-y} dy,
$$
(where $P$ denotes a principal part integral) and its  Fourier transform is 
$$
\widetilde{{\mathcal H}f}(\omega) = i {\rm sgn}(\omega)\tilde f,
$$
where $\tilde f$ is the Fourier transform of $f$.
Now if $\omega_1+\omega_2+\omega_3=0$ we have
$$
{\rm sgn}(\omega_1){\rm sgn}(\omega_2)+ {\rm sgn}(\omega_2){\rm sgn}(\omega_3)+{\rm sgn}(\omega_3){\rm sgn}(\omega_1)=-1,
$$
(Two of the frequencies must have one sign and the remaining one the other sign, so two of the terms are negative and one positive.)
Now one can use this expression in the  convolution theorem to get the result that  $$
\int \frac{u_2(z)}{x-z} \left(\int \frac{u_1(y)}{z-y}\,dy\right) dz 
+\int \frac{u_1(y)}{x-y} \left(\int \frac{u_2(z)}{y-z}\,dz\right) dy
+ \int \frac{u_1(y)}{y-x}\,dy \int \frac{u_2(z)}{x-z}\, dz = -\pi^2 u_1(x)u_2(x). 
$$
Principal part integrals are implied in all these terms.
We can rewrite   this (without doing anything more dangerous than pulling constants under the integration signs)  as 
$$
\int\left(\int\frac{u_1(y)u_2(z)}{(x-z)(z-y)}\, dy\right) dz+ \int\left(\int \frac{u_1(y)u_2(z)}{(x-y)(y-z)}dz\right)dy
+ \int \left(\int\frac{u_1(y)u_2(z)}{(y-x)(x-z)} dy\right)dz 
 = -\pi^2 u_1(x)u_2(x).$$
$P$ integrals are still implied whenever we try to divide by zero.
Now we can (again innocently) combine  the integrands of the first and third terms using 
$$
\frac{1}{(y-z)(z-x)} +\frac{1}{(z-x)(x-y)}= - \frac1{(x-y)(y-z)}
$$
to get 
$$
-\int_{-\infty}^\infty \left(\int_{-\infty}^\infty\frac{u_1(y)u_2(z)}{(x-y)(y-z))}\, dy\right) dz+ \int_{-\infty}^\infty \left(\int_{-\infty}^{\infty} \frac{u_1(y)u_2(z)}{(x-y)(y-z)}dz\right)dy
  = -\pi^2 u_1(x)u_2(x),
$$
which shows that it is not allowed to interchange the order of integrations.
I'm going to think about the Matsubara sums.
