Evaluating the following integral: $\int_{0}^{\infty }\frac{x\cos(ax)}{e^{bx}-1}\ dx$

We first start by converting our integral into a sum. We have that (by the geometric series): $$\frac{1}{e^{bx}-1}=\sum_{n=1}^{\infty} e^{-bxn}$$ Additionally, we have that (as can be found by simply finding an antiderivative): $$\int_0^{\infty} x\cos(ax)e^{-bxn}~dx=\frac{-a^2+b^2 n^2}{(a^2+b^2 n^2)^2}$$ Therefore, the integral in question is equal to the following sum: $$\int_{0}^{\infty }\frac{x\cos(ax)}{e^{bx}-1}\ dx=\sum_{n=1}^{\infty} \frac{-a^2+b^2n^2}{(a^2+b^2 n^2)^2}=-a^2 S_1+b^2 S_2 \tag{1}$$ Therefore, it suffices to compute the following sums: $$S_1:=\sum_{n=1}^{\infty} \frac{1}{(a^2+b^2 n^2)^2},\quad S_2:=\sum_{n=1}^{\infty} \frac{n^2}{(a^2+b^2 n^2)^2}$$ To do this, we start from the well-known result that (which can be derived using real analysis, as shown by several answers): $$\sum_{n=0}^\infty\frac{1}{c^2+n^2}=\frac{1+c\pi\coth (c\pi)}{2c^2}$$ From this, by using the substitution $c:=a/b$ one can easily derive the more general result that (note the difference in the lower limit of the sum): $$\sum_{n=1}^{\infty} \frac{1}{a^2+b^2 n^2}=\frac{-b+a\pi \coth(a\pi/b)}{2a^2 b} \tag{2}$$ Now, we can compute closed forms for the sums $S_1$ and $S_2$ by differentiating $(2)$ with respect to $a$ and $b$ respectively. We hence have that: $$S_1=\frac{\pi^2 a^2 \operatorname{csch}^2(a\pi /b)+\pi a b \coth(a\pi /b)-2b^2}{4a^4 b^2}$$ $$S_2=\frac{-\pi^2 a^2 \operatorname{csch}^2(a\pi /b)+\pi a b \coth(a\pi /b)}{4a^2 b^4}$$ Therefore, by $(1)$, the integral is given by: $$\bbox[5px,border:2px solid #C0A000]{\int_{0}^{\infty }\frac{x\cos(ax)}{e^{bx}-1}\ dx=\frac{1}{2} \left(\frac{1}{a^2}-\frac{\pi ^2 \operatorname{csch}^2(a\pi/b)}{b^2}\right)}$$

Let $$ I(a)=\int_{0}^{\infty}\frac{\sin(ax)}{e^{bx}-1}\ dx. $$ It is easy to see $$ I'(a)=\int_{0}^{\infty}\frac{x\cos(ax)}{e^{bx}-1}\ dx. $$ Since \begin{eqnarray} I(a)&=&\int_{0}^{\infty}\frac{e^{-bx}\sin(ax)}{1-e^{-bx}}\ dx=\int_{0}^{\infty}\sum_{n=0}^\infty e^{-b(n+1)x}\sin(ax)\ dx\\ &=&\sum_{n=0}^\infty\int_{0}^{\infty} e^{-b(n+1)x}\sin(ax)\ dx=\sum_{n=0}^\infty\frac{a}{a^2+b^2(n+1)^2}. \end{eqnarray} Using the result from @projectilemotion, one has \begin{eqnarray} I(a)&=&\sum_{n=0}^\infty\frac{a}{a^2+b^2(n+1)^2}=\frac{-b+a\pi \coth(a\pi/b)}{2a b}=\frac12\left(-\frac{1}{a}+\frac\pi b\coth(\frac{a\pi}{b})\right) \end{eqnarray} and hence $$ I'(a)=\frac12\left(\frac1{a^2}-\frac{\pi^2}{b^2}\text{csch}^2(\frac{a\pi}{b})\right).$$