thermal field theory - Matsubara frequencies as poles of distribution function

Is there any deeper meaning to why the bosonic/fermionic Matsubara frequencies appear as poles of their corresponding distribution functions (with an additional $i$)?

For example in the bosonic case we have: $$\omega_n=\frac{2n\pi}{\beta}$$ which are the poles of $$n_B=\frac{1}{e^{\beta \hbar \omega}-1}$$ if we say that $\omega = i\omega_n$

I know this is sometimes exploited when evaluating sums of functions of $\omega_n$.

See Matsubara frequency.

Answer

The point is that to achieve the sum over Matsubara frequencies $$\sum_{n} g(i\omega_n)$$ we can use a contour integral $$\oint_C g(z) f(z)$$ with the contour described in fig 1 here, so long as we choose an $f(z)$ with simple poles exactly at the Matsubara frequencies $\omega_n$. This determines $f(z)$ to be proportional to the Bose-Einstein distribution. In fact, if we think about the sum as computed in a correlation function, we could have derived the integral above instead by considering computing the expectation value in the thermal density matrix $\int dE n(\beta,E) |E\rangle\langle E|$ determined by the Bose-Einstein distribution $n(\beta,E)$. I think this is an equally good starting point of the logic, then deriving the Matsubara frequencies as the poles of the distribution function.

I'd say the whole point of the Matsubara frequencies is that they are the poles of the distribution function.