With some differences I will point out along the way, the development of the solution of the above 2D heat equation follows the solution of the 1D heat equation for the semi-infinite rod. If you don't mind, I will develop the solution from the beginning, with my own notation which I will define below. The result is not pretty, but is at least in a calculable form.
The Equation
The problem is to solve the heat equation in the region $x^2+y^2 \gt 1$, where the temperature $u$ in the unit circle is held fixed at unity. Because there is no angular dependence, the equation becomes 1D in the radial coordinate $r=\sqrt{x^2+y^2}$:
$$\frac{\partial^2 u}{\partial r^2} + \frac{1}{r} \frac{\partial u}{\partial r} = \frac{\partial u}{\partial t}$$
$$u(1,t)=1$$
$$\lim_{r \to \infty} u(r,t)=0$$
$$u(r,0)=0 \quad r >1$$
The solution
The above equation lends itself to a Laplace transform solution. Define
$$\hat{u}(r,s) = \int_0^{\infty} dt \, e^{-s t} \, u(r,t)$$
The Laplace transform turns the PDE above into an ODE, as follows:
$$\frac{d^2 \hat{u}}{d r^2} + \frac{1}{r} \frac{d \hat{u}}{d r} - s \hat{u}=0$$
$$\hat{u}(1,s) = \frac{1}{s}$$
$$\lim_{r \to \infty} \hat{u}(r,s)=0$$
The general solution to the first equation is
$$\hat{u}(r,s) = A I_0(\sqrt{s} r) + B K_0(\sqrt{s} r)$$
where $I_0$ and $K_0$ are the modified Bessel functions of the first and second kind, respectively. Applying the boundary conditions, we find that the solution to the above equation becomes
$$\hat{u}(r,s) = \frac{K_0(\sqrt{s} r)}{s K_0(\sqrt{s})}$$
as shown in the problem statement. The solution to the heat equation is obtained upon performing an inverse Laplace transform:
$$u(r,t) = \frac{1}{i 2 \pi} \int_{c-i \infty}^{c+i \infty} \frac{ds}{s} e^{s t} \frac{K_0(\sqrt{s} r)}{K_0(\sqrt{s})}$$
The inverse transform
As the OP pointed out, the above is not a very useful way to express the solution. We therefore will use Cauchy's theorem and define a contour in the complex plane over which we may express the above ILT in terms of the value of an integral over other parts of the contour.
Consider the following integral in the complex $z$ plane:
$$\oint_C \frac{dz}{z} e^{t z} \frac{K_0(\sqrt{z} r)}{K_0(\sqrt{z})}$$
where $C$ is the following contour:
We will define $\text{Arg}{z} \in (-\pi,\pi]$, so the branch is the negative real axis. There are $6$ pieces to this contour, $C_k$, $k \in \{1,2,3,4,5,6\}$, as follows.
$C_1$ is the contour along the line $z \in [c-i R,c+i R]$ for some large value of $R$.
$C_2$ is the contour along a circular arc of radius $R$ from the top of $C_1$ to just above the negative real axis.
$C_3$ is the contour along a line just above the negative real axis between $[-R, -\epsilon]$ for some small $\epsilon$.
$C_4$ is the contour along a circular arc of radius $\epsilon$ about the origin.
$C_5$ is the contour along a line just below the negative real axis between $[-\epsilon,-R]$.
$C_6$ is the contour along the circular arc of radius $R$ from just below the negative real axis to the bottom of $C_1$.
We will show that the integral along $C_2$and $C_6$ vanish in the limits of $R \rightarrow \infty$ and $\epsilon \rightarrow 0$.
First of all, please note that the Bessel function in the integrand has the following asymptotic behavior:
$$K_0(w) \sim \sqrt{\frac{\pi}{2 w}} e^{-w} \quad (w \to \infty)$$
On $C_2$, the real part of the argument of the exponential is
$$R t \cos{\theta} - (r-1)\sqrt{R} \cos{\frac{\theta}{2}}$$
where $\theta \in [\pi/2,\pi)$. Clearly, $\cos{\theta} < 0$ and $\cos{\frac{\theta}{2}} > 0$, so that the integrand exponentially decays as $R \rightarrow \infty$ and therefore the integral vanishes along $C_2$.
On $C_6$, we have the same thing, but now $\theta \in (-\pi,-\pi/2]$. This means that, due to the evenness of cosine, the integrand exponentially decays again as $R \rightarrow \infty$ and therefore the integral also vanishes along $C_6$.
Thus, we are left with the following by Cauchy's integral theorem (i.e., no poles inside $C$):
$$\left [ \int_{C_1} + \int_{C_3} + \int_{C_4} + \int_{C_5}\right] \frac{dz}{z} \: e^{z t} \frac{K_0(\sqrt{z} r)}{K_0(\sqrt{z})} = 0$$
On $C_4$, note that the integral is finite and nonzero as $\epsilon$ in the limit $\epsilon \rightarrow 0$:
$$\int_{C_4} \frac{dz}{z} \exp{(t z)} \frac{K_0(\sqrt{z} r)}{K_0(\sqrt{z})} \sim i \int_{\pi}^{-\pi} d\phi = -i 2 \pi $$
On $C_3$, we parametrize by $z=e^{i \pi} x$ and the integral along $C_3$ becomes
$$\int_{C_3} \frac{dz}{z} \exp{(t z)} \frac{K_0(\sqrt{z} r)}{K_0(\sqrt{z})} = \int_{\infty}^0 \frac{dx}{x} \: e^{-x t} \frac{K_0(i \sqrt{x} r)}{K_0(i \sqrt{x})}$$
On $C_5$, however, we parametrize by $z=e^{-i \pi} x$ and the integral along $C_5$ becomes
$$\int_{C_5} \frac{dz}{z} \exp{(t z)} \frac{K_0(\sqrt{z} r)}{K_0(\sqrt{z})} = \int_0^{\infty} \frac{dx}{x} \:e^{-x t} \frac{K_0(-i \sqrt{x} r)}{K_0(-i \sqrt{x})}$$
We may now write
$$-1 +\frac{1}{i 2 \pi} \int_0^{\infty} \frac{dx}{x} \: e^{- x t} \left [ \frac{K_0(-i \sqrt{x} r)}{K_0(-i \sqrt{x})} - \frac{K_0(i \sqrt{x} r)}{K_0(i \sqrt{x})} \right ] + \frac{1}{i 2 \pi} \int_{c-i \infty}^{c+i \infty} \frac{ds}{s} \: e^{s t} \frac{K_0( \sqrt{s} r)}{K_0(\sqrt{s})} = 0$$
Therefore, the ILT of $\hat{u}(r,s) = \frac{K_0(\sqrt{s} r)}{s K_0(\sqrt{s})}$ is given by
$$\begin{align}\frac{1}{i 2 \pi} \int_{c-i \infty}^{c+i \infty} \frac{ds}{s} \: e^{s t} \frac{K_0(\sqrt{s} r)}{K_0(\sqrt{s})} &= 1-\frac{1}{i 2 \pi} \int_0^{\infty} \frac{dx}{x} \: e^{- x t} \left [\frac{K_0(-i \sqrt{x} r)}{K_0(-i \sqrt{x})} - \frac{K_0(i \sqrt{x} r)}{K_0(i \sqrt{x})} \right ]\\ &=1- \frac{2}{\pi} \int_{0}^{\infty} \frac{du}{u} \: \,e^{-t u^2} \frac{J_0(u) Y_0(r u)-J_0(r u) Y_0(u)}{J_0^2(u)+Y_0^2(u)} \end{align}$$
where $J_0$ and $Y_0$ are the Bessel functions of the first and second kind, respectively. The last step is a result of the identity
$$K_0(i b) = -\frac{\pi}{2} [Y_0(b)+i J_0(b)]$$
How do I know that the above integral even converges? Using the properties of these Bessels near the origin, I get that the integrand behaves as $2 \log{r}/(u \log^2{u})$ as $u \to 0^+$. Therefore, while there is a singularity at the origin, that singularity is indeed integrable. (This non-trivial result delayed this answer by several hours.)
This is about as far as I can take the integral so far. I will post additional results if I can simplify it.
BONUS
The integral form derived for the solution takes the form of something called an inverse Weber transform (see this PDF reference, for example, on the second page). Essentially, by using a Laplace transform approach, I derived the expression for the inverse Weber transform. This of course is no help in figuring out an analytical expression for the solution, but it provides a further sanity check.
Here are some plots of the temperature $u(r,t)$ for $r \in (1,4]$ and $t \in \{0.01,0.1,0.5,1.0,10\}$:
These were generated in Mathematica directly from the integral expression above.
EDIT
I forgot to mention that in the above, it should be clear that I assumed that $t \gt 0$. When $t \lt 0$, we must close the contour $C$ to the right; $c$ is defined such that the line $C_1$ lies to the right of any poles or branch cuts. Thus, the ILT is zero when $t \lt 0$.
Best Answer
We may start with the formal power series of $\frac{1}{\sqrt{1-x}}$: $$\frac{1}{\sqrt{1-x}}=\sum_{n\geq 0}\frac{1}{4^n}\binom{2n}{n}x^n$$ in order to get $$\frac{1}{\sqrt{1-e^{-\sqrt{t}}}}=\sum_{n\geq 0}\frac{1}{4^n}\binom{2n}{n}e^{-n\sqrt{t}}$$ and $$\mathcal{L}\left(\frac{1}{\sqrt{1-e^{-\sqrt{t}}}}\right)(s) = \sum_{n\geq 0}\frac{1}{4^n}\binom{2n}{n}\int_{0}^{+\infty}2u e^{-(nu+su^2)}\,du $$ which equals $$\sum_{n\geq 0}\frac{1}{4^n}\binom{2n}{n}\left( \frac{1}{s}-\frac{\sqrt{\pi } n e^{\frac{n^2}{4 s}} \text{erfc}\left(\frac{n}{2 \sqrt{s}}\right)}{2 s^{3/2}}\right)$$ with a similar behaviour on $\mathbb{R}^+$ to the one of $s^{-3/4}$.