Exact expression of a trigonometric integral Let $a>2$ be a real number and consider the following integral
$$
I(a)=\int_0^\pi\int_0^\pi \frac{\sin^2(x)\sin^2(y)}{a+\cos(x)+\cos(y)} \mathrm{d}x\,\mathrm{d}y
$$

My question. Does there exist a closed-form expression of $I(a)$?

Some comments. Since $a-2<a+\cos(x)+\cos(y)<a+2$ and $\int_0^\pi \int_0^\pi \sin^2(x)\sin^2(y)\ \mathrm{d}x\, \mathrm{d}y=\frac{\pi^2}{4}$, we have the following bounds
$$
\frac{\pi^2}{4(a+2)} < I(a) < \frac{\pi^2}{4(a-2)},
$$ 
however I didn't manage to find an exact expression for $I(a)$. Any help is welcome!
 A: With CAS help:
$$\int _0^{\pi }\int _0^{\pi }\frac{\sin ^2(x) \sin ^2(y)}{a+\cos (x)+\cos (y)}dydx=\\\mathcal{L}_q\left[\int _0^{\pi }\int _0^{\pi }\mathcal{L}_a^{-1}\left[\frac{\sin ^2(x) \sin ^2(y)}{a+\cos (x)+\cos
   (y)}\right](q)dydx\right](a)=\\\mathcal{L}_q\left[\int_0^{\pi } \left(\int_0^{\pi } e^{-q (\cos (x)+\cos (y))} \sin ^2(x) \sin ^2(y) \, dx\right) \, dy\right](a)=\\\mathcal{L}_q\left[\int_0^{\pi }
   \frac{e^{-q \cos (y)} \pi  I_1(q) \sin ^2(y)}{q} \, dy\right](a)=\\\mathcal{L}_q\left[\frac{\pi ^2 I_1(q){}^2}{q^2}\right](a)=\\\frac{a \pi ^2}{2}-\frac{2}{3} a \pi  E\left(\frac{4}{a^2}\right)-\frac{1}{6}
   a^3 \pi  E\left(\frac{4}{a^2}\right)-\frac{2}{3} a \pi  K\left(\frac{4}{a^2}\right)+\frac{1}{6} a^3 \pi  K\left(\frac{4}{a^2}\right)=\\\frac{\pi ^2 \, _3F_2\left(\frac{1}{2},1,\frac{3}{2};2,3;\frac{4}{a^2}\right)}{4 a}$$
for: $a>2$
$$\frac{1}{6} \pi  \left(-2 \left(a^2-4\right) K\left(\frac{a^2}{4}\right)-2 \left(a^2+4\right) E\left(\frac{a^2}{4}\right)+3 \pi  a\right)$$
for: $a<2$
where: $K$,$E$ gives the elliptic integral of the first kind  and second kind.
Mathematica code:
HoldForm[Integrate[(Sin[x]^2*Sin[y]^2)/(a + Cos[x] + Cos[y]), {x, 0, Pi}, {y, 0, Pi}] == (a \[Pi]^2)/2 - 2/3 a \[Pi] EllipticE[4/a^2] - 1/6 a^3 \[Pi] EllipticE[4/a^2] - 2/3 a \[Pi] EllipticK[4/a^2] + 1/6 a^3 \[Pi] EllipticK[4/a^2] == Pi^2/(4 a)*HypergeometricPFQ[{1/2, 1, 3/2}, {2, 3}, 4/a^2]] // TeXForm
Plot a solution:
f[a_?NumericQ] := NIntegrate[(Sin[x]^2*Sin[y]^2)/(a + Cos[x] + Cos[y]), {x, 0, Pi}, {y,0, Pi}]; 
g[a_] := (a \[Pi]^2)/2 - 2/3 a \[Pi] EllipticE[4/a^2] - 
1/6 a^3 \[Pi] EllipticE[4/a^2] - 2/3 a \[Pi] EllipticK[4/a^2] + 
1/6 a^3 \[Pi] EllipticK[4/a^2]; Plot[{f[a], g[a]}, {a, 2, 20}, 
PlotStyle -> {Red, {Dashed, Black}}, PlotLabels -> {"integral", "Analytic solution"}]


