Following this link, Zacky calculated the integral $$\int_0^\infty x{J}_0(rx) \,\frac{e^{-z\sqrt{x^2+a^2}}}{\sqrt{x^2+a^2}}\, {\rm d}x=\frac{e^{-a\sqrt{r^2+z^2}}}{\sqrt{r^2+z^2}} \,.$$ It is a special case of the more general Bessel type integrals found in Watson.
I'm wondering if also a closed form expression for $$\int_0^\infty {J}_0(rx) \,\frac{e^{-z\sqrt{x^2+a^2}}}{\sqrt{x^2+a^2}} \, {\rm d}x$$ can be found.
I tried the method Zacky applied, but it doesn't seem to work, because the appearance of error-functions complicate matters. I also tried the approach from Watson, but only converted the improper integral into a definite integral. Equation 2f (Watson, p.416) gives, while omitting the factor $t^{\mu+1}$, for $\nu=1/2$ and $\mu=0$ $$\sqrt{\frac{\pi}{2a}}\int_0^\infty J_0(bt) \, \frac{e^{-a\sqrt{t^2+z^2}}}{\sqrt{t^2+z^2}} \, {\rm d}t = \frac12 \int_0^\infty {\rm d}u \int_0^\infty {\rm d}t \, J_0(bt) \, u^{-3/2}\, e^{-\frac{a}{2}\left(u+\frac{t^2+z^2}{u}\right)} \\ =\frac14 \sqrt{\frac{2\pi}{a}} \int_0^\infty \frac{{\rm d}u}{u} \, I_0\left(\frac{b^2u}{4a}\right) \, e^{-\frac{(2a^2+b^2)u}{4a}-\frac{z^2a}{2u}} $$ where the $t$-integral can be found here. Substituting $v=\frac{b^2u}{4a}$ then gives $$=\frac14 \sqrt{\frac{2\pi}{a}} \int_0^\infty \frac{{\rm d}v}{v} \, I_0(v) \, e^{-\frac{(2a^2+b^2)v}{b^2}-\frac{b^2z^2}{8v}} \, .$$ Using the representation $$I_0(v)=\frac1\pi \int_0^\pi {\rm d}\theta \, e^{-v\cos(\theta)}$$ finally yields $$=\frac{1}{\sqrt{8\pi a}} \int_0^\pi {\rm d}\theta \int_0^\infty \frac{{\rm d}v}{v} \, e^{-\frac{\left(a^2+b^2\cos^2(\theta/2)\right)2v}{b^2}-\frac{b^2z^2}{8v}} \\ \stackrel{w=\frac{\left(a^2+b^2\cos^2(\theta/2)\right)2v}{b^2}}{=} \frac{1}{\sqrt{8\pi a}} \int_0^\pi {\rm d}\theta \int_0^\infty \frac{{\rm d}w}{w} \, e^{-w-\frac{\left(z\sqrt{a^2+b^2\cos^2(\theta/2)}\right)^2}{4w}} \\ = \frac{1}{\sqrt{2\pi a}} \int_0^\pi {\rm d}\theta \, K_0\left(z\sqrt{a^2+b^2\cos^2(\theta/2)}\right)$$ resorting to this integral representation of the Bessel-K-function.
From Bateman's table of transforms, Vol. 2, page 31, #20, or Gradshteyn & Rhyzik, #6.637, we have $$\int_0^\infty \operatorname{J}_0(rx) \,\frac{e^{-z\sqrt{x^2+a^2}}}{\sqrt{x^2+a^2}} \, {\rm d}x =\operatorname{I}_{0}\left(\tfrac{a}{2}\left[\sqrt{r^2+z^2}-z\right]\right) \operatorname{K}_0\left(\tfrac{a}{2}\left[\sqrt{r^2+z^2}+z\right]\right). $$ I have tested it numerically and it works. Integrals of this type appear in the paper by Dixon & Ferrar, Integrals for the product of two Bessel functions (II) by A. L. Dixon & W. L. Ferrar, The Quarterly Journal of Mathematics, Volume os-4, Issue 1, 1933, Pages 297–304, at https://academic.oup.com/qjmath/article-abstract/os-4/1/297/1524326?redirectedFrom=fulltext.
Similar integrals, eqn. (4.7), also occur in the paper Analogue of a Fock-Type integral arising from electromagnetism and its applications in number theory by Atul Dixit and Arindam Roy, at https://arxiv.org/pdf/1907.03650.pdf.
The integral from Gradshteyn & Rhyzik, #6.645, is of the form $$ \int_{1}^{\infty}(x^2-1)^{-1/2}\,{\rm e}^{-\alpha x}J_{\nu}\left(\beta\sqrt{x^2-1}\right)\,{\rm d}x=I_{\nu/2}\left(\frac{1}{2}\left(\sqrt{\alpha^2+\beta^2}-\alpha\right)\right) K_{\nu/2}\left(\frac{1}{2}\left(\sqrt{\alpha^2+\beta^2}+\alpha\right)\right). $$ If we put $x=\cosh t,\,\, \alpha=\xi-z,\,\,\beta=2(z\xi)^{1/2},\,\,\nu\mapsto 2\nu,$ the integral reduces to the form of Dixon & Ferrar, $$ \operatorname{I}_{\mspace{1mu}\nu}(z)\operatorname{K}_{\nu}(\xi)= \int_{0}^{\infty}\operatorname{J}_{\mspace{1mu}2\nu}\left( 2(z\mspace{2mu}\xi)^{1/2}\operatorname{\sinh} t \right)\,{\rm e}^{-(\xi-z\mspace{1mu})\operatorname{\cosh}t}\,{\rm d}t\, . $$
Edit
The required integral is $$ I=\int_{0}^{\infty}J_{0}(rx)\,\frac{{\rm e}^{-z\sqrt{x^2+a^2}}}{\sqrt{x^2+a^2}}\,{\rm d}x\, , $$ where, as far as the integration is concerned, $r,z,a$ are considered to be constants.
We can replace the exponential in the integrand by using Poisson's integral, which is discussed at the end of $\S25$ in the paper by J.W.L.$\!$ Glaisher, On Riccati’s Equation and its Transformations, and on some Definite Integrals which satisfy them, in Phil. Trans. R. Soc. London, Page 792, Vol.$\!$ 172, 1881 https://royalsocietypublishing.org/doi/epdf/10.1098/rstl.1881.0018 or on page 186, $\S6.32,$ of Watson's Treatise, 2nd Edition, 1966. It can also be found on math.stackexchange: How to evaluate $\int_{0}^{+\infty}\exp(-ax^2-\frac b{x^2})\,dx$ for $a,b>0$ .
The Poisson integral representation is $$ \int_{0}^{\infty}{\rm e}^{-\alpha^2 q^2-\frac{\beta^2}{q^2}}\,{\rm d}q=\frac{\sqrt{\pi}}{2\alpha}{\rm e}^{-2\alpha\beta}. $$ This allows the integrand to be represented as $$ \frac{{\rm e}^{-z\sqrt{x^2+a^2}}}{\sqrt{x^2+a^2}}=\frac{1}{\sqrt{\pi}} \int_{0}^{\infty}{\rm e}^{-\frac{\left(x^2+a^2\right)q^2}{4}-\frac{z^2}{q^2}}\,{\rm d}q\, , $$ and the required integral becomes $$ I=\frac{1}{\sqrt{\pi}} \int_{0}^{\infty}{\rm d}q\,{\rm e}^{-\frac{a^2 q^2}{4}-\frac{z^2}{q^2}} \int_{0}^{\infty}{\rm d}x\,J_{0}(rx)\,{\rm e}^{-\frac{x^2 q^2}{4}}. $$ From #$10.22.52$ of DLMF, Digital Library of Mathematical Functions, or Watson, $\S13.3, \#5,$ $$ \int_{0}^{\infty}J_{0}(bt)\,{\rm e}^{-p^2t^2}\,{\rm d}t=\frac{\sqrt{\pi}}{2p}{\rm e}^{-\frac{b^2}{8p^2}}\,I_{0}\!\left(\frac{q^2 b^2}{8p^2}\right), $$ where $I_{0}(x)$ is the modified Bessel function of the first kind. This gives $$ I=\int_{0}^{\infty}\frac{{\rm d}q}{q}\,{\rm e}^{-\frac{a^2q^2}{4}-\frac{2z^2+r^2}{2q^2}}\,I_{0}\!\left(\frac{r^2}{2q^2}\right)\, . $$ Now put the integral representation ($\#32.10.1,$ DLMF) of the modified Bessel function of the first kind, $$ I_{0}\!\left(\frac{r^2}{2q^2}\right)=\frac{1}{2\pi}\int_{0}^{2\pi} {\rm e}^{\frac{r^2}{2q^2}\cos(\theta)}\,{\rm d}\theta\, , $$ normally defined over the interval $0$ to $\pi,$ but here represented over the interval $0$ to $2\pi,$ into the integrand to give $$ I=\frac{1}{2\pi}\int_{0}^{2\pi}\!\!\!{\rm d}\theta\int_{0}^{\infty} \frac{{\rm d}q}{q}\,{\rm e}^{-\frac{a^2q^2}{4}-\frac{2z^2+r^2-r^2\cos(\theta)}{2q^2}}. $$ Substituting $v=aq/2$ then gives $$ I=\frac{1}{2\pi}\int_{0}^{2\pi}{\rm d}\theta\int_{0}^{\infty} \frac{{\rm d}v}{v}\,{\rm e}^{ -v^2-\frac{\frac{a^2}{2}\left( 2z^2+r^2-r^2\cos(\theta)\right) }{4v^2} }. $$ The $v$ integral represents the Macdonald modified Bessel function of the second kind, and from Watson, $\S6.22, \#15,$ its integral representation is $$ K_{0}(\varpi)=\frac{1}{2}\int_{0}^{\infty}{\rm e}^{\tau-\frac{\varpi^2}{4\tau}}\,\frac{{\rm d}\tau}{\tau}\stackrel{\,\tau=p^2}{=} \int_{0}^{\infty}{\rm e}^{p^2-\frac{\varpi^2}{4p^2}}\,\frac{{\rm d}p}{p}\,. $$ Putting $ \varpi=\sqrt{\frac{a^2}{2}\Bigl(2z^2+r^2-r^2\cos(\theta)\Bigr)}, $ the integral then assumes the form $$ I=\frac{1}{2\pi}\int_{0}^{2\pi}{\rm d}\theta\, K_{0}\left(\sqrt{\frac{a^2}{2}\Bigl(2z^2+r^2-r^2\cos(\theta)\Bigr)} \,\, \right)\,. $$ If we introduce new constants defined by $$ \lambda=\frac{a}{2}\left(\sqrt{r^2+z^2}-z\right),\quad \mu=\frac{a}{2}\left(\sqrt{r^2+z^2}+z\right),\tag{1} $$ then we have $$ \sqrt{\vphantom{Y^a}\lambda\mu}=ar/2,\quad \mu-\lambda=az, $$ and this allows the quantity under the square root sign to be written as $$ \varpi=\sqrt{\frac{a^2}{2}\Bigl(2z^2+r^2-r^2\cos(\theta)\Bigr)}= \sqrt{\lambda^2+\mu^2-2\lambda\mu\cos(\theta)}\,. $$ By Graf's addition formula, Watson $\S11.3,$ $\#8,$ and A definite integral on a circle with Bessel functions, we can then write the integral as $$ I=\frac{1}{2\pi}\int_{0}^{2\pi}{\rm d}\theta\, K_{0}\left(\varpi\right)=\frac{1}{2\pi}\int_{0}^{2\pi}{\rm d}\theta\, \sum_{m=-\infty}^{\infty}I_{m}\left(\lambda\right)K_{m}\left(\mu\right) {\rm e}^{{\rm i}\mspace{1mu}m\mspace{1mu}\theta}\,. $$ The integration over $\theta\mspace{1mu}$ leaves only the term with $m=0,$ and finally, by substituting for the Graf triangle constants, $\lambda$ and $\mu$ from $(1),$ we are left with the final form $$ I=I_{0}\left( \frac{a}{2}\left(\sqrt{r^2+z^2} -z \right)\right) K_{0}\left( \frac{a}{2}\left(\sqrt{r^2+z^2} +z \right)\right)\,. $$ This derivation is similar to @Diger's, except that he replaced the exponential term in the required integral by a $K_{1/2}(z)=\frac{{\rm e}^{-z}}{z}\sqrt{\frac{\pi z}{2}}$ Bessel function integral representation. It may be possible to redefine the $a$'s and $b$'s in his last line into the $\lambda$'s and $\mu$'s.