3.471.9
I will start with 3.471.9, since it is much simpler. The key ingredient
Slater theorem, saying that
$$
\int_0^\infty t^{\alpha-1} f(t) g\left(\frac{x}{t}\right) \mathrm{d} t = \frac{1}{2 \pi i} \int_{\gamma - i \infty}^{\gamma + i \infty} \mathcal{M}_f(s+\alpha) \mathcal{M}_g(s) x^{-s} \mathrm{d} s = \mathcal{M}^{-1} \left\{ \mathcal{M}_f(s+\alpha) \mathcal{M}_g(s) \right\}(x)
$$
where $\mathcal{M}_f(s) = \int_0^\infty t^{s-1} f(t) \mathrm{d} t$ is the Mellin transform of $f$, valid within some strip $\min < \Re(s) < \max $, and $\gamma$ is a real constant, such that $\gamma$ is within the strip of validity of the Mellin transform of $g$, and $\gamma + \Re(\alpha)$ in within the strip of the Mellin transform of $f$.
The proof is easy:
$$ \begin{eqnarray}
\int_0^\infty t^{\alpha-1} f(t) g\left(\frac{x}{t}\right) \mathrm{d} t &=& \int_0^\infty t^{\alpha-1} f(t) \left( \frac{1}{2 \pi i} \int_{\gamma - i \infty}^{\gamma+\infty} \mathcal{M}_g(s) x^{-s} t^s \mathrm{d} s \right) \mathrm{d} t \\ &=& \frac{1}{2 \pi i} \int_{\gamma - i \infty}^{\gamma+\infty} \left( \int_0^\infty t^{s+\alpha-1} f(t) \mathrm{d} t \right) \mathcal{M}_g(s) x^{-s} \mathrm{d} s \\
&=& \frac{1}{2 \pi i} \int_{\gamma - i \infty}^{\gamma+\infty} \mathcal{M}_f(s+\alpha) \mathcal{M}_g(s) x^{-s} \mathrm{d} s
\end{eqnarray}
$$
We now apply this theorem to 3.471.9, using Cahen-Mellin integral $\mathcal{M}_{\exp(-\bullet)}(s) = \Gamma(s) $:
$$\begin{eqnarray}
\int_0^\infty t^{\alpha-1} \exp\left(-y t\right) \exp\left( - \frac{x}{t}\right) \mathrm{d} t &=& \frac{1}{2 \pi i} \int_{\gamma - i \gamma}^{\gamma + i \infty} \Gamma(s+\alpha) \Gamma(s) y^{-\alpha -s} x^{-s} \mathrm{d}s \\ &=& y^{-\alpha} \left( \frac{1}{2 \pi i} \int_{\gamma - i \gamma}^{\gamma + i \infty} \Gamma(s+\alpha) \Gamma(s) \left( x y \right)^{-s} \mathrm{d}s \right) \\ &=& y^{-\alpha} \left( 2 (x y)^{\alpha/2} K_{\alpha}(2 \sqrt{x y} ) \right) \\
&=& 2 \left( \frac{x}{y} \right)^{\alpha/2} K_\alpha\left(2 \sqrt{x y}\right)
\end{eqnarray}
$$
where DMLF 10.43.19 was used.
6.635.3
A formula of kin is worked out in Bateman, Erdelyi et al, "Higher Transcendental Functions", chapter 7, section 7.7.6 on Macdonald's and Nochilson's formulas, formula (37). It says:
$$
\int_0^\infty \exp\left(-\frac{t}{2} - \frac{x^2+X^2}{2 t} \right) I_\nu \left( \frac{x X}{t} \right) \frac{\mathrm{d} t}{t} = \cases{2 I_\nu\left( x \right) K_\nu\left( X \right) & x < X \\ 2 I_\nu\left( X \right) K_\nu\left( x \right) & x > X }
$$
By doing an analytic continuation $x \to \mathrm{e}^{i \pi/2} x$ of the (the only possible) first branch, we get
$$
\int_0^\infty \exp\left(-\frac{t}{2} - \frac{X^2-x^2}{2 t} \right) J_\nu \left( \frac{x X}{t} \right) \frac{\mathrm{d} t}{t} = 2 J_\nu(x) K_\nu(X)
$$
Now, performing reparameterization $X = \sqrt{2 \alpha \left( \sqrt{\beta^2 + \gamma^2}+ \beta \right)}$, $x = \sqrt{2 \alpha \left( \sqrt{\beta^2 + \gamma^2} - \beta \right)}$, accompanied by a change of variables $t \to \left(2 \alpha t\right)^{-1}$ we arrive at 6.635.3.
Derivation of the point of departure formula relies on two auxiliary results:
$$
\int_0^\infty J_{\nu}(x u) J_{\nu}(X u) \exp\left(-\frac{t u^2}{2}\right) u \mathrm{d} u = \frac{1}{t} \exp\left(-\frac{x^2+X^2}{2t}\right) I_\nu\left( \frac{x X}{t}\right)
$$
which is obtained by expanding the exponential in series and integrating term-wise. Then, using the above integral representation in the formula under consideration, and carrying out simple integration with respect to $t$, we get
$$
\int_0^\infty \frac{\mathrm{e}^{-t/2}}{t} \exp\left(-\frac{x^2+X^2}{2t}\right) I_\nu\left( \frac{x X}{t}\right) \mathrm{d} t =
\int_0^\infty J_\nu(x u) J_\nu(X u) \frac{2 u \mathrm{d} u}{1+u^2}
$$
The resulting integral is of Sonine-Gegenabuer-type, and equals to $2 I_\nu( \min(x,X)) K_\nu(\max(x,X))$.
By exploiting the integral representation for the $J_0$ function:
$$ J_0(z) = \frac{1}{2\pi}\int_{0}^{2\pi}e^{iz\cos t}\,dt \tag{1}$$
we have:
$$ \int_{0}^{+\infty}\sin(ax)\, J_0(b\sqrt{1+x^2})\,dx =\frac{1}{2\pi}\int_{0}^{+\infty}\int_{0}^{2\pi}\sin(ax)e^{ib\sqrt{x^2+1}\cos t}\,dt\,dx \tag{2}$$
where:
$$ \sqrt{x^2+1}\cos t = \cos(t-t_0)-x\sin(t-t_0),\qquad x=\tan t_0 \tag{3}$$
so:
$$\begin{eqnarray*} I &=& \frac{1}{2\pi}\int_{0}^{+\infty}\int_{0}^{2\pi}\sin(ax)\,e^{ib\cos t-ibx\sin t}\,dt\,dx\\&=&\frac{1}{2\pi}\int_{0}^{2\pi}\frac{a}{a^2-b^2\sin^2 t}e^{ib\cos t}\,dt\\&=&\frac{1}{2\pi}\int_{0}^{2\pi}\frac{a}{a^2-b^2\sin^2 t}\,\cos(b\cos t)\,dt\tag{4} \end{eqnarray*}$$
but since:
$$ \frac{1}{2\pi}\int_{0}^{2\pi}\sin^{2n}(t)\cos(b\cos t)\,dt = \frac{ (2n-1)!!}{b^n}\, J_n(b) \tag{5}$$
we have:
$$\begin{eqnarray*} I &=& \frac{1}{a}\sum_{n\geq 0}\left(\frac{b}{a}\right)^{2n} \frac{1}{2\pi}\int_{0}^{2\pi}\sin^{2n}(t)\cos(b\cos t)\,dt \\&=&\frac{1}{a}\sum_{n\geq 0}\frac{(2n)!}{n!}\left(\frac{b}{2a^2}\right)^n J_n(b)\\&=&\frac{1}{a}\sum_{m=0}^{+\infty}\frac{(-1)^m}{m!}\left(\frac{b}{2}\right)^m\sum_{n\geq 0}\left(\frac{b}{2a}\right)^{2n}\frac{(2n)!}{n!(n+m)!}\\&=&\frac{1}{\sqrt{a^2-b^2}}+\frac{1}{a}\int_{0}^{1}\sum_{m\geq 1}\frac{(-1)^m}{m!(m-1)!}\left(\frac{b}{2}\right)^m(1-u)^{m-1}\sum_{n\geq 0}\binom{2n}{n}\left(\frac{b^2}{4a^2}\right)^{n}u^n\,du\\&=&
\frac{1}{\sqrt{a^2-b^2}}+\int_{0}^{1}\sum_{m\geq 1}\left(\frac{b}{2}\right)^m \frac{(-1)^m}{m!(m-1)!}\frac{(1-u)^{m-1}}{\sqrt{a^2-b^2 u}}\,du\\&=&
\frac{1}{\sqrt{a^2-b^2}}-\sqrt{\frac{b}{2}}\int_{0}^{1}\frac{J_1(\sqrt{2b(1-u)})}{\sqrt{(a^2-b^2 u)(1-u)}}\,du\\&=&\frac{1}{\sqrt{a^2-b^2}}-\sqrt{2b}\int_{0}^{1}\frac{J_1(t\sqrt{2b})}{\sqrt{(a^2-b^2)+b^2 t^2}}\,dt.\tag{6}\end{eqnarray*}$$
Best Answer
Note that for $r>0$ one has integral representation $$J_0(r)=\frac{1}{2\pi}\int_0^{2\pi}e^{ir\cos\phi}d\phi$$ Hence $$I=\int_0^{\infty}J_0\left(\alpha\sqrt{x^2+z^2}\right)\cos \beta x\,dx= \frac{1}{4\pi}\int_0^{2\pi}\int_{-\infty}^{\infty}e^{i\alpha\sqrt{x^2+z^2}\cos\phi}\cos\beta x \, d\phi \, dx.\tag{1}$$ On the other hand, $$\sqrt{x^2+z^2}\cos\phi=z\cos(\phi-\phi_0)+x\sin(\phi-\phi_0),$$ where $\tan\phi_0=-\frac{x}{z}$. Exchanging the order of integration in (1) and shifting $\phi$ by $\phi_0$, we arrive at $$I=\frac{1}{4\pi}\int_0^{2\pi}\int_{-\infty}^{\infty}e^{i\alpha(z\cos\phi+x\sin\phi)}\cos\beta x \, d\phi \, dx.$$ Finally, using that $\displaystyle\frac{1}{2\pi}\int_{-\infty}^{\infty}e^{i\omega x}dx=\delta(\omega)$ we obtain $$I=\frac{1}{4}\int_0^{2\pi}e^{i\alpha z\cos\phi}\Bigl[\delta\left(\alpha\sin\phi+\beta\right)+\delta\left(\alpha\sin\phi-\beta\right)\Bigr]d\phi$$ It remains to use $\delta(f(x))=\sum\limits_{\text{zeros of }f}\frac{1}{|f'(x_k)|}\delta(x-x_k)$ and compute the two contributions coming from each of the two delta-functions.