Integral involving exponentials of cosh functions

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)=intlimits_0^infty dx e^{-px} K_0left[p(x^2+1)right] $$

Using the asymptotic expansion for $K_0$

$$ K_0(z)simsqrt{frac{pi}{2z}}e^{-z}sumlimits_{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)} sumlimits_{n=0}^infty c_np^{-n} $$

But this simplifies because the sum is known

$$ K_0(p)simsqrt{frac{pi}{2p}}e^{-p}sumlimits_{n=0}^infty c_np^{-n} , p to infty $$

So we get

$$ I(p)sim intlimits_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 sqrt{frac{pi}{p}}e^{p/4}operatorname{erfc}left(frac{sqrt{p}}{2} right) K_0(p) , p to infty $$

Here is a plot of the approximation versus numeric result for 'large' values of $p$ less than $1$:

By Fubinis theorem: $$I_p = intlimits_0^{+infty},intlimits_{-infty}^{+infty},e^{-p,[c+(c^2+1)cosh x]},mathrm{d}x,mathrm{d}c =intlimits_{-infty}^{+infty},intlimits_{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: $$intlimits_{0}^{+infty},e^{-p,[c+(c^2+1)cosh x]},mathrm{d}c = -dfrac{sqrt{{pi}}mathrm{e}^{frac{p}{4coshleft(xright)}-pcoshleft(xright)}left(operatorname{erf}left(frac{sqrt{p}}{2sqrt{coshleft(xright)}}right)-1right)}{2sqrt{pcoshleft(xright)}}$$ I am afraid that we can't solve the second integral analytically.