Integral involving exponentials of cosh functions I recently encountered the following integral
$$I_p\equiv \int\limits_0^{+\infty}\,\int\limits_{-\infty}^{+\infty}\,e^{-p\,[c+(c^2+1)\cosh x]}\,\mathrm{d}x\,\mathrm{d}c\,,$$
with $p>0$. Can $I_p$ be calculated in terms of elementary functions?
Update:
Using the following integral representation for the Bessel function
$$
\int\limits_0^{+\infty}e^{-a\cosh x}\;\mathrm{d}x=K_0(a)\,,
$$
one can bring $I_p$ to the following form:
$$
I_p=2 \int\limits_0^{+\infty}e^{-c\,p}\;K_0\left[\left(c^2+1\right) p\right]\;\mathrm{d}c\,.
$$
Which shows that $I_p$ can be seen as a Laplace transform. However, I haven't managed to do the last integral either.
 A: By Fubinis theorem:
$$I_p = \int\limits_0^{+\infty}\,\int\limits_{-\infty}^{+\infty}\,e^{-p\,[c+(c^2+1)\cosh x]}\,\mathrm{d}x\,\mathrm{d}c =\int\limits_{-\infty}^{+\infty}\,\int\limits_{0}^{+\infty}\,e^{-p\,[c+(c^2+1)\cosh x]}\,\mathrm{d}c\,\mathrm{d}x.$$
We can solve the integral with respect to $c$ analytically:
$$\int\limits_{0}^{+\infty}\,e^{-p\,[c+(c^2+1)\cosh x]}\,\mathrm{d}c = -\dfrac{\sqrt{{\pi}}\mathrm{e}^{\frac{p}{4\cosh\left(x\right)}-p\cosh\left(x\right)}\left(\operatorname{erf}\left(\frac{\sqrt{p}}{2\sqrt{\cosh\left(x\right)}}\right)-1\right)}{2\sqrt{p\cosh\left(x\right)}}$$
I am afraid that we can't solve the second integral analytically.
A: Given the mess that results in doing just one of the integrals, I think it would be shocking if your integral had a closed form in elementary functions. We will find the large $p$ behavior. Let
$$
I(p)=\int\limits_0^\infty dx \ e^{-px} K_0\left[p(x^2+1)\right]
$$
Using the asymptotic expansion for $K_0$
$$
K_0(z)\sim\sqrt{\frac{\pi}{2z}}e^{-z}\sum\limits_{n=0}^\infty c_n z^{-n} \ \ , \ z \to \infty
$$
We will not need the form of the $c_n$ here. When we substitute $z=p(x^2+1)$, we will get a Gaussian $e^{-p(x^2+1)}$; we will use (essentially) Laplace's method to do the integral for large $p$. When we do, all the factors of $(x^2+1)$ outside the exponent will be evaluated at the maxima $x=0$. We will be left with
$$
\sqrt{\frac{\pi}{2p}}e^{-p(x^2+1)} \sum\limits_{n=0}^\infty c_np^{-n}
$$
But this simplifies because the sum is known
$$
K_0(p)\sim\sqrt{\frac{\pi}{2p}}e^{-p}\sum\limits_{n=0}^\infty c_np^{-n} \ \ , \ \ p \to  \infty
$$
So we get
$$
I(p)\sim \int\limits_0^\infty dx \ e^{-px} e^{p-p(x^2+1)} K_0(p) \ \ , \ \ p \to \infty
$$
The integral is an error function
$$
I(p) \sim \frac{1}{2}\sqrt{\frac{\pi}{p}}e^{p/4}\operatorname{erfc}\left(\frac{\sqrt{p}}{2} \right) K_0(p) \ \ , \ \ p \to \infty
$$
We can simplify further (thanks @Gary). Expand $K_0$ and $\operatorname{erfc}$ for large arguments, yielding
$$
I(p) \sim \sqrt{\frac{\pi}{2}} p^{-3/2}e^{-p} \quad , \quad p \to \infty
$$
Here is a plot of the approximations versus numeric result for 'large' values of $p$ less than $1$:

