A possible solution is to use Poisson's integral $$\newcommand{\sgn}{\operatorname{sgn}}J_2(z)=\frac{z^2}{3\pi}\int_0^\pi\cos(z\cos\theta)\sin^4\theta\,d\theta$$ and switch the integrations (which is easy, though not entirely trivial, to justify): $$f(r,t)=\frac{r^2}{12}\int_0^\pi\big(\sgn(t+r\cos\theta)+\sgn(t-r\cos\theta)\big)\sin^4\theta\,d\theta.$$
Assume $r,t$ are positive. If $r<t$, the sum of $\sgn$'s is $2$ for all $\theta$, and we get your result.
Otherwise, introducing $\phi=\arcsin(t/r)$, we get (simplifiable if desired) $$f(r,t)=\frac{r^2}{6}\int_{\pi/2-\phi}^{\pi/2+\phi}\sin^4\theta\,d\theta=\frac{r^2}{96}(12\phi+8\sin2\phi+\sin4\phi).$$
UPDATE:
A generalization of this integral appears on page 429 of the textbook A Treatise on the Theory of Bessel Functions, namely
$$\small \int_{0}^{\infty} \frac{x}{x^{2}+\alpha^{2}} \, J_{\mu} (\beta x) \left(\cos \left(\frac{\pi(\mu-\nu)}{2} \right)J_{\nu}(rx) + \sin \left(\frac{\pi(\mu-\nu)}{2} \right)Y_{\nu}(rx) \right) \, \mathrm dx = I_{\mu}(\beta \alpha) K_{\nu}(r \alpha) $$ where $r \ge \beta >0$, and $\mu$ and $\nu$ are nonnegative real parameters such that $\mu > \nu -2$.
As in my previous answer, we can exploit properties of the Hankel function of the first kind.
Let $H_{0}^{(1)}(z)$ be the Hankel function of first kind of order zero defined as $$H_{0}^{(1)}(z) = J_{0}(z) + i Y_{0}(z), $$ where $Y_{0}(z)$ is the Bessel function of the second kind of order zero.
The principal branch of $H_{0}^{(1)}(z)$ has a branch cut on the negative real axis.
Since $$Y_{0}(xe^{i \pi})= Y_{0}(x) + 2i J_{0}(x), \quad x >0,$$ it follows that
$$H_{0}^{(1)}(xe^{i \pi}) = - J_{0}(x) + iY_{0}(x), \quad x >0. $$
And since $\frac{x}{x^{2}+\alpha^{2}}$ is an odd function, we have $$\int_{0}^{\infty} \frac{x}{x^{2}+\alpha^{2}} \, J_{0}(x) J_{0}(rx) \, \mathrm dx = \frac{1}{2} \, \Re \int_{-\infty}^{\infty} \frac{x}{x^{2}+\alpha^{2}} \, H_{0}^{(1)} (x) J_{0}(rx) \, \mathrm dx,$$ where the integration form $-\infty$ to $0$ is done on the upper side of the branch cut.
Now let's integrate the function $$f(z) = \frac{z}{z^{2}+\alpha^{2}} \, H_{0}^{(1)}(z) J_{0}(rz), \quad 0 < r \le 1 \, , $$ around a contour consisting of the upper side of the branch cut from $-R$ to $-\epsilon$, a small clockwise- oriented semicircle about the origin, the real axis from $\epsilon$ to $R$, and the upper half of the circle $|z|=R$.
Since $\lim_{z \to 0} f(z) = 0$, the contribution from the small semicircle about the origin vanishes as $\epsilon \to 0$.
And as $|z| \to \infty$ in the upper half plane, $|f(z)|$ is asymptotic to $$\frac{1}{\pi \sqrt{r}} \frac{e^{(r-1)\Im(z)}}{|z|^{2}}. $$
(See here.)
Since $0 < r \le 1$, the integral vanishes on the upper half of the circle $|z|=R$ as $R \to \infty$ (by the estimation lemma), and we have $$ \begin{align} \int_{-\infty}^{\infty} \frac{x}{x^{2}+\alpha^{2}} \, H_{0}^{(1)} (x) J_{0}(rx) \, \mathrm dx &= 2 \pi i \operatorname*{Res}_{z=i \alpha} f(z) \\ &= 2 \pi i \lim_{z \to i \alpha} \frac{z}{z+i \alpha} \, H_{0}^{(1)}(z) J_{0}(rz) \\ &= i \pi \, H_{0}^{(1)}(i \alpha)J_{0}(i r \alpha) \\ &\overset{(\spadesuit)}{=} i \pi \left(\frac{2K_{0}(\alpha)}{i \pi} \right) I_{0}(r \alpha) \\&= 2 K_{0}(\alpha) I_{0}(r \alpha). \end{align}$$
Equating the real parts on both sides of the equation, we have $$\int_{0}^{\infty} \frac{x}{x^{2}+ \alpha^{2}} \, J_{0}(x) J_{0}(rx) \, \mathrm dx = K_{0}(\alpha) I_{0}(r \alpha), \quad 0 < r \le 1. $$
This result also holds for $r=0$.
$\spadesuit$ https://dlmf.nist.gov/10.27#E8
To show that $$\int_{0}^{\infty} \frac{x}{x^{2}+ \alpha^{2}} \, J_{0}(x) J_{0}(rx) \, \mathrm dx = I_{0}(\alpha) K_{0}(r \alpha), \quad r \ge 1, $$ integrate $$g(z) = \frac{z}{z^{2}+\alpha^{2}} \, J_{0}(z) H_{0}^{(1)}(rz) $$ around the same contour.
Best Answer
It's easier to start with Bessel's Integral and apply residue theorem. For simplicity make a change of variable to ensure $r = 1$. Then, we have
\begin{align*}I = \int_0^\infty e^{-(a-it)x} J_0(x)\,dx &= \frac{1}{2\pi}\int_0^{\infty} \int_{-\pi}^{\pi} e^{-(a-it)x} e^{ix\sin \tau}\,d\tau\,dx \\&= \frac{1}{2\pi} \int_{-\pi}^{\pi} \int_0^{\infty} e^{-(a-i(t+\sin \tau))x} \,dx \,d\tau \\&= \frac{1}{2\pi} \int_{-\pi}^{\pi} \frac{1}{a-i(t+\sin \tau)} \,d\tau \\&= \frac{1}{2\pi} \int_{|z| = 1} \frac{1}{a-i\left(t+\frac{1}{2i}\left(z - \frac{1}{z}\right)\right)} \,\frac{dz}{iz}\tag{$\ast$}\\&= \frac{i}{\pi} \int_{|z| = 1} \frac{1}{z^2 - 2(a-it)z - 1} \,dz\end{align*}
where, in step $(\ast)$ we made the change of variable $z = e^{i\tau}$.
The function has two simple poles at $\alpha^{\pm} = a-it \pm \sqrt{(a-it)^2 + 1}$ (where we are using the principle branch of square root since $a > 0$) and only $\alpha^{-}$ is inside the unit disk. Therefore, the integral is $$I = \frac{i}{\pi}\frac{2\pi i}{\alpha^- - \alpha^+} = \frac{1}{\sqrt{(a-it)^2 + 1}} = \frac{1}{\sqrt{a^2 - t^2 + 1 - 2iat}}.$$
Now, comparing imaginary parts on both sides should give the result.