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