Well, I can reduce the problem to a 1-dimensional integral for you. Perhaps someone from the community can take it all the way. I use some basic results of stochastic geometry to reduce to the 1-D case. I'm going to ramble a bit because I tried several things and I think they are each interesting in themselves.
First note: The average of the cosine as described in THIS post, which is related to the current one, is wrong. The solution given in the current post by user31373 is the correct one. Installero was missing a $\sin\theta_1\sin\theta_2$ in the integral. Here is what it should look like:
The mean of the cosine of the angle between random rays within the cone of aperture $2x$ is given by (for reference, the dot product between two unit vectors given in spherical coordinates is provided in THIS post, notably the answer given by Dan Shved.): $$\frac{\int_0^{2\pi}\int_0^{2\pi}\int_0^{x}\int_0^{x}\left (\cos\theta_1\cos\theta_2+\sin\theta_1\sin\theta_2\cos(\phi_1-\phi_2) \right ) \sin\theta_1\sin\theta_2d\theta_2d\theta_2d\phi_1d\phi_2}{\int_0^{2\pi}\int_0^{2\pi}\int_0^{x}\int_0^{x}\sin\theta_1\sin\theta_2d\theta_2d\theta_2d\phi_1d\phi_2}\\
=\frac{\pi^2\sin^4x}{16\pi^2\sin^4(x/2)}\\
=\cos^4(x/2)$$
Again, the solution originally provided by Installero, namely $\frac{\sin^2x}{x^2}$, is incorrect. Here, the integral in the denominator is necessary to normalize the distribution function given in the integral in the numerator.
Let us first consider the cosine case. As I recently described in THIS post, since we are interested in the mean of a function of a single quantity, in this case the cosine of the angle between random rays in the cone, we can instead consider a one-dimensional integral (if nothing else, a one dimensional integral numerically integrates much faster than a 3-dimensional one): $$\int_0^{2x}\cos(\theta)p(x,\theta)d\theta,$$ where $p(x,\theta)$ represents the appropriate density function that returns the solution $\cos^4(x/2)$; $p$ represents the probability density of the angle $\theta$ between randomly chosen rays within the cone of aperture $2x$. As a probability density, it must satisfy the following equation: $$\int_0^{2x}p(x,\theta)d\theta=1 $$We now proceed to find this probability density function $p$ and then use it to find (or at least try to find) the mean of any other function of $\theta$ over the cone, including sine.
First (failed) attempt - integral equations:
We are going to try to leverage the known solution for the cosine case in order to find other solutions. We have:$$\int_0^{2x}\cos(\theta)p(x,\theta)d\theta=\cos^4(x/2).$$It is worth noting that this looks very much like a Volterra integral equation of the first kind except we are looking to solve for the kernel. We now proceed to change this integral equation into a differential equation that Mathematica can solve. First, we must assume that $p$ is separable such that $$p(x,\theta)=a(x)b(\theta)$$ (unfortunately, as seen later, it is not actually separable, at least not obviously). Making careful use of the derivative rules given on THIS Wikipedia page, we can change the above integrals into a system of coupled differential equations (and then decouple them). Let $f(x)\equiv\cos^4(x/2)$. From the normalization equation we find $$b(2x)=\frac{-a'(x)}{2a(x)^2},$$ and from the other equation we find $$a'(x)\frac{f(x)}{a(x)}+a(x)\left[ 2\cos(2x)b(2x) \right]=f'(x),$$which gives $$f'(x)a(x)=a'(x)[f(x)-\cos(2x)].$$ Mathematica solves this equation: $$a(x)=\frac{(5 + 7 \cos x)^{1/21}}{\sin^{2/3}(x/2)}.$$ As a side note, I think this solution may not be good beyond the hemisphere case ($x=\pi/2$) and there may be a clever way to simplify it but I couldn't find one. Then, from the normalization differential equation we have:$$b(2x)=\frac{2 \cos^3(x/2)}{(5 + 7 \cos[x])^{22/21} \sin^{1/3}(x/2)}\\
\Rightarrow b(\theta)=\frac{2 \cos^3(\theta/4)}{(5 + 7 \cos[\theta/2])^{22/21} \sin^{1/3}(\theta/4)}.$$ Again, perhaps there is a clever simplification. We now have the full $p(x,\theta)$. First, we can check that it is correct by numerical integration and comparison with the known solution:
$a(x)\int_0^{2x}\cos(\theta) b(\theta)d\theta= \cos^4(x/2)?$ Lets plot and check:
.
Okay looks good. Now we try with the sine. Mathematica gives an "analytic" solution: $$a(x)\int_0^{2x}\sin(\theta) b(\theta)d\theta\\
=\frac{\tan(x/2)}{70} (-175 \cos^2[x/2] + 35 \cos^2[x/2] \cos[x] +
54\times 2^{17/42}\times 3^{20/21}\times F_1[5/6, -(20/21), 1/2, 11/6, (7/6) \sin^2(x/2),
\sin^2(x/2)]\times \sqrt{1 + \cos x} \times [5 + 7 \cos x]^{1/21} -
12\times 2^{17/42}\times 3^{20/21}\times
F_1[5/6, 1/21, 1/2, 11/6, (7/6) \sin^2(x/2), \sin^2(x/2)]\times \sqrt{
1 + \cos x}\times [5 + 7 \cos x]^{1/21}),$$
where $F_1(a,b_1,b_2,c;x,y)$ is the Appell $F_1$ hypergeometric series. This is what it looks like:
Unfortunately, this is wrong. The assumption that p is separable is incorrect. As seen by Monte Carlo methods, the average of the sine should look something like this (similar but not the same):
Second (partially successful) attempt - stochastic geometry:
The first attempt to find $p(x,\theta)$ was unsuccessful because the assumption of separability (at least that particular form of separability) was bad. Instead we will follow the basic method laid out by García-Pelayo (Ricardo García-Pelayo 2005 J. Phys. A: Math. Gen. 38 3475). The key here is to note that on a unit spherical cap, the geodesic distance between two points is equal to the angle between the vectors directed to those points from the center of the sphere. In this case we are trying to find the probability density for that angle and therefore an equivalent problem is to find the probability density of the geodesic distance between two points on the unit cap.
The trick now is to understand how to parameterize a tilted unit cap (for reasons which will become clear later). I did it in cylindrical coordinates but spherical coordinates would have been fine too. First, we parameterize the un-tilted cap and then rotate it about the y-axis to get the tilted case:
Untilted:
$$x=\sqrt{1-T^2}\cos t,y=\sqrt{1-T^2}\sin t,z=T,0\leq t< 2\pi,\cos m \leq T \leq 1,$$ where $m$ is denotes the half-aperture of the cone (that is, $x$ as defined at the beginning of this post). Now we rotate around the y-axis by angle $\vartheta$:
Tilted:
$$x=\sqrt{1-T^2}\cos t \cos\vartheta+T\sin\vartheta,y=\sqrt{1-T^2}\sin t,z=T\cos\vartheta-\sqrt{1-T^2}\cos t \sin\vartheta,\\0\leq t< 2\pi,\cos m \leq T \leq 1,$$
Now, to describe it in the original parameter system we must also have
$$x=\sqrt{1-\zeta^2}\cos t',y=\sqrt{1-\zeta^2}\sin t',z=\zeta,\\-a\leq t'<a,\cos (m+\vartheta) \leq \zeta \leq 1.$$
The trick now is to find that $a$. This is accomplished by first finding T in terms of $\zeta$, $t'$, $m$, and $\vartheta$ and then using the inequality $\cos m \leq T \leq 1$. We find $$T=\frac{\sqrt{1-\zeta^2}\cos t' + \zeta \cot \vartheta}{\cos\vartheta\cot\vartheta+\sin\vartheta}. $$
Careful application of the inequality gives $$a=\text{Re}\left ( \arccos \left [ \frac{(\cos m -\zeta\cos\vartheta)\csc\vartheta}{\sqrt{1-\zeta^2}} \right ] \right ). $$
Lets check ($m=\pi/2$ and $\vartheta=\pi/4$):
Looks good. Now we apply the method given by García-Pelayo. In equation (1) he notes that $$p(l)=\frac{\int_Dd\vec{y}\int_Dd\vec{x}\delta(|\vec{x}-\vec{y}|-l)}{\int_Dd\vec{y}\int_Dd\vec{x}}, $$where $p(l)$ is the probability density for the distance between two points in the domain $D$. In equation (2) he notes that the numerator can be rewritten as $$\int_Dd\vec{x}O(\vec{x},l),$$ where in his case $O(\vec{x},l)$ represents the length of the overlap of a circle of radius $l$ with the domain $D$. In our case, $O(\vec{x},l)$ represents the overlap of of a thin spherical segment of geodesic-distance-from-the-north-pole $l$ with the spherical cap which has been tilted by angle $\vartheta \leq m \leq \pi/2$. These are nothing more than the latitude lines in the left picture above. We note again that $l=\theta$ on the unit cap (this is the spherical coordinate $\theta$ not the tilt angle $\vartheta$). We also note that $\int_Dd\vec{x}=2\pi\int_0^m\sin(x) dx$, which is the surface area of the cap. We now have two cases, $\theta<m$ and $\theta>m$, each of which can be split up into two more case. In the $\theta<m$ case we have $x+\theta<m$ and $x+\theta>m$. In the $x+\theta<m$ case we have $O(\vec{x},\theta)=2\pi\sin\theta$. For $x+\theta>m$ we have $O(\vec{x},\theta)=2\sin\theta a(x,\theta)$, where $a$ is the one from the parameterized cap, above. Noting that $\zeta=\cos\theta$ and $\vartheta=x$ we find $$a=\arccos([\cos m -\cos\theta\cos x]\csc x\csc \theta).$$
In the $\theta>m$ case we have $x+m<\theta$ and $x+m>\theta$. For $x+m<\theta$ we have $O(\vec{x},\theta)=0$. For $x+m>\theta$ we have $O(\vec{x},\theta)=2\sin\theta a(x,\theta)$. Putting this all together gives $p(m,\theta):$
$$p(m,\theta) = 4\pi\sin(\theta) f(m) \times\\ \left\{
\begin{array}{lr}
\pi\int_{0}^{m-\theta}\sin(x)dx + \int_{m-\theta}^{m}\arccos([\cos m -\cos\theta\cos x]\csc x\csc \theta)\sin(x)dx & :0 \leq \theta \leq m\\
\int_{\theta-m}^{m}\arccos([\cos m -\cos\theta\cos x]\csc x\csc \theta)\sin(x)dx & : m \leq \theta \leq 2m
\end{array}
\right. $$
Here, $f(m)=(2\pi(1-\cos m))^{-2}$ is the normalization factor represented by the denominator in García-Pelayo's equation (1). The integral with the inverse cosine is a serious pain. Mathematica spits out a solution but it's horrible. Here is my attempt at simplification: $$p(m,\theta) =\frac{\cos(\theta/2)\sin(\theta/2)\csc(m/2)^4}{4\pi}\left( \pi \\ -
2 \arctan\left [
1/2 \cos(\theta/2) (\cos[m] - \cos[\theta]) \sqrt{(-\cos[2 m] + \cos[\theta])/(
1 + \cos[\theta])} \csc(m - \theta/2) \csc([m + \theta]/2)^2 \right ] \\ +
\arctan\left [\frac{(\cos[m] + \cos[\theta]) \sqrt{(-\cos[2 m] + \cos[\theta])/(1 + \cos[\theta]}
\cot(\theta/2))}{\cos(m) + \cos(\theta) + \sin(m)^2} \right ]\\ -
4 \arccos[\cot(m) \tan(\theta/2)] \cos[m] \right ).$$ Clearly it is still quite ugly. I am willing to bet that the two arctans and the arccos can be combined into one but I have had trouble doing that and still keeping things on the correct branch. Anyway, lets check it against Monte Carlo. The following is the probability density for the angle $\theta$ between two randomly chosen vectors within a cone of half-aperture $m=3 \pi/8$:
Looks good. Now the question is, does $\int_0^{2x}\cos(\theta)p(x,\theta)d\theta=\cos(x/2)^4$? The answer is yes. Mathematica can do the integral if you are careful about limits and things. However, Mathematica cannot do the sine case. I have the feeling that if one could come up with a nice and simplified form of $p(m,\theta)$ that Mathematica might have a chance of succeeding.
So here is my challenge to the community: Find a better expression for $p(m,\theta)$ and try to compute $\int_0^{2x}\sin(\theta)p(x,\theta)d\theta$ with it. One alternate route might be to leverage the cosine solution now that we know $p$, that is, do some trick with derivatives under the integral sign and get a nice differential equation, but I couldn't find a good way. Also, I assume this result for $p(m,\theta)$ is in the literature from 50 years ago somewhere and they probably have a nicer expression. Find it!
Best Answer
Start from
$$\left\langle\Psi_{2p1} | \Psi_{2px}\right\rangle = \sqrt{1 \over{2^{11}\pi^2 a^{10}}}\int\limits_{0}^\infty \int\limits_{0}^\pi \int\limits_{0}^{2\pi}r^4e^{-r \over a} \sin^3\theta \cos\phi e^{i\phi}~ \mathrm d\phi~\mathrm d\theta~\mathrm dr $$ (which, note, is not the same as the expression in the question statement; thanks @joriki.)
Expand the $e^{i\phi}$ term using Euler's formula to be $\cos\phi + i \sin\phi$. If you multiply this against the $\cos\phi$ factor already in the expression, and integrate from 0 to $2\pi$, you see that the term $\sin\phi\cos\phi$ integrate to zero (why?) and that $\cos^2\phi$ integrate to some constant (why? And I'll leave it to you to compute that constant yourself).
Now the innermost integral is taken care of, you can integrate the term $\sin^3\theta$ from $0$ to $\pi$. This gives you another constant (what is it?).
Lastly, you need to evaluate the $r$ integral. Note that after the previous two steps you are left with something that looks like $$ \text{Constant}\cdot \frac{1}{a^5} \int_0^\infty r^4 e^{-\frac{r}{a}} \mathrm{d}r $$ Now, you can rewrite it as $$ \text{Constant}\cdot \int_0^\infty \left(\frac{r}{a}\right)^4 e^{-\frac{r}{a}} \mathrm d \left(\frac{r}{a}\right) $$ so doing the change of variables $\rho = r/a$, your integral becomes $$ \text{Constant}\cdot \int_0^\infty \rho^4 e^{-\rho} \mathrm{d}\rho $$ This you can solve simply by repeated integration by parts, or by appealing to the Gamma function (whose values at positive integers are explicitly known).