This $y = \frac{<x_1, x_2>}{\|x_1\|\|x_2\|}$ is the distribution of the $cos \theta$ where $\theta$ is known as the canonical angle/principal angle of two random vectors $x_1,x_2$. $cos\theta$ is known as the canonical correlation between $X_1,X_2$ since we know their joint distribution $(X_1,X_2)$ from independence.

Usually we discussed when $(X_1,X_2)$ are assumed to be eigenvectors corresponding to a random matrix. And the derivation actually does not require normality assumption and can be extended to higher $r$, see [Anderson] Chap.10 for example.

The derivation can be done as shown in [Anderson] using iterative optimization, however, there is a general treatment [Anderson2][Hsu] which derives the exact distribution of eigenvectors and therefore the moments associated under normality assumption of the sample.

Update: I have answer a updated version of this question using directional statistics here, now it can be addressed by using projected normal distribution: https://stats.stackexchange.com/questions/263896/moment-mgf-of-cosine-of-two-random-vectors

[Anderson]An Introduction to Multivariate Statistical Analysis,1958,Wiley

[Anderson2]Multiple discoveries: Distribution of roots of determinantal
equations
T.W. Anderson

[Bartlett]Bartlett, M. S. "The general canonical correlation distribution." The Annals of Mathematical Statistics (1947): 1-17.http://projecteuclid.org/download/pdf_1/euclid.aoms/1177730488

[Hsu]Hsu, P. L. "On the distribution of roots of certain determinantal equations." Annals of Human Genetics 9.3 (1939): 250-258.