Find the leading singular behavior of a parameter-dependent double integral The integral
$$S(\eta)=\int^1_{-1}\int^1_{-1} \sqrt{x^2+\eta y^2 +\eta^2 +y^4}dx dy$$
is not real-analytic at $\eta\to 0^+$. The question is to find out its leading singular behavior at $\eta\to 0^+$. Numerical results seem to suggest the behavior $S(\eta) \propto \eta^{5/2}$, i.e. $S'''(\eta)$ diverges as $1/\sqrt{\eta}$ when $\eta\to 0^+$.
 A: There is no problem for the inner integral
$$\int^1_{-1} \sqrt{x^2+\eta y^2 +\eta^2 +y^4}\,dx =$$
$$\sqrt{\eta ^2+y^4+\eta  y^2+1}+\frac 12 (\eta ^2+y^4+\eta  y^2)\log \left(\frac{\sqrt{\eta ^2+y^4+\eta  y^2+1}+1}{\sqrt{\eta ^2+y^4+\eta 
   y^2+1}-1}\right)$$ For the outer antiderivative no problem except a bunch of elliptic integrals of all kinds. The resulting integral is just a monster.
The only thing which is simple is to evaluate the value of the definite integral for $\eta=0$; it is
$$k=\frac{2}{15} \left(-8 \,
   _2F_1\left(-\frac{3}{4},\frac{1}{2};\frac{1}{4};-1\right)+15 \sqrt{2}+3 \sinh
   ^{-1}(1)+\frac{6 \sqrt{2 \pi } \Gamma \left(\frac{5}{4}\right)}{\Gamma
   \left(\frac{7}{4}\right)}\right)\sim 2.66132$$ But is seems to me that, from numerical calculations, using
$$S(\eta)=k +a\, \eta^b$$  it seems to me that $b\sim  \frac 54$
To allow a double check, I give below a few values
$$\left(
\begin{array}{cc}
\eta & S(\eta)- k \\
 0.0 & 0.00000 \\
 0.1 & 0.13385 \\
 0.2 & 0.32357 \\
 0.3 & 0.55336 \\
 0.4 & 0.81345 \\
 0.5 & 1.09701 \\
 0.6 & 1.39905 \\
 0.7 & 1.71580 \\
 0.8 & 2.04440 \\
 0.9 & 2.38266 \\
 1.0 & 2.72884
\end{array}
\right)$$
I have been unable to expand $S(\eta)$ as a series around $\eta=0$.
