Due to the simple poles of $\tan(ax)$ and the slow decay of $\frac{x}{x^2+b^2}$, the integral should be interpreted in a very careful way. So here is one way to make sense of the integral:
Let $(R_k)_{k\geq1}$ be an increasing sequence in $(0, \infty)$ satisfying the following conditions:
- $R_k \to +\infty$ as $k \to \infty$, and
- there exists $\delta > 0$ such that $\lvert R_k - \frac{\pi}{2a}(2n-1) \rvert \geq \delta $ for any $k, n \geq 1$. That is, each $R_k$ is at least at a distance $\delta$ from any of the poles of $\tan(ax)$. (For example, $R_k = \frac{\pi k}{a}$ works.)
Then
$$ \mathcal{I} := \lim_{k\to\infty} \operatorname{PV}\!\!\int_{0}^{R_k} \frac{x \tan(ax)}{x^2+b^2} \, \mathrm{d}x = \frac{\pi}{e^{2ab} + 1}. $$
For the interested user, here is the graph of $R \mapsto \operatorname{PV}\!\!\int_{0}^{R} \frac{x \tan(ax)}{x^2+b^2} \, \mathrm{d}x $ when $a = b = 1$:
Sketch of the 1st Proof. We find that $2\mathcal{I}$ can be written as
\begin{align*}
2\mathcal{I}
&= \lim_{k\to\infty} \operatorname{PV}\!\!\int_{-a R_k}^{a R_k} \frac{y \tan y}{y^2+(ab)^2} \, \mathrm{d}x \tag{$y = ax$}\\
&= \lim_{N\to\infty} \sum_{k=-N}^{N} \int_{-\pi/2}^{\pi/2} \frac{1}{2} \left( \frac{1}{y + k\pi + iab} + \frac{1}{y + k\pi - iab} \right) \tan y \, \mathrm{d}y
\end{align*}
Then by utilizing the partial fraction decomposition of the cotangent,
$$ \cot z = \lim_{N\to\infty} \sum_{k=-N}^{N} \frac{1}{z - \pi k}, $$
we get
$$ 2\mathcal{I}
= \int_{-\pi/2}^{\pi/2} \frac{\cot(y + iab) + \cot(y - iab)}{2} \tan y \, \mathrm{d}y. $$
By utilizing trigonometric identities, we find that the integrand simplifies as follows:
$$ \frac{\cot(y + iab) + \cot(y - iab)}{2} \tan y = \frac{\sin^2 y}{\cosh^2(ab) - \cos^2 y}. $$
Therefore
\begin{align*}
2\mathcal{I}
&= \int_{-\pi/2}^{\pi/2} \frac{\sin^2 y}{\cosh^2(ab) - \cos^2 y} \, \mathrm{d}y \\
&= \int_{-\infty}^{\infty} \left( \frac{1}{t^2+1} - \frac{\sinh^2(ab)}{\cosh^2(ab) t^2 + \sinh^2(ab)} \right) \, \mathrm{d}y \tag{$t = \tan y$} \\
&= \pi - \pi \tanh(ab)
= \frac{2\pi}{e^{2ab} + 1}.
\end{align*}
Therefore the desired equality follows.
Sketch of the 2nd Proof. Let me only sketch the proof of this equality. We will instead study
$$ 2\mathcal{I} = \lim_{k\to\infty} \operatorname{PV}\!\!\int_{-R_k}^{R_k} \frac{x \tan(ax)}{x^2+b^2} \, \mathrm{d}x. $$
To this end, let $0 < \epsilon < b < T$ and $R > 0$. If $\Gamma$ is a counter-clockwise rectangular contour with the four corners $\pm R + i\epsilon$ and $ \pm R + iT$, then
$$ \int_{\Gamma} \frac{z \tan(az)}{z^2+b^2} \, \mathrm{d}z = 2\pi i \, \underset{z=ib}{\mathrm{Res}} \, \frac{z \tan(az)}{z^2+b^2} = -\pi \tanh(ab). $$
Moreover, as $R \to \infty$ along the sequence $(R_k)_{k \geq 1}$ and $\epsilon \to 0^+$,
$$ \int_{-R+i\epsilon}^{R+i\epsilon} \frac{z \tan(az)}{z^2+b^2} \, \mathrm{d}z \quad \longrightarrow \quad 2\mathcal{I} - \pi i \lim_{k\to\infty} \sum_{\substack{\xi : |\xi| < R; \\ \tan(a\xi) = 0 }} \underset{z=\xi}{\mathrm{Res}} \, \frac{z \tan(az)}{z^2+b^2} = 2\mathcal{I}, $$
since
$$ \sum_{\substack{\xi : |\xi| < R; \\ \tan(a\xi) = 0 }} \underset{z=\xi}{\mathrm{Res}} \, \frac{z \tan(az)}{z^2+b^2}
= \sum_{\substack{\xi : |\xi| < R; \\ \tan(a\xi) = 0 }} \frac{a \xi \sec^2(a\xi)}{\xi^2+b^2}
= 0 $$
by the symmetry. Finally, consider the integral
$$ \int_{\gamma} \frac{z \tan(az)}{z^2+b^2} \, \mathrm{d}z $$
along the polygonal path $\gamma$ from $R+i\epsilon$ to $R+iT$ to $-R+iT$ to $-R+i\epsilon$. To study the behavior of this integral, we make the following observations:
$ \tan(z) \to i $ uniformly as $\operatorname{Im}(z) \to +\infty$.
The argument change along $\gamma$, viewed from any given point, is $\pi + o(1)$ as $R, T \to \infty$ and $\epsilon \to 0^+$
So we expect that
$$ \int_{\gamma} \frac{z \tan(az)}{z^2+b^2} \, \mathrm{d}z \approx \frac{i}{2} \int_{\gamma} \left( \frac{1}{z+ib} + \frac{1}{z-ib} \right) \, \mathrm{d}z \approx \frac{i}{2} (i\pi + i\pi) = -\pi $$
and the error in the above approximate equalities will vanish along the limit. Therefore
$$ 2\mathcal{I} - \pi = \lim_{\substack{R_k\to\infty \\ T \to \infty \\ \epsilon \to 0^+}} \int_{\Gamma} \frac{z \tan(az)}{z^2+b^2} \, \mathrm{d}z = -\pi \tanh(ab), $$
and solving this for $\mathcal{I}$ yields the desired equality.
I'm not sure how you showed the two integrals are equivalent, but the following is an evaluation of $$\int_{0}^{\infty} e^{-x} \cos (x) \tanh(x) \, \frac{\mathrm dx}{x}.$$
For $\Re(s) >0$, we have
$$ \begin{align} \int_{0}^{\infty} \tanh(t) e^{-st} \, \mathrm dt &= \int_{0}^{\infty} \frac{1-e^{-2t}}{1+e^{-2t}} \, e^{-st} \, \mathrm dt \\ &= \int_{0}^{\infty} (1-e^{-2t})e^{-st} \sum_{n=0}^{\infty} (-1)^{n}e^{-2tn} \, \mathrm dt \\ &= \sum_{n=0}^{\infty} (-1)^{n}\int_{0}^{\infty}\left(e^{-(2n+s)t} -e^{-(2n+s+2)t} \right) \, \mathrm dt \\ &= \sum_{n=0}^{\infty} \frac{(-1)^{n}}{2n+s} - \sum_{n=0}^{\infty} \frac{(-1)^{n}}{2n+s+2} \\ &= \frac{1}{2} \left(\sum_{n=0}^{\infty} \frac{(-1)^{n}}{n+\frac{s}{2}} - \sum_{n=0}^{\infty} \frac{(-1)^{n}}{n+\frac{s}{2}+1} \right) \\ &\overset{(1)}{=} \frac{1}{4} \left(\psi \left(\frac{s}{4} + \frac{1}{2}\right)- \psi \left( \frac{s}{4}\right) - \psi \left(\frac{s}{4}+1 \right) + \psi \left(\frac{s}{4}+\frac{1}{2} \right)\right) \\ &\overset{(2)}{=} \frac{1}{2} \left( \psi \left(\frac{s}{4}+ \frac{1}{2} \right) - \psi \left(\frac{s}{4} \right) - \frac{2}{s} \right). \end{align}$$
Therefore, $ \begin{align} \int_{0}^{\infty} e^{-x} \cos (x) \tanh(x) \, \frac{\mathrm dx}{x} &= \Re\int_{0}^{\infty} e^{-(1+i)x} \frac{\tanh (x)}{x} \, \mathrm dx \\ &= \Re \int_{0}^{\infty}e^{-(1+i)x} \tanh(x) \int_{0}^{\infty} e^{-xt} \, \mathrm dt \, \mathrm dx \\ &= \Re \int_{0}^{\infty} \int_{0}^{\infty}\tanh(x) e^{-(1+i+t)x} \, \mathrm dx \, \mathrm dt \\ &= \Re \int_{0}^{\infty} \frac{1}{2}\left(\psi \left(\frac{1+i+t}{4}+ \frac{1}{2} \right) - \psi \left(\frac{1+i+t}{4} \right) - \frac{2}{1+i+t} \right) \, \mathrm dt \\ &= \Re \left(2 \ln\Gamma \left(\frac{3+i+t}{4} \right) -2 \ln \Gamma \left(\frac{1+i+t}{4} \right) -\ln (1+i+t)\right) \Bigg|_{0}^{\infty} \\ &\overset{(3)}= \Re \left(-2 \ln(2) - 2 \ln \Gamma \left(\frac{3+i}{4} \right)+2 \ln \Gamma \left(\frac{1+i}{4} \right) + \ln(1+i)\right)\\ &= - \frac{3 \ln (2)}{2} + 2\Re \left(\ln \Gamma \left(\frac{1+i}{4} \right)- \ln \Gamma \left(\frac{3+i}{4} \right) \right). \end{align}$
$(1)$ For $\Re(z) > 0$, $\sum_{k=0}^{\infty} \frac{(-1)^{k}}{k+z} = \frac{1}{2} \left(\psi \left(\frac{z+1}{2} \right) - \psi \left(\frac{z}{2} \right) \right).$
$(2)$ Recurrence relation of the digamma function
$(3)$ For $a >0$, $\ln \Gamma(z_{1}+ax)- \ln \Gamma(z_{2}+ax) \sim (z_{1}-z_{2}) \ln(ax) + \mathcal{O} \left(\frac{1}{x}\right)$ as $x \to \infty$.
Best Answer
$$I = \int_{0}^{\infty} \frac{x \cos (a x)}{x^2+1} \coth (s x) \, dx$$
Recall the Mittag-Leffler pole expansion of $\cot (z)$: $$\cot (z) = \frac{1}{z} + 2z \sum_{k=1}^{\infty} \frac{1}{z^2 - (k \pi)^2}$$ Letting $z \to i s x$ and multiplying by $i$ gives: $$\coth (s x) = \frac{1}{s x} + 2s x \sum_{k=1}^{\infty} \frac{1}{(s x)^2 + (k \pi)^2}$$ Substituting this expression in, we have then: $$I = \frac{1}{s} \int_{0}^{\infty} \frac{\cos (a x)}{x^2+1} \, dx + \sum_{k=1}^{\infty} \int_{0}^{\infty} \frac{2sx^2 \cos (a x)}{(x^2+1)((s x)^2 + (k \pi)^2)} \, dx$$ We can separate the integrals because each converges individually. The first integral has a well-known result: $$\frac{1}{s} \int_{0}^{\infty} \frac{\cos (a x)}{x^2+1} \, dx = \frac{\pi e^{-a}}{2s}$$
The second integral is a little less trivial, but it can be done.
Mathematica gives: $$\int_{0}^{\infty} \frac{2sx^2 \cos (a x)}{(x^2+1)((s x)^2 + (k \pi)^2)} \, dx=\frac{\pi^2 k e^{-\frac{\pi a k}{s}}-\pi e^{-a}s}{\pi^2k^2-s^2}$$
Furthermore, Mathematica gives: $$\sum_{k=1}^{\infty} \frac{\pi^2 k e^{-\frac{\pi a k}{s}}-\pi e^{-a}s}{\pi^2k^2-s^2} = \frac{1}{2}e^{-\frac{\pi a}{s}}\left(\Phi\left(e^{-\frac{a\pi}{s}},1,1-\frac{s}{\pi}\right)+\Phi\left(e^{-\frac{a\pi}{s}},1,1+\frac{s}{\pi}\right)\right)-\frac{\pi e^{-a}}{2s}+\frac{\pi}{2} e^{-a}\cot(s)$$ Where $\Phi$ is the Lerch transcendent.
Thus summing the two results, we have:
$$I = \frac{1}{2}e^{-\frac{\pi a}{s}}\left(\Phi\left(e^{-\frac{a\pi}{s}},1,1-\frac{s}{\pi}\right)+\Phi\left(e^{-\frac{a\pi}{s}},1,1+\frac{s}{\pi}\right)\right)+\frac{\pi}{2} e^{-a}\cot(s)$$
However, this expression does not hold for $s$ being an integer multiple of $\pi$ due to $\cot (s)$ causing it to be an indeterminate when one takes the limit, so there needs to be some finesse for evaluating these cases. I'm not sure if there is simplification possible for the result I obtained for $I$ involving $\Phi$.