Consider the following integral for $x > 0$: $$ F(x) \ = \ \mathscr{P} \int_{0}^{\infty} \frac{d\omega}{e^{\omega} - 1} \frac{\omega}{\omega^{2} - x^{2}} $$
The $\mathscr{P}$ is there because there is a singularity at $\omega = x$, so I suppose we can write the above function more carefully as: $$ F(x) \ = \ \lim_{\epsilon \to 0^{+}} \left\{ \int_{0}^{x - \epsilon} \frac{d\omega}{e^{\omega} - 1} \frac{\omega}{\omega^{2} - x^{2}} \ + \ \int_{x + \epsilon}^{\infty} \frac{d\omega}{e^{\omega} - 1} \frac{\omega}{\omega^{2} - x^{2}} \right\} $$
It looks like $F(0) = \int_{0}^{\infty} \frac{d\omega}{e^{\omega} - 1} \frac{1}{\omega}$ which is divergent (the integrand looks like $\frac{1}{\omega^2}$ near the $\omega =0$ in this case).
I am interested in obtaining asymptotics for $F(x)$ as we approach $x \to 0$. I am unsure how to do this due to the presence of the principal value $\mathscr{P}$. Is there a method for dealing of asymptotics of such integrals?
I've figured it out. I've even got an analytic result for this integral! $$ F(x) := \mathscr{P} \int_{0}^{\infty} \frac{d\omega}{e^{\omega} - 1} \frac{\omega}{\omega^2 - x^2} = \frac{1}{2} \ln\left( \frac{x}{2\pi} \right) - \frac{1}{2} \mathrm{Re}\left[ \psi^{(0)}\left( \frac{ix}{2\pi} \right) \right] $$
And as we take $x \to 0^{+}$ the above looks like: $$ F(x) = \frac{1}{2} \ln\left( \frac{x}{2\pi} \right) \ + \ \frac{\gamma}{2} \ + \ \mathscr{O}(x^2) $$
I've done some numerical tests in Mathematica comparing the integral using the principal value method, the full formula, and then the asymptotic form and everything works.
Proof: Note the following formula (see this post) which is derived by differentiating Binet's Second Formula: $$ G(y) : = \int_{0}^{\infty} \frac{d\omega}{e^{\omega} - 1} \frac{\omega}{\omega^2+y^2} \ = \ \frac{1}{2} \ln\left( \frac{y}{2\pi} \right) - \frac{\pi}{2 y} - \frac{1}{2} \psi^{(0)}\left( \frac{y}{2\pi} \right) $$
We take an analytical continuation by considering $\lim\limits_{\epsilon \to 0^{+}}G(\epsilon+i x)$ for $x>0$. We find: $$ \lim\limits_{\epsilon \to 0^{+}}G(\epsilon+i x) = \lim_{\epsilon \to 0^{+}} \int_{0}^{\infty} \frac{d\omega}{e^{\omega} - 1} \frac{\omega}{\omega^2+(\epsilon + i x)^2}\to\lim_{\epsilon \to 0^{+}} \int_{0}^{\infty} \frac{d\omega}{e^{\omega} - 1} \frac{\omega}{\omega^2 - x^2 + i \epsilon} $$
Then we can use the identity $\lim\limits_{\epsilon \to 0^{+}}\frac{1}{z+i \epsilon} = \mathscr{P}\frac{1}{z} - i \pi \delta(z)$. Taking $\epsilon \to 0^{+}$ on both sides gives: $$ G(i x) = \mathscr{P}\int_{0}^{\infty} \frac{d\omega}{e^{\omega} - 1} \frac{\omega}{\omega^2 - x^2} \ - \ i \pi \int_{0}^{\infty} d\omega \frac{\omega}{e^{\omega} - 1} \delta\big(\omega^2 - x^2 \big) $$
Where we can evaluate the imaginary part easily, and the real part is the function we're after: $$ G(i x) = F(x) \ - \ i \frac{\pi}{2} \frac{1}{e^{x} - 1} $$
This tells us that: $$ F(x) = \mathrm{Re}\left[ G(i x) \right] = \mathrm{Re}\left[ \frac{1}{2} \ln\left( \frac{ix}{2\pi} \right) - \frac{\pi}{2 ix} - \frac{1}{2} \psi^{(0)}\left( \frac{ix}{2\pi} \right) \right] $$
Which for $x>0$ simplifies to the stated result.
$\ $
BONUS: Using the above, we find that for $\alpha >0$: $$ \mathrm{Im}\left[ \psi^{(0)}\left( i \alpha \right) \right] \ = \ \frac{1}{2 \alpha} + \frac{\pi}{2} \coth\left( \pi \alpha \right) $$ I don't know if this is a known formula, but I certainly couldn't find it anywhere (UPDATE: It is known, I found it here on DLMF!). It works numerically, and seems to also work for $\alpha <0$. This allows us to write $F(x)$ in the form: $$ F(x) = \frac{1}{2} \ln\left( \frac{x}{2\pi} \right) + \frac{i \pi}{2 x } + \frac{i \pi}{4} \coth\left( \frac{x}{2} \right) - \frac{1}{2} \psi^{(0)}\left( \frac{ix}{2\pi} \right) $$ which is purely real, despite the appearance of the $i$'s.