In this paper by VNP Anghel, a result obtained by Glasser is generalized (eq. 32):
\begin{equation}
I_{\frac{m+1}{2}}(b)=\left( \frac{b}{2\pi \sin^m\alpha}\right)^{1/2} \int_0^\pi \exp(b\cos\alpha\cos\theta)I_{\frac m2}(b\sin\alpha\sin\theta)(\sin\theta)^{\frac{m+2}{2}}\,d\theta \tag{1}\label{eq1}
\end{equation}
Then, in the case $A=B=b$, by choosing $m=0$, we have directly
\begin{align}
J(b,b)&=\int_0^{2 \pi} \int_0^{\pi} e^{b \cos \phi \cos \theta} \, I_0(b \sin \phi \sin \theta) \sin \theta \, d\theta \, d\phi\\
&=\int_0^{2 \pi} \, d\phi \left( \frac{2\pi}{b} \right)^{1/2}I_{1/2}(b)\\
&=4\pi\frac{\sinh b}{b}
\end{align}
The result \eqref{eq1} may also be used in the case $A\ne B$ to express the integral as a series. We remark first that $A$ can be chosen to be positive (parity of the integral wrt $A$ is clear from the substitution $\phi\to \pi-\phi$ in the integral). By the multiplication theorem,
\begin{equation}
I_{\nu}\left(\lambda z\right)=\lambda^{\nu}\sum_{k=0}^{\infty}
\frac{(\lambda^{2}-1)^{k}(\frac{1}{2}z)^{k}}{k!}I_{\nu+ k}\left(z
\right)
\end{equation}
with $\nu=0,\lambda=B/A,z=A\sin\phi\sin\theta$, one can express
\begin{equation}
I_0\left(B\sin\phi\sin\theta\right)=\sum_{k=0}^{\infty}
\frac{\left( B^2-A^2 \right)^k}{2^kk!A^k}\sin^k\phi\sin^k\theta \,I_{ k}\left(A\sin\phi\sin\theta\right)
\end{equation}
and thus, by interverting integral and summation,
\begin{align}
J(A,B)&=\int_0^{2 \pi} \int_0^{\pi} e^{A \cos \phi \cos \theta} \, I_0(B \sin \phi \sin \theta) \sin \theta \, d\theta \, d\phi\\
&=\sum_{k=0}^{\infty}\frac{\left( B^2-A^2 \right)^k}{2^kk!A^k}\int_0^{2 \pi}\sin^k\phi\,d\phi
\int_0^{\pi} e^{A \cos \phi \cos \theta}I_{ k}\left(A\sin\phi\sin\theta\right) \sin^{k+1}\theta \, d\theta
\end{align}
With $m=2k$ in eq. \eqref{eq1}, it may be expressed as
\begin{align}
J(A,B)&=\sum_{k=0}^{\infty}\frac{\left( B^2-A^2 \right)^k}{2^kk!A^k}\int_0^{2 \pi}\left( \frac{2\pi\sin^{2k}\phi}{A} \right)^{1/2}\sin^k\phi \,I_{k+1/2}(A)\,d\phi\\
&=\sqrt{\frac{2\pi}{A}}\sum_{k=0}^{\infty}\frac{\left( B^2-A^2 \right)^k}{2^kk!A^k}I_{k+1/2}(A)\int_0^{2 \pi}\sin^{2k}\phi\,d\phi\\
&=\frac{(2\pi)^{3/2}}{\sqrt{A}}\sum_{k=0}^{\infty}\frac{(2k)!}{2^{3k}(k!)^3}\frac{\left( B^2-A^2 \right)^k}{A^k}I_{k+1/2}(A)
\end{align}
Numerical experiments show that this series converges very quickly, which is related to the asymptotic expansion of the modified Bessel function for large orders (DLMF).
I shall discuss the asymptotic properites of the integral when $r$ becomes large. Using the relations
$$\tag{1}
K'_0 (r) = - K_1 (r),\quad K'_1 (r) = -K_0 (r) - \frac{{K_1 (r)}}{r}
$$
and performing integration by parts, we find
\begin{align*}
\int_r^{ + \infty } {xK_0 (x){\rm e}^{ - a^2 x^2 } {\rm d}x} & = \frac{1}{{2a^2 }}{\rm e}^{ - a^2 r^2 } K_0 (r) - \frac{1}{{2a^2 }}\int_r^{ + \infty } {{\rm e}^{ - a^2 x^2 } K_1 (x){\rm d}x}
\\ & = \frac{1}{{2a^2 }}{\rm e}^{ - a^2 r^2 } K_0 (r) - \frac{1}{{4a^4 }}\frac{1}{r}{\rm e}^{ - a^2 r^2 } K_1 (r) \\ & \quad\,+ \frac{1}{{4a^4 }}\int_r^{ + \infty } {\frac{1}{x}{\rm e}^{ - a^2 x^2 } K_0 (x){\rm d}x} + \frac{1}{{2a^4 }}\int_r^{ + \infty } {\frac{1}{{x^2 }}{\rm e}^{ - a^2 x^2 } K_1 (x){\rm d}x} .
\end{align*}
Continuing in this way, it is suspected that the integral has an asymptotic expansion of the form
$$\tag{2}
\int_r^{ + \infty } {xK_0 (x){\rm e}^{ - a^2 x^2 } {\rm d}x} \sim \frac{1}{{2a^2 }}{\rm e}^{ - a^2 r^2 } \left[ {K_0 (r)\sum\limits_{n = 0}^\infty {\frac{{A_n (a)}}{{r^{2n} }}} - \frac{1}{{2a^2 }}\frac{{K_1 (r)}}{r}\sum\limits_{n = 0}^\infty {\frac{{B_n (a)}}{{r^{2n} }}} } \right]
$$
as $r\to +\infty$, where $A_n (a)$ and $B_n (a)$ are polynomials in $a^{-2}$. Differentiating both sides of $(2)$, using $(1)$, and equating like powers of $r$, we find that $
A_0 (a) = B_0 (a) = 1$ and
$$
A_n (a) = - \frac{{n - 1}}{{a^2 }}A_{n - 1} (a) + \frac{1}{{4a^4 }}B_{n - 1} (a),\quad B_n (a) = A_n (a) - \frac{n}{{a^2 }}B_{n - 1} (a)
$$
for $n\geq 1$.
To prove the asymptotic character of the expansion $(2)$, one may introduce a remainder term and show that it has the right order of magnitude. I am not going to pursue this here, for there is a simpler asymptotic expansion for the integral. Indeed, introducing the known large-$r$ asymptotic series of $K_{0,1}(r)$ in $(2)$ shows that there should be an asymptotic expansion of the form
$$\tag{3}
\int_r^{ + \infty } {xK_0 (x){\rm e}^{ - a^2 x^2 } {\rm d}x} \sim \frac{1}{{2a^2 }}\sqrt {\frac{\pi }{{2r}}} {\rm e}^{ - a^2 r^2 -r} \sum\limits_{n = 0}^\infty {\frac{{C_n (a)}}{{r^n }}}
$$
as $r\to +\infty$, where the $C_n (a)$ are polynomials in $a^{-2}$. Differentiating both sides of $(3)$, using the asymptotic expansion of $K_0(r)$, and equating like powers of $r$, we find that
$$
C_0 (a) = 1,\quad C_1 (a) = - \frac{1}{{2a^2 }} - \frac{1}{8}
$$
and
$$
C_n (a) = ( - 1)^n \frac{{(2n)!^2 }}{{32^n n!^3 }} - \frac{1}{{2a^2 }}C_{n - 1} (a) - \frac{{2n - 3}}{{4a^2 }}C_{n - 2} (a)
$$
for $n\geq 2$. It is seen that $C_n(a)$ is a polynomial in $a^{-2}$ of degree $n$.
Finally, I shall demonstrate the asymptotic property of the expansion in $(3)$. For each $N\geq 0$ and $r \gg 0$, we write
$$\tag{4}
\int_r^{ + \infty } {xK_0 (x){\rm e}^{ - a^2 x^2 } {\rm d}x} = \frac{1}{{2a^2 }}\sqrt {\frac{\pi }{{2r}}} {\rm e}^{ - a^2 r^2 - r} \sum\limits_{n = 0}^{N - 1} {\frac{{C_n (a)}}{{r^n }}} + R_N (a,r)
$$
We know that, for each $N\geq 0$,
$$\tag{5}
K_0 (r) = \sqrt {\frac{\pi }{{2r}}} {\rm e}^{ - r} \sum\limits_{n = 0}^{N - 1} {( - 1)^n \frac{{(2n)!^2 }}{{32^n n!^3 r^n }}} + \varepsilon _N (r),
$$
with
$$
\varepsilon _N (r) = \frac{{{\rm e}^{ - r} }}{{\sqrt r }}\mathcal{O}\!\left( {\frac{1}{{r^N }}} \right)
$$
for $r\gg 0$. Differentiating both sides of $(4)$, using $(5)$, and simplifying the expression gives
\begin{align*}
\sqrt {\frac{{2r}}{\pi }} \frac{1}{r}{\rm e}^{a^2 r^2 + r} R'_N (a,r) = & - \sqrt {\frac{{2r}}{\pi }} {\rm e}^r \varepsilon _N (r) + \frac{1}{{2a^2 }}\frac{{C_{N - 1} (a)}}{{r^N }} + \frac{1}{{4a^2 }}\frac{{(2N - 3)C_{N - 2} (a)}}{{r^N }} \\ & + \frac{1}{{4a^2 }}\frac{{(2N - 1)C_{N - 1} (a)}}{{r^{N + 1} }} = \mathcal{O}\!\left( {\frac{1}{{r^N }}} \right)
\end{align*}
as $r\to +\infty$. It is clear from $(4)$ that $\lim_{r\to+\infty}R_N (a,r)=0$, whence
$$
R_N (a,r) = \mathcal{O}(1)\int_r^{ + \infty } { \sqrt t {\rm e}^{ - a^2 t^2 - t} \frac{1}{{t^N }}\mathrm{d}t}
$$
as $r\to +\infty$. By L'Hôpital's rule
$$
\int_r^{ + \infty } {\sqrt t {\rm e}^{ - a^2 t^2 - t} \frac{1}{{t^N }}\mathrm{d}t} \sim \frac{1}{{2a^2 }}\frac{{{\rm e}^{ - a^2 r^2 - r} }}{{\sqrt r }}\frac{1}{{r^N }}
$$
and therefore
$$
R_N (a,r) = \frac{{{\rm e}^{ - a^2 r^2 - r} }}{{\sqrt r }}\mathcal{O}\!\left( {\frac{1}{{r^N }}} \right)
$$
as $r\to +\infty$. This proves that for each $N\geq 0$,
$$
\int_r^{ + \infty } {xK_0 (x){\rm e}^{ - a^2 x^2 } {\rm d}x} = \frac{1}{{2a^2 }}\sqrt {\frac{\pi }{{2r}}} {\rm e}^{ - a^2 r^2 - r}\! \left(\sum\limits_{n = 0}^{N - 1} {\frac{{C_n (a)}}{{r^n }}} +\mathcal{O}\!\left( {\frac{1}{{r^N }}} \right)\right)
$$
as $r\to +\infty$. $\quad \square$
Best Answer
These two integrals are classical results of integral involving Bessel functions. You may find methods in an very old but nice book, Chapter 13 "A treatise on the theory of Bessel functions" by Prof G.N.Watson.
In particular, the first one is more complicated than the second one.
Have fun.