I have a vector $v_{0} \in \mathbb{R}^{N}$ on the unit sphere. The vector is constructed by sampling each element from a standard normal distribution and normalizing the final result.
Now what I want to do is generate other vectors $v \in \mathbb{R}^{N}$ on the unit sphere such that the "angle" between the vectors is uniformly distributed, that is $cos^{-1}(v \cdot v_{0})$ is uniformly distributed. To do this, I've been generating random rotation matrices in $SO(N)$ using scipy.stats.special_ortho_group and rotating $v_{0}$ with these matrices to get $v$.
I'm observing a very curious phenomenon with this procedure, which is that the resulting distribution of $cos^{-1}(v \cdot v_{0})$ is normally distributed! And it appears to be centered at 90 degrees. I've attached a histogram of $cos^{-1}(v \cdot v_{0})$ for 3000 samples of this procedure with $N = 22$.
Could someone provide intuition why the result is normally distributed? How could I generate a result which is uniformly distributed from $0$ to $2\pi$?
What you get by applying to a vector $v$ a matrix from a uniform distribution of rotation matrices is a uniform distribution of vectors. Notice that the vast majority of vectors form wide-angles with respect to any given vector.
There's simply more space at the equator of the sphere than there is at the top or bottom. So more vectors tend to form wide-angles with the top vector (or any other vector).
You do however get a uniform distribution of angles for $N=2$. To get a uniform distribution of angles with respect to $v_{0}$ at higher dimensions, you need to form new vectors from your $v$ vectors: $$v_{1} = v- \langle v,v_{0} \rangle v_{0}\\v_{2} = \frac{v_{1}}{|v_{1}|}$$
You will then have $v_{0}\perp v_{2}$. Then you select a random number from $\theta \in [0,2\pi]$ and compute $$v_{3} = v_{0} \cos(\theta) + v_{2} \sin(\theta)$$ The $v_{3}$ should give you uniform angles with respect to the original vector $v_{0}$.