Number of random unit vectors which are less than theta apart Given $n$ unit vectors which are uniformly distributed on a unit sphere, what the expected number of groups of $k$ vectors which are within an angle $\theta$ of one another
For example, if I have $n=50$ unit vectors, what is the expected number of triplets of vectors ($k=3$) that are within $\theta=20^\circ$ of one another? From simulation it seems like ~13.5

As Robert Israel says, this is probably difficult to compute exactly for $k\gt3$, but can be done for $k$ up to $3$.
For $k=1$, the answer is trivially $n$.
For $k=2$, the probability for two vectors to be within $\theta$ of each other is the surface area of a spherical cap with angle $\theta$, which is $2\pi(1-\cos\theta)$. There are $\binom n2$ pairs of vectors, so the expected number of pairs within $\theta$ is $n(n-1)\pi(1-\cos\theta)$.
For $k=3$, we can consider two points $\alpha$ apart and determine the probability that the third point is at distance at most $\alpha$ from both of them. The area in which the third point must lie is a sort of lune, but it’s not a spherical lune, since it’s formed by small circles, whereas a spherical lune is formed by great circles.
To compute the area of this lune, imagine one of the two points at the north pole $(0,0,1)$, with a spherical cap at angle $\alpha$ around it, and the second point on the boundary of the cap at $(\sin\alpha,0,\cos\alpha)$. Now rotate this arrangement by $\frac\alpha2$ about the positive $y$ axis, so that the first and second point are now at $(-\sin\frac\alpha2,0,\cos\frac\alpha2)$ and $(\sin\frac\alpha2,0,\cos\frac\alpha2)$, respectively. The points, and thus the lune, are now symmetric about the $x$-$z$ plane. One half of the lune is the part of the cap that has retained a positive $x$ coordinate after the rotation. The $x$ coordinate of the point $(\sin\beta\cos\phi,\sin\beta\sin\phi,\cos\beta)$ after the rotation is $\cos\frac\alpha2\sin\beta\cos\phi-\sin\frac\alpha2\cos\beta$. This is positive for
$$
\cos\phi\gt\frac{\tan\frac\alpha2}{\tan\beta}\;.
$$
Thus the area of the lune is
$$
4\int_\frac\alpha2^\alpha\sin\beta \mathrm d\beta\arccos\frac{\tan\frac\alpha2}{\tan\beta}\;.
$$
This integral can be computed using Sage:
alpha,beta = var('alpha,beta')
assume(alpha > 0)
integral (arccos(tan(alpha/2) / tan (beta)) * sin (beta),beta,alpha/2,alpha)

Simplifying the output yields
$$
\int_\frac\alpha2^\alpha\sin\beta\,\mathrm d\beta\arccos\frac{\tan\frac\alpha2}{\tan\beta}
=\arctan\frac{\sqrt{\cos^2\frac\alpha2 + \sin^2\alpha - 1}}{\sin\frac\alpha2}-\cos\alpha\arccos\frac{\tan\frac\alpha2}{\tan\alpha}
\;.
$$
Any two of the three points are equally likely to be the pair at maximal distance, and these three events are mutually exclusive. Thus the probability that the maximal distance among the three points is $\alpha$ is $3$ times the probability that two particular points are $\alpha$ apart and the third point is closer to both of them, and that probability is the area just calculated over the total area $4\pi$. The density for two points to be $\alpha$ apart is $\frac12\sin\alpha$. Thus the desired probability that the three points are not more than $\theta$ apart is
$$
\frac3{2\pi}\int_0^\theta\sin\alpha\,\mathrm d\alpha\left(\arctan\frac{\sqrt{\cos^2\frac\alpha2 + \sin^2\alpha - 1}}{\sin\frac\alpha2}-\cos\alpha\arccos\frac{\tan\frac\alpha2}{\tan\alpha}
\right)\;.
$$
Somewhat miraculously, Wolfram|Alpha can solve this monster of an integral. After a lot of simplification, the result is
$$
\frac1{8\pi}\Bigg(
2(\cos\theta-4)\sqrt{2\cos\theta+1} + 3\cos2\theta\arccos\left(\frac{\tan\frac\theta2}{\tan\theta}\right)
-12\cos\theta\arctan\left(\frac{\sqrt{\cos\theta-\cos2\theta}}{\sqrt2\sin\frac\theta2}\right)-16\arcsin\left(\frac1{2\cos\frac\theta2}\right)-10\arctan\sqrt{2\cos\theta+1}-\frac{3\sqrt2\sqrt{\cos\theta-\cos2\theta}}{2\cos\theta+1}\left(\sin\frac\theta2+\cot\frac\theta2\left(2\sqrt{3-\tan^2\frac\theta2}\arcsin\left(\frac12\sqrt{3-\tan^2\frac\theta2}\right)-3\cos\frac\theta2\right)\right)+13\pi\Bigg)\;.
$$
This is valid for $\theta\le\frac{2\pi}3$. For $\alpha\gt\frac{2\pi}3$, the circles  around the two points at angle $\alpha$ no longer intersect (as indicated by the fact that the argument of the arccosine becomes $-1$), so there’s no longer a lune and it’s not possible for a point to be at distance more than $\alpha$ from both points. Thus we just have to subtract the area of the complementary spherical caps with angle $2\pi-\alpha$ around the antipodes of the two points, where a point would be at distance more than $\alpha$ from one (but not the other) of the two points:
$$
4\pi-2\cdot2\pi(1-\cos(\pi-\alpha))=4\pi\cos(\pi-\alpha)\;,
$$
so the probability for the third point to be at most $\alpha$ away from both of the two points is simply $\cos(\pi-\alpha)$. Multiplying by $3$ as above and integrating with the density $\frac12\sin\theta$ from $\theta$ to $\pi$ yields the desired probability for $\frac{2\pi}3\le\theta\le\pi$:
$$
1-3\cdot\frac12\int_{\theta}^\pi\sin\alpha\,\mathrm d\alpha\cos(\pi-\alpha)=1-\frac34\sin^2\theta\;.
$$
Here’s a plot of the probability for the entire range $0\le\theta\le\pi$ (generated with Sage):

Here are some specific values of the probability:
\begin{array}{c|c}
\theta&\text{probability}\\\hline
\frac\pi3&\frac1{8\pi}\left(5\pi-\sqrt2-12\arctan\sqrt2-\frac32\arccos\frac13\right)\\
\frac\pi2&\frac14-\frac1{4\pi}\\
\frac{2\pi}3&\frac7{16}\\
\end{array}
For your example $\theta=20^\circ=\frac\pi9$, the probability is about $5.37\cdot10^{-4}$. The expected number of triples in which the vectors are within $\alpha$ of each other is the number of triples, $\binom n3$, times this probability. For $n=50$, this is about $10.5$, which is not in agreement with your result of $13.5$. I checked the result with simulations using this Java code. So it seems that either I misunderstood the question or there’s a bug in your simulation.
