An elementary way of simplifying a trigonometric triple integral? By stressing my manipulative powers and a bit of help from Mathematica I was able to show that the triple integral
$$\frac1{8\pi}\int_0^{2\pi}\int_0^\pi\int_0^\pi\,(2\cos\theta\cos\varphi-\sin\theta\sin\varphi\cos\psi)^{2k}\sin\theta\sin\varphi \,\mathrm{d}\theta\mathrm{d}\varphi\mathrm{d}\psi $$
where k is a positive integer is equivalent to
$${\left(\frac{2^k}{2k+1}\right)^2}\sum_{j=0}^k\frac{\binom{2j}{j}}{\binom{2k}{k}}$$
by making use of the binomial expansion and then using a hypergeometric identity that degenerates to a ratio of gamma functions.
For some reason I think there was no need for Mathematica to invoke the Gauss hypergeometric function ${}_2 F_1\left(a,b;c;z\right)$ in performing the triple integration, so I'm wondering if there might be a way to show that this identity holds, using nothing more complicated than a gamma function.
 A: Here is a hint to an other approach: Search for the key words "convolution formula" "Legendre polynomial" "addition formula" and you will be en-lighted... 
--- A little explanation and some spherical harmonics ---
@J. M. You are totally right in particular when I was wrong in my suggestion! :/
What I thought about was convolution on the double coset space $X=K\backslash G/K$ of $G= SO(3)$ where $K=SO(2)$. The space $X$ may be identified with $[-1,1]$ and the convolution (as induced from the translation operator on $G$) is then given by
$$f*g(x)= \frac{1}{4\pi}\int_{-1}^1\int_{0}^{2\pi}f(y)g(xy+\sqrt{1-x^2}\sqrt{1-y^2}\cos\theta) d\theta dy.$$
An interesting fact is that $f*g=g*f$ even though $G$ is not Abelian.
Choosing $f(x)=1$ and $g(x)=x^{2k}$ we get
$$\frac{1}{8\pi}\int_{-1}^1\int_{-1}^1\int_{0}^{2\pi}(xy+\sqrt{1-x^2}\sqrt{1-y^2}\cos\theta)^{2k} d\theta dydx=$$
$$\int_{-1}^1f*g(x)dx/2 = \int_{-1}^1 g*f(x)dx/2=\int_{-1}^1x^{2k}dx/2= 1/(2k+1).$$
In our case the integrand is (after changing variables) 
$$(2xy + \sqrt{1-x²}\sqrt{1-y²}\cos\theta)^{2k}=$$
$$=\sum_{j=0}^{2k}{2k \choose j}(xy)^j(xy+\sqrt{1-x²}\sqrt(1-y²)\cos\theta)^{2k-j}$$
Using this we get an other convolution - but I can not see how this can be used directly. :(
Legendre polynomials might be good since they are the characters on the space $X$.
A: Okay ... here we go ...
To save typing, I'm going to use $A$, $B$, and $C$ for $\theta$, $\phi$, and $\psi$. So, our integral starts as this ...
$$
\int\int\int \left(2 \cos A \cos B - \sin A \sin B \cos C\right)^{2k} \sin A \sin B
$$
... and expands to this ...
$$\int\int\int \sum_{j=0}^{2k} {2k \choose j} \; 2^j (-1)^{2k-j} \cos^j A \; \sin^{2k-j+1} A \; \cos^j B \; \sin^{2k-j+1} B \; \cos^{2k-j} C$$
The integral over the sum will be a sum of integrals over the terms; and, since the $A$, $B$, and $C$ factors are independent, each integrated term is a product of integrals. Moreover, as mentioned in my comment on the original question, we're integrating $\cos^{2k-j}C$ over the full period of cosine, so terms with odd powers of $\cos C$ --those corresponding to odd values of $j$-- vanish. So, we take $j=2h$ and for convenience define $m := k-h$, then we write ...
$$\sum_{h=0}^{k} 2^{2h} {2k \choose 2h} \; \left( \int \cos^{2h} A \; \sin^{2m+1} A \right) \left( \int \cos^{2h} B \; \sin^{2m+1} B \right) \left( \int \cos^{2m} C \right)$$
The first two integrals are identical, so we're really looking at the following:
$$\sum_{h=0}^{k} 2^{2h} {2k \choose 2h} \; \left( \int \cos^{2h} A \; \sin^{2m+1} A \right)^2 \left( \int \cos^{2m} C \right)$$
Consider the pieces ...
By the reduction formula for the integral of powers of cosine, we have 
$$\int_{0}^{2\pi} \cos^{2m}C \; dC = \frac{(2m-1)(2m-3)\cdots(1)}{(2m)(2m-2)\cdots(2)}\int_{0}^{2\pi} dC = \frac{2\pi(2m-1)!!}{2^{m} m!} = \frac{2\pi(2m)!}{2^{2m} m!^2}$$
where $x!!$ is the double factorial.
As for the other factor ...
$$\int_{0}^{\pi} \cos^{2h} A \; \sin^{2m+1} A \; dA = \int_{0}^{\pi} \cos^{2h} A \; \left( 1 - \cos^{2}A\right)^m \; \sin A \; dA = \int_{-1}^{1} u^{2h} \left( 1 - u^2 \right)^m \; du$$
$$= \int_{-1}^{1} u^{2h} \sum_{n=0}^{m} {m \choose n} (-1)^n u^{2n} \; du = \sum_{n=0}^{m} (-1)^n {m \choose n} \int_{-1}^{1} u^{2n+2h} \; du = \sum_{n=0}^{m} {m \choose n} \frac{2 (-1)^n}{2n+2h+1} $$
Mathematica evaluates the sum into a form we manipulate into something we can use: 
$$\frac{\Gamma(h+\frac{1}{2}) m!}{\Gamma(m+h+\frac{3}{2})}  = \frac{\sqrt{\pi} (2h-1)!!}{2^h}\frac{2^{k+1}m!}{\sqrt{\pi}(2k+1)!!}= \frac{(2h)!}{2^{2h}h!}\frac{2^{k+1}m!}{(2k+1)}\frac{2^{k}k!}{(2k)!}= \frac{2^{2m+1}m!}{2k+1}\frac{(2h)!}{h!}\frac{k!}{(2k)!}$$
This leads to a simplification bonanza in the formula for the integral ...
$$\sum_{h=0}^{k} \;\; 2^{2h} \;\;\;\; \frac{(2k)!}{(2h)!(2m)!} \;\;\;\;  \frac{2^{4m+2}m!^2}{(2k+1)^2}\frac{(2h)!^2}{h!^2}\frac{k!^2}{(2k)!^2} \;\;\;\; \frac{2\pi(2m)!}{2^{2m} m!^2}$$
$$= 8 \pi \sum_{h=0}^{k} \frac{2^{2(m+h)}}{(2k+1)^2}\frac{(2h)!}{h!^2}\frac{k!^2}{(2k)!}$$
$$= 8 \pi  \frac{2^{2k}}{(2k+1)^2} \frac{k!^2}{(2k)!} \sum_{h=0}^{k}\frac{(2h)!}{h!^2}$$
$$= 8 \pi  \left(\frac{2^{k}}{2k+1}\right)^2 {2k \choose k}^{-1} \;\; \sum_{h=0}^{k}{2h \choose h}$$
... as expected (well, "expected" from your having previously forced it out of Mathematica).
