The $I_2(a)$ integral can't be reduced to a say, a Bessel function with an argument that is a function of $a.$ I say this because I worked out the asymptotics to be
$$ I_2(a) \sim \sqrt{\pi a/2} \ \frac{1}{1+e^a}, \quad a \to \infty $$
$$ I_2(a) \sim \log(2) + \frac{a^2}{4}(\kappa + \log{a}) +..., a \to 0$$
where
$$ \kappa = \gamma - 1/2 + \frac{5}{6} \ \log{2} + 12\ \zeta'(-1) = -1.3302154...$$
Putting the constant question aside, $I_2(a)$ has a functional form of $a \ K_1(a) $ for $a \to 0$ and $a K_{1/2}(a) $ for $a \to \infty,$ where $a$ is the McDonald function, a fundamental obstacle to seeking a simple form. It may be possible to get something as an infinite series of Bessel functions, but the integral is nice, so why bother.
Added later, sketch of how to find asymptotic formulas:
For $a \to \infty, \exp(\sqrt{x^2+a^2)} \sim \exp(a) \ \exp(x^2/(2a)). $ Then
$$ I_2(a) \sim \int_0^\infty dx \sum_{k=1}^\infty (-1)^{k-1}\Big(e^{-a} e^{-x^2/(2a)} \Big)^k = -\sum_{k=1}^\infty \big(-e^{-a} \big)^k \int_0^\infty dx \ e^{-x^2 k /(2a)} $$
$$ \quad = -\sqrt{\pi a/2} \sum_{k=1}^\infty \frac{ \big(-e^{-a} \big)^k}{\sqrt{k}} = -\sqrt{\pi a/2} \, \text{Li}_{1/2}(-e^{-a}) $$
The extra complexity of the polylogarithm is probably not worth the effort; in fact $e^{-a}$ is sufficient, but numerically a better appx to $ -\text{Li}_{1/2}(-e^{-a})$ is $1/(1+e^{-a}).$ By happenstance, the $e^{-a}$ form is better than the polylogarithm form for $1<a<3.$
For $a \to 0, \exp(\sqrt{x^2+a^2)} \sim \exp(x) \ \exp(a^2/(2x)) $ Then
$$ I_2(a) \sim \int_0^\infty dx \sum_{k=1}^\infty (-1)^{k-1}\Big(e^{-x}
\ \exp(-a^2/(2x)) \Big)^k =\sqrt{2}a \sum_{k=1}^\infty (-1)^{k-1}K_1(\sqrt{2} k a)$$ where again $\Sigma$ and $\int$ has been interchanged and the integral done in closed form in terms of the McDonald function. Now expand the McDonald function as $a \to 0.$ Then
$$ \sqrt{2} a K_1(\sqrt{2}ka) \sim \frac{1}{k} + k a^2\Big(\gamma - 1/2 +
\big(\log{a} + \log{(k/\sqrt{2}}) \big) \Big) + ... $$
The first sum, $\sum_{k=1}^\infty (-1)^{k-1}/k = \log{2},$ is easy. The second sum does not converge, but it can be made to have meaning within the context of zeta function regularization. I don't want to spend the time trying to justify the arguments. Instead, I numerically calculated to see that the result gave good agreement.