I going to manipulate the integral into a form that I can analyze the method of stationary phase. Let $k=\cos{\theta}$ and the integral becomes
$$\begin{align}I(x) &= \int_0^{\pi} d\theta \, \sin^2{\theta} \, e^{i x \cos{\theta}}\\ &= \int_0^{\pi} d\theta \, e^{i x \cos{\theta}} - \int_0^{\pi} d\theta \, \cos^2{\theta} \, e^{i x \cos{\theta}}\\ &= \int_0^{\pi} d\theta \, e^{i x \cos{\theta}} + \frac{d^2}{dx^2} \int_0^{\pi} d\theta \, e^{i x \cos{\theta}}\end{align}$$
Now, note that
$$\int_0^{\pi} d\theta \, e^{i x \cos{\theta}} = \int_0^{\pi/2} d\theta \, e^{i x \cos{\theta}} + \int_{\pi/2}^{\pi} d\theta \, e^{i x \cos{\theta}} = 2 \Re{\left [ \int_0^{\pi/2} d\theta \, e^{i x \cos{\theta}}\right]}$$
Now, we may apply stationary phase. The stationary point of the integrand is at $\theta=0$; there, we may approximate the argument of the exponential by its Taylor expansion. Further, because of the oscillatory cancllations, we may simply draw the upper limit of the integral to infinity to first order:
$$\int_0^{\pi/2} d\theta \, e^{i x \cos{\theta}} \sim e^{i x} \int_0^{\infty} d\theta \, e^{-(i x/2) \theta^2} = \frac12 e^{i(x-\pi/4)} \sqrt{\frac{2 \pi}{x}} \quad (x\to\infty)$$
Therefore
$$\int_0^{\pi} d\theta \, e^{i x \cos{\theta}} \sim \sqrt{\frac{2 \pi}{x}} \cos{\left ( x-\frac{\pi}{4}\right)}\quad (x\to\infty)$$
To get the asymptotic behavior of $I(x)$, it looks like we need to take the second derivative of the above result. It turns out that
$$\frac{d^2}{dx^2} \left[x^{-1/2} \cos{\left ( x-\frac{\pi}{4}\right)}\right ]= \left (\frac{3}{4} x^{-5/2} - x^{-1/2}\right ) \cos{\left ( x-\frac{\pi}{4}\right)} + x^{-3/2}\sin{\left ( x-\frac{\pi}{4}\right)} $$
Note that the $x^{-1/2}$ piece drops out when combined with the original integral, and the $x^{-5/2}$ piece is subdominant. Thus, the leading asymptotic behavior of the integral $I(x)$ is finally
$$\int_{-1}^1 dk \, \sqrt{1-k^2} \, e^{i k x} \sim \sqrt{2 \pi} x^{-3/2} \sin{\left ( x-\frac{\pi}{4}\right)}\quad (x\to\infty)$$
Here is a plot illustrating this behavior against the exact value of the integral (the lower plot is the asymptotic result derived here):
ADDENDUM
One objection to the above derivation might be that I neglected the next term in the asymptotic expansion of the first integral of $I(x)$, which we know goes as $A x^{-3/2} \sin{(x-\pi/4)}$. How could I ignore that? It turns out $I$ is the sum of this term and its second derivative, and the only term that remains $O(x^{-3/2})$ vanishes because it is a simple sine term, the sum of itself and its second derivative vanishing. So the derivation of the asymptotic approximation above remains valid.
Using the saddle point method with $a=0$, $b=+\infty$, $t_0=\mathrm{e}^{-1}$, $z=\mathrm{e}^n$, $p(t)=t\log t$ and $q(t)=1$, we find
\begin{multline*}
\int_0^{ + \infty } {x^{ - x} \mathrm{e}^{nx} \mathrm{d}x} = \mathrm{e}^n \int_0^{ + \infty } {\mathrm{e}^{ - \mathrm{e}^n t\log t} \mathrm{d}t} \\ \sim \exp \left( {\tfrac{{n - 1}}{2} +\mathrm{e}^{n - 1} } \right)\sqrt {2\pi } \left( {1 + \sum\limits_{k = 1}^\infty {\sqrt {\frac{{2\mathrm{e}}}{\pi }} \Gamma \left( {k + \tfrac{1}{2}} \right)\frac{{b_{2k} }}{{\mathrm{e}^{nk} }}} } \right),
\end{multline*}
where the coefficients are given as complex residues:
$$
b_{2k} = \frac{1}{2}\mathop{\operatorname{Res}}\limits_{t = \mathrm{e}^{ - 1} }\left[ {(t\log t + \mathrm{e}^{ - 1} )^{ - k - 1/2} } \right] =
\frac{1}{2}\mathrm{e}^{k - 1/2} \mathop{\operatorname{Res}}\limits_{s = 0} \left[ {((1 + s)\log (1 + s) - s)^{ - k - 1/2} } \right].
$$
For example,
$$
b_2 = - \frac{1}{{12}}\sqrt {\frac{\mathrm{e}}{2}} ,\quad b_4 = - \frac{{23\mathrm{e}}}{{864}}\sqrt {\frac{\mathrm{e}}{2}} ,
$$
whence
$$
\int_0^{ + \infty } {x^{ - x} \mathrm{e}^{nx} \mathrm{d}x} \sim \exp \left( {\tfrac{{n - 1}}{2} + \mathrm{e}^{n - 1} } \right)\sqrt {2\pi } \left( {1 - \frac{1}{{24}}\frac{1}{{\mathrm{e}^{n - 1} }} - \frac{{23}}{{1152}}\frac{1}{{\mathrm{e}^{2(n - 1)} }} + \cdots } \right).
$$
Best Answer
Let $g(z) = f(z) - \sqrt{2}$, and consider the substitution
$$1-s = \frac{x}{\sqrt{x^2+z^2}}, \qquad 1-t = \frac{y}{\sqrt{y^2+z^2}}, \qquad w = 1-\frac{1}{\sqrt{z^2+1}}.$$
Then from the computation
$$\mathrm{d}x=-\frac{z}{s^{3/2}(2-s)^{3/2}} \mathrm{d}s, \qquad \mathrm{d}y=-\frac{z}{t^{3/2}(2-t)^{3/2}} \mathrm{d}t, $$
we obtain the following integral representation:
\begin{align*} g = g(z) &= z^2 \int_{w}^{1} \int_{w}^{1} \frac{\sqrt{s+t-st} + \sqrt{2-(s+t-st)} - \sqrt{2}}{(st)^{3/2} (2-s)^{3/2}(2-t)^{3/2}} \, \mathrm{d}s\mathrm{d}t \\ &= 2 z^2 \int_{w}^{1} \int_{t}^{1} \frac{\sqrt{s+t-st} + \sqrt{2-(s+t-st)} - \sqrt{2}}{(st)^{3/2} (2-s)^{3/2}(2-t)^{3/2}} \, \mathrm{d}s\mathrm{d}t. \end{align*}
Now by noting that $w \sim \frac{z^2}{2}$ as $z \to 0$, we show that $g \sim c\sqrt{w}\log w$ as $w \to 0^+$ for some constant $c \neq 0$. Indeed,
\begin{align*} &\lim_{w \to 0^+} \frac{g}{\sqrt{w}\log w} \\ &= \lim_{w \to 0^+} \frac{4}{w^{-1/2}\log w} \int_{w}^{1} \int_{t}^{1} \frac{\sqrt{s+t-st} + \sqrt{2-(s+t-st)} - \sqrt{2}}{(st)^{3/2} (2-s)^{3/2}(2-t)^{3/2}} \, \mathrm{d}s\mathrm{d}t \\ &= \lim_{w \to 0^+} \frac{8}{w^{-3/2}\log w} \int_{w}^{1} \frac{\sqrt{s+w-sw} + \sqrt{2-(s+w-sw)} - \sqrt{2}}{(sw)^{3/2} (2-s)^{3/2}(2-w)^{3/2}} \, \mathrm{d}s \\ &= \lim_{w \to 0^+} \frac{8}{\log w} \int_{1}^{1/w} \frac{\sqrt{w(r+1-wr)} + \sqrt{2-w(r+1-wr)} - \sqrt{2}}{\sqrt{w} r^{3/2} (2-wr)^{3/2}} \, \mathrm{d}r, \end{align*}
where the L'Hospital's Rule is applied in the second step and the substitution $s=wr$ is utilized in the last step. Now it is not hard to show that the last limit is $-1$, and therefore,
$$ g(z) \sim -\sqrt{w}\log w \sim -\sqrt{2}z\log z \qquad \text{as} \qquad z \to 0^+. $$