A: Partial answer:
As the integrand is an even function, write $$I(a)=\frac14\int_{-\pi}^\pi\int_{-\pi}^\pi\frac{\sin^2x\sin^2y}{a+\cos x+\cos y}\,dx\,dy$$ and let $w=e^{ix}$. Then for $C_1=\{w:|w|=1\}$, we have \begin{align}\int_{-\pi}^\pi\frac{\sin^2x}{a+\cos x+\cos y}\,dx&=\oint_{C_1}\frac{\frac{(w-1/w)^2}{-4}}{b+\frac{w+1/w}2}\frac{dw}{iw}=\frac i2\oint_{C_1}\frac{(w^2-1)^2}{w^2(w^2+2bw+1)}\,dw\end{align} where $b:=a+\cos y$. Inside $C_1$, there is a double pole at $0$, so letting $f(w)$ be the integrand, \begin{align}\operatorname{Res}(f,0)=\lim_{w\to0}\frac d{dw}(w^2f(w))=-2b.\end{align} Solving the quadratic in the denominator of $f$ gives $w=-b\pm\sqrt{b^2-1}$ but only the positive root lies in $C_1$. As this is a simple pole, \begin{align}\operatorname{Res}(f,-b+\sqrt{b^2-1})=\lim_{w\to-b+\sqrt{b^2-1}}(w+b-\sqrt{b^2-1})f(w)=2\sqrt{b^2-1}.\end{align} Hence \begin{align}\int_{-\pi}^\pi\frac{\sin^2x}{a+\cos x+\cos y}\,dx=\frac i2\cdot2\pi i(2\sqrt{b^2-1}-2b)=-2\pi(\sqrt{b^2-1}-b)\end{align} so that \begin{align}I(a)&=-\frac\pi2\int_{-\pi}^\pi\sin^2y\left(\sqrt{(a+\cos y)^2-1}-a-\cos y\right)\,dy\\&=-\frac\pi2\int_{-\pi}^\pi\sin^2y\sqrt{(a+\cos y)^2-1}\,dy+\frac{a\pi^2}2.\end{align} Now let $z=e^{iy}$. Then for $C_2=\{z:|z|=1\}$, we have \begin{align}\int_{-\pi}^\pi\sin^2y\sqrt{(a+\cos y)^2-1}\,dy&=\oint_{C_2}\left(\frac{z-1/z}{2i}\right)^2\sqrt{\left(a+\frac{z+1/z}2\right)^2-1}\cdot\frac{dz}{iz}\\&=\frac i8\oint_{C_2}\frac{(z^2-1)^2}{z^4}\sqrt{(z^2+2az+1)^2-4z^2}\,dz.\end{align} Evidently, there is a pole of order $4$ at $0$ which lies in $C_2$. However, due to the presence of $\sqrt\cdot$, branch points also need to be considered, and solving the equation $z^2+2az+1=\pm2z$ reveals that only $z=\pm1-a+\sqrt{a^2\mp2a}$ are the two branch points in $C_2$.
A: We use the identities
$$
2\sin \phi \sin \psi=\cos(\phi-\psi)-\cos(\phi+\psi)
$$
and
$$
\cos \phi+\cos \psi=2\cos\left(\frac{\phi+\psi}{2}\right)\cos\left(\frac{\phi-\psi}{2}\right)
$$
to get
$$
I=\int^{\pi}_{0}\int^{\pi}_{0}\frac{\sin^2 x \sin^2 y}{a+\cos x+\cos y}dxdy=
$$
$$
\frac{1}{4}\int^{\pi}_{0}\int^{\pi}_{0}\frac{(\cos(x-y)-\cos(x+y))^2}{a+2\cos\left(\frac{x+y}{2}\right)\cos\left(\frac{x-y}{2}\right)}dxdy
$$
Now we use $\cos \phi=2\cos^2(\phi/2)-1$ and arrive to
$$
I=\int^{\pi}_{0}\int^{\pi}_{0}\frac{(\cos^2\left(\frac{x-y}{2}\right)-\cos^2\left(\frac{x+y}{2}\right))^2}{a+2\cos\left(\frac{x+y}{2}\right)\cos\left(\frac{x-y}{2}\right)}dxdy=
$$
$$
=a^{-1}\int^{\pi}_{0}\int^{\pi}_{0}\frac{(\cos^2 A-\cos^2 B)^2}{1+2/a \cos A\cos B}dAdB.
$$
Then set $t=\cos A$, $s=\cos B$ and use
$$
\frac{d\left(\cos^{(-1)}(x)\right)}{dx}=-\frac{1}{\sqrt{1-x^2}}
$$
to get
$$
I=1/a\int^{1}_{-1}\int^{1}_{-1}\frac{(t^2-s^2)^2}{1+(2/a)ts}\frac{1}{\sqrt{1-t^2}}\frac{1}{\sqrt{1-s^2}}dtds
$$
If $|2/a|\leq  1$ and $|t|,|s|<1$, we expand $(1+(2/a)ts)^{-1}$ into $\sum^{\infty}_{k=0}(-1)^k(2/a)^kt^ks^k$ we can split the two integrals (that of $t$ and that of $s$). We set
$$
I_k:=1/a\int^{1}_{-1}\int^{1}_{-1}(t^2-s^2)^2(-1)^k(2/a)^kt^ks^k\frac{1}{\sqrt{1-t^2}}\frac{1}{\sqrt{1-s^2}}dtds.
$$
Then $I_{2k+1}=0$  and
$$
I_{2k}=\frac{2^{2k}(1/a)^{2k+1}\pi\Gamma\left(\frac{2k+1}{2}\right)\Gamma\left(\frac{2k+3}{2}\right)}{\Gamma\left(k+2\right)\Gamma(k+3)}.
$$
Hence summing
$$
I=I(a)=\sum^{\infty}_{k=0}I_{2k}=\frac{a\pi^2}{2}-\frac{a\pi^2}{2}{}_2F_1\left(\frac{-1}{2},\frac{1}{2};2;\frac{4}{a^2}\right)=
$$
$$
=\frac{a\pi^2}{2}-\frac{a(a^2+4)\pi}{6}E\left(4a^{-2}\right)+\frac{a(a^2-4)\pi}{6}K(4a^{-2}),
$$
where $K(x)$,$E(x)$ are the complete elliptic integrals of the first and second kind. Examples of evaluations can be given if we consider the singular modulus $k_r$. For example with $a=2\sqrt{2}$, we get $k_1=1/2$ and we have
$$
I\left(2\sqrt{2}\right)=\sqrt{2}\pi^2-2\sqrt{2\pi}\Gamma\left(\frac{3}{4}\right)^2-\frac{8}{3}\sqrt{2\pi}\Gamma\left(\frac{5}{4}\right)^2
$$
...etc
see elliptic integrals singular valus
