Multivariate Elliptic Integral I want to find a "closed form" solution to the following integral:
$$I(a)=\int_{-1}^1\int_{-1}^1\int_{-1}^1\sqrt{\frac{a+ xy+yz+xz }{(1-x^2)(1-y^2)(1-z^2)}}\, \text{d}x\,\text{d}y\,\text{d}z$$
where $a>1$ to keep the integrand real valued.
I hope that the symmetry of the integrand might help.
I have tried integrating one variable after the other, but this involves Integrals of elliptic integrals that I am not familiar with:
Another form of the integral is:
$$I(a)=\int_{0}^\pi\int_0^\pi\int_0^\pi\sqrt{ a+ \cos(x)\cos(y)+\cos(y)\cos(z)+\cos(x)\cos(z) } \,\text{d}x\,\text{d}y\,\text{d}z$$
Evaluating the first integral, Mathematica gives:
$I(a)=\int_0^\pi\int_0^\pi 2\sqrt{ a+\cos(y)+\cos(y)\cos(z)+ \cos(z) }E\left[\frac{2\cos(y)+2\cos(z) }{ a+\cos(y)+\cos(y)\cos(z)+ \cos(z) }\right]  \,\text{d}y\,\text{d}z$
Where $E$ is the Elliptic Integral of second kind. I dont know how to proceed.
Is there maybe a nice coordinate system I could transform into? Maybe I have to try a completely different approach?
Any insight would be helpful.
Edit:
I have made some progress by using the FCC Lattice Green's function.
It is defined by:
$$G(a)=\frac{1}{\pi^3}\int_{0}^\pi\int_0^\pi\int_0^\pi\frac{1}{ 1 - \frac{a}{3}( \cos(x)\cos(y)+\cos(y)\cos(z)+\cos(x)\cos(z)) } \,\text{d}x\,\text{d}y\,\text{d}z$$
and has a known closed form solution in terms of Complete Elliptic Functions that can be found here (equation 2.20)
I will add more info soon.
 A: I was able to reduce the integral to a one-dimensional integral. I felt this was to long for an edit to my question.
I will consider a slightly different Integral than I originally stated: 
$$J(a)=\frac{1}{\pi^3}\int_{0}^\pi\int_0^\pi\int_0^\pi\sqrt{ 1-\frac{a}{3} (\cos(x)\cos(y)+\cos(y)\cos(z)+\cos(x)\cos(z)) } \,\text{d}x\,\text{d}y\,\text{d}z$$
but it can easily be transformed into the original by considering $I(a) = \sqrt{a}J(-\frac{3}{a})$
I assume $\lvert a\rvert < 1$, but this restriction can probably be lifted to a $-3 < a < 1$.
Using the series expansion of the integrand and by exchanging integrals and summation we can rewrite J(a) as
$$J(a) = \sum_{n=0}^{\infty} \binom{\frac{1}{2}}{n} c_n (-a)^n $$
where the coefficients $c_n$ are defined by
$c_n = \frac{1}{\pi^3}\int_{0}^\pi\int_0^\pi\int_0^\pi \frac{1}{3^n}( \cos(x)\cos(y)+\cos(y)\cos(z)+\cos(x)\cos(z))^n$
Both functions $\sum_n^{\infty} \binom{\frac{1}{2}}{n}(-a)^nx^n = \sqrt{1-ax}$ and $\sum_n^{\infty} c_n x^n =W_1(x)$ are known.
$$W_1(x) =\frac{3\sqrt{3}}{(3+x)^{\frac{3}{2}}}(2-\sqrt{1-x})[ \frac{2}{\pi} K(\frac{1}{2}-\frac{2\sqrt{3}x}{(3+x)^{\frac{3}{2}}} -\frac{\sqrt{3}(3-x)\sqrt{1-x}}{(3+x)^{\frac{3}{2}}})]^2$$
is the Lattice Greens function of the FCC lattice. $K(k)$ is the complete Elliptic Integral of the first kind. A link to a paper with this closed form can be found here. The formula I used is Nr. (2.26).
The Integral is then the Hadamard product of both Series $W_1(x)$ and $\sqrt{1-ax}$ at x = 1:
$$J(a) = \frac{1}{2\pi}\int_0^{2\pi}\sqrt{1-ae^{-it}}W_1(e^{it}) \text{d}t$$
This can be rewritten as a contour integral:
$$J(a) = \frac{1}{2\pi i}\oint_{\lvert z \rvert=1}\frac{1}{z}\sqrt{1-\frac{a}{z}}W_1(z)$$
Now $W_1(z)$ has its Branch points at $-3,1, \infty$ so they wont bother us. But $\sqrt{1-\frac{a}{z}}$ has its branch points at $0$ and $a$. So we can modify the circular contour to a dogbone contour around the branch cut from 0 to a. With a circle around 0 with infinitesimal radius $\epsilon$, a circle around $a$ with infinitesimal radius $r$ and the back and forth integrals on the Real line which are complex conjugates of each other.
$$\oint_{\lvert z\rvert=1} = \oint_{\lvert z \rvert=\epsilon}  +\oint_{z = a + ry, \lvert y\rvert=1 } - 2i \text{Im}[\int_\epsilon^{a-r} ] $$
The Integral around $a$ just gives $0$ for $r\xrightarrow{}0$. The Integral around $0$ is trickier since it diverges for $\epsilon = 0$ . This can be solved by first integrating for infinitesimal $\epsilon$, then adding the integrals together, and then letting $\epsilon$ got to 0. This way all divergences cancel. In the neighbourhood of $z=0$ the Lattice Greens function is just $1$. Therefore the integral gives
$$\frac{1}{2\pi i} \oint_{\lvert z \rvert=\epsilon}\frac{1}{z}\sqrt{1-\frac{a}{z}}W_1(z)\text{d}z \approx \frac{1}{2\pi i} \oint_{\lvert z \rvert =\epsilon}\frac{a}{z^{3/2}} \text{d}z + \mathcal{O}(\epsilon)= \frac{2}{\pi}\sqrt{\frac{a}{\epsilon}}=  \frac{2}{\pi}+\frac{1}{\pi}\int_\epsilon^a\frac{\sqrt{a}}{z^\frac{3}{2}}\text{d}z$$
After adding this to the remaining integral we arrive at the 1D Integral form:
$$J(a)= \frac{2}{\pi}+\frac{1}{\pi}\int_0^a\frac{\sqrt{az}-\sqrt{az-z^2}W_1(z)}{z^2}\text{d}z = \frac{2}{\pi}+\frac{1}{\pi}\int_0^1\frac{\sqrt{ z}-\sqrt{ z-z^2}W_1(az)}{z^2}\text{d}z $$
I have checked the equality numerically for values of $-3<a<1$.
Unfortunately I am stuck here. This is quite a lot nicer than the original expression but sadly still far away from any final solution. Maybe after a nice substitution one can use contour integration again. But the Branch structure of $W_1(x)$ might cause a lot of Problems. Maybe someone can handle the Series expansion of the Integral. The closed form of $c_n$ can be found here in equation (13) and (17).
Edit:
An even nicer form of the resulting integral is:
$$J(a)=1-\frac{1}{\pi}\int_0^1\frac{\sqrt{ z-z^2}}{z^2}(W_1(az)-1)\text{d}z$$
If one applies the Euler substitution $z-> \frac{1}{1+t^2}$ one can get rid of the square root
$$J(a)=1-\frac{2}{\pi}\int_0^\infty \frac{t^2}{1+t^2}(W_1(\frac{a}{1+t^2})-1)\text{d}t$$
