Asymptotic behaviour of an integral whose integrand cannot be expanded How does one evaluate how a limit is approached in the following case? 
Consider the function
$$
f(z) = \int_{0}^1 d x \int_{0}^1 d y \left( \sqrt{1 - \frac{x y}{\sqrt{x^2 + z^2}\sqrt{y^2 + z^2}}} + \sqrt{1 + \frac{x y}{\sqrt{x^2 + z^2}\sqrt{y^2 + z^2}}} \right).
$$
This has limiting values 
$$
\lim_{z \to 0}f(z) = \sqrt{2} \quad \text{and} \quad \lim_{z \to \infty}f(z) = 2.
$$
and interpolates monotonically between them for intermediate values. I am interested in establishing the expansions of $f(z)$ about these limits. Specifically I want to know the leading order behaviour of 
$$
f(z) - \sqrt{2} \quad \text{as} \quad z \to 0, \qquad \text{and} \qquad f(z) - 2 \quad \text{as} \quad z \to \infty.
$$
In the case of $z \to \infty$ the expansion about $z = \infty$ is easily established by Taylor expansion:
$$
f(z) = \int_{0}^1 d x \int_{0}^1 d y \left( 2 - \frac{x^2y^2}{4z^4} + \mathrm{O}(z^{-6}) \right) = 2 - \frac{1}{36 z^4} + \mathrm{O}(z^{-6}).
$$
However the same trick fails in the limit $z \to 0$
$$
f(z) = \int_{0}^1 d x \int_{0}^1 d y \left( \sqrt{2} + \frac{z\sqrt{x^2+y^2}}{\sqrt{2} x y} + \mathrm{O}(z^2) \right)
$$
as the second term in the integral does not converge. This makes sense as the integrand is not Taylor expandable on the lines $x=0,y=0$ when $z=0$. 
How does one determine how the limit is approached in this case?
Numerically it appears that $(f(z) - f(0)) \sim z \log z$ as $z \to 0$ but I have been unable to show this formally
 A: 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^+. $$
