How do I calculate the triple integral $\int_{0}^{2\pi}\int_0^1 \int_0^1 xy\sqrt{x^2 + y^2 -2xy\cos(\theta)} dx \text{ }dy \text{ } d\theta$?

I have tried a couple substitutions which didn't pan out, and I tried differentiating under the integral sign on $$ I(\theta) = \int_0^1 \int_0^1 xy\sqrt{x^2 + y^2 -2xy\cos(\theta)} dx \text{ }dy $$ to get $$ \frac{dI}{d\theta} = -\sin(\theta) \int_0^1 \int_0^1 \frac{x^2y^2}{\sqrt{x^2 + y^2 -2xy\cos(\theta)}} dx \text{ }dy $$ but nothing has worked so far. Can any of you help me?

By symmetry, we have: $$J=\int_{0}^{2\pi}I(\theta)\,d\theta = 2\int_{0}^{2\pi}\iint_{D}xy\sqrt{x-ye^{i\theta}}\sqrt{x-ye^{-i\theta}}\,dx\,dy\,d\theta $$ where $D$ is the region: $$ D = \{(x,y)\in[0,1]^2 : y\leq x\},$$ so: $$\begin{eqnarray*} J &=& 2\int_{0}^{2\pi}\int_{0}^{1}\int_{0}^{x}xy\sqrt{x-ye^{i\theta}}\sqrt{x-ye^{-i\theta}}\,dy\,dx\,d\theta\\&=&2\int_{0}^{2\pi}\int_{0}^{1}x^4\int_{0}^{1}t\sqrt{1-te^{i\theta}}\sqrt{1-te^{-i\theta}}\,dt\,dx\,d\theta\\&=&\frac{2}{5}\int_{0}^{1}t\int_{0}^{2\pi}\sqrt{1-te^{i\theta}}\sqrt{1-te^{-i\theta}}\,d\theta\,dt.\end{eqnarray*}$$ Since for any $|z|<1$ we have: $$\sqrt{1-z}=\sum_{n\geq 0}\binom{\frac{1}{2}}{n}(-1)^n z^n $$ it happens that: $$ \int_{0}^{2\pi}\sqrt{1-te^{i\theta}}\sqrt{1-te^{-i\theta}}\,d\theta =2\pi\sum_{n\geq 0}\binom{\frac{1}{2}}{n}^2 t^{2n} $$ so: $$ J = \frac{4\pi}{5}\int_{0}^{1}\sum_{n\geq 0}\binom{\frac{1}{2}}{n}^2 t^{2n+1}\,dt = \frac{4\pi}{5}\sum_{n\geq 0}\binom{\frac{1}{2}}{n}^2\frac{1}{2n+2}=\frac{4\pi}{5}\cdot\frac{16}{9\pi}=\color{red}{\frac{64}{45}}.$$ To prove the last identity, notice that: $$\sum_{n\geq 0}\binom{\frac{1}{2}}{n}^2\frac{1}{2n+2}=\sum_{n\geq 0}\frac{\frac{1}{16^n}\binom{2n}{n}^2}{(2n-1)^2(2n+2)}$$ can be computed from: $$ \sum_{n\geq 0}\frac{\frac{1}{4^n}\binom{2n}{n}}{(2n-1)^2}x^{2n}=\sqrt{1-x^2}+x\arcsin x,$$ $$ \sum_{n\geq 0}\frac{\frac{1}{4^n}\binom{2n}{n}}{(2n+2)}x^{2n}=\frac{1-\sqrt{1-x^2}}{x^2}$$ through a contour integral, as already guessed by the OP.

This can be tackled also in an alternative and very elegant way.

If we consider the extremely special function:

$$ B(\lambda)=\sum_{n\geq 0}\left(\frac{\binom{2n}{n}}{(2n-1)4^n}\right)^2\lambda^{2n}=\frac{1}{2\pi}\int_{0}^{2\pi}\sqrt{1+\lambda^2-2\lambda\cos\theta}\,d\theta=\phantom{}_2 F_1\left(-\frac{1}{2},-\frac{1}{2};1;\lambda^2\right)$$ it is easy to check from the power series representation that $B$ satisfies the ODE: $$ \lambda B = (\lambda^2+1) B' + (\lambda-\lambda^3) B'' $$ but from that differential equation it is also easy to check, by integration by parts, that: $$ \int_{0}^{1} \lambda\, B(\lambda)\,d\lambda = \frac{16}{9\pi} $$ since $B(1)=\frac{4}{\pi}=2\cdot B'(1)$. The previous integral is everything we need to solve our problem.