Can the integral $\int_0^\infty\ \frac{1}{\sqrt{x^2+a^2}} \cdot \frac{1}{1+e^\sqrt{x^2+a^2}}\,dx$ be solved analytically? I'm trying to solve the following integrals:

$$I_1(a)=\int_0^\infty\ \frac{1}{\sqrt{x^2+a^2}} \cdot \frac{1}{1+e^\sqrt{x^2+a^2}}\,dx$$
and
$$I_2(a)=\int_0^\infty\ \frac{1}{1+e^\sqrt{x^2+a^2}}\,dx.$$

They seem similar to the two integrals solved in this post, so I wonder if there is something in the Bessel function family that will do the trick. But so far, substitutions like $x=a \sinh(t)$ and $u=\sqrt{x^2+a^2}$ are not giving me integrals that can be manipulated into one of the Bessel function forms on Wikipedia or Wolfram.
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&lt;a&lt;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.
yes and very easy
let u=sqrt(x^2+a^2) 
then the integral becomes
int(2/(1+exp(u))
multiply numerator and denominator by exp(-u).
the derivates of denominator appear in numerator. and so on