On the Wikipedia page for Matsubara frequencies, the following formula is presented, $$ \sum_{i\omega_n} \frac{(i\omega_n)^2}{(i\omega_n)^2 - \xi^2} = -\frac{\xi}{2}\Big(1 - 2 N_{\text{FD}}(\xi)\Big), $$ where $\omega_n = (2n+1)\pi/\beta$ are fermionic Matsubara frequencies and $N_{\text{FD}}(x):= (e^{\beta x}+1)^{-1}$ is the Fermi-Dirac distribution function.
I am familiar with the residue theorem and I can derive `easier' results involving Matsubara summations. However, in this case I don't see how to derive the result. I guess it is necessary to include the usual converging factor $e^{i\omega_n\eta}$ with $\eta\rightarrow0$ at the end of the computation, otherwise the sum would not even converge. Any help or reference would be appreciated.
The reason I am asking this, is because I actually want to evaluate a Matsubara sum of the form $$ \sum_{i\omega_n} \frac{(i\omega_n)^2}{(i\omega_n - \xi_1)(i\omega_n - \xi_ 2)}. $$