Let $X_1$ and $X_2$ be uniform on $n$-spheres. What is the distribution of $\| X_1+X_2\|$? Suppose we have two  independent random variables $X_1$ and $X_2$ distribution on $n-1$-sphere of radius $r_1 $ and radius $r_2$, respectivly. Assume $r_1>r_2$. 
Recall, that the  $n-1$-sphere of radius $r$ is defined as 
\begin{align}
S_{n-1}= \{ x \in \mathbb{R}^n : \|x\|=r \}.
\end{align} 
We have to find the distribution of 
\begin{align}
U=X_1+X_2
\end{align} 
We can see that $U$ will be distributed on an annulus 
\begin{align}
A=\{ x:  r_1-r_2 \le \| x\|\le r_1+r_2  \}
\end{align} 
It is not difficult to see that $U$ has a uniform spherical angle.
Therefore, the question is what is the distribution of the magnitude of $U$ that is $\| U\|$?
This question is an extension of the question previously asked here .
For the bounty: I would like to see the exact expression for the distribution of $U$. 
 A: It is known that the uniform distribution on the unit $(n-1)$-sphere can be represented as a standard multivariate Gaussian divided by its norm. Therefore,
\begin{align*}
X_1 \overset{\mathcal{D}}{=} r_1\frac{Z_1}{\|Z_1\|} \qquad \text{and} \qquad X_2 \overset{\mathcal{D}}{=} r_2\frac{Z_2}{\|Z_2\|}
\end{align*}
where $Z_1, Z_2 \overset{\text{iid}}{\sim} N(\mathbf{0}_n, I_{n\times n})$. Hence,
\begin{align*}
\|U\| \overset{\mathcal{D}}{=} \left\|r_1\frac{Z_1}{\|Z_1\|} + r_2\frac{Z_2}{\|Z_2\|}\right\| = \sqrt{r_1^2 + r_2^2 + 2r_1r_2 \frac{Z_1^\intercal Z_2}{\|Z_1\|\|Z_2\|}}
\end{align*}
The next step is to find the distribution of $P \overset{\text{def}}{=} \frac{Z_1^\intercal Z_2}{\|Z_1\|\|Z_2\|}$, which takes the form of Pearson's correlation coefficient, except we are not subtracting out the sample mean in the variance/covariance calculation. You can actually show that
\begin{align*}
T \overset{\text{def}}{=}  \frac{P}{\sqrt{1-P^2}} \sim t_{n-1}/\sqrt{n-1}
\end{align*}
where $t_n$ is the t-distribution with $n$ degrees of freedom. This follows from the proof in Hotelling's "New Light on the Correlation Coefficient and its Transforms" (1953), changed from $n-2$ to $n-1$ because of not needing to estimate the mean.
A: Continuing on from Tom Chen's answer, let $X\sim t_{n-1}$ and
$$
f(x) = \frac x{(1-x^2)^{1/2}}.
$$
Then $f(P)=_d (n-1)^{-1/2}X$, so that $\|U\|=_d g(X)$, where
$$
g(x) = \sqrt{r_1^2 + r_2^2 + 2r_1r_2 f^{-1}((n-1)^{-1/2}x)}\,.
$$
Note that $X$ has density
$$
\varphi(x) = C_n\left({1 + \frac{t^2}{n-1}}\right)^{-n/2},
$$
where
$$
C_n = \frac{\Gamma(n/2)}{\sqrt{(n-1)\pi}\,\Gamma((n-1)/2)}.
$$
Thus, the density of $\|U\|$ is
$$
h(y) = \left({\frac d{dy}(g^{-1}(y))}\right)\varphi(g^{-1}(y)),
$$
defined for $r_1-r_2\le y\le r_1+r_2$. If we let
$$
\gamma_y = \frac{y^2 - r_1^2 - r_2^2}{2r_1r_2},
$$
then $g^{-1}(y)=\sqrt{n-1}\,f(\gamma_y)$, so that
\begin{align}
\frac d{dy}(g^{-1}(y))
&= \sqrt {n-1}\,f'(\gamma_y)\left({\frac y{r_1r_2}}\right)\\
&= \frac{\sqrt {n-1}}{r_1r_2}\frac y{(1 - \gamma_y^2)^{3/2}}.
\end{align}
Also,
\begin{align}
\varphi(g^{-1}(y))
&= C_n(1 + f^2(\gamma_y))^{-n/2}\\
&= C_n\left({\frac1{1 - \gamma_y^2}}\right)^{-n/2}\\
&= C_n(1 - \gamma_y^2)^{n/2}.
\end{align}
Thus,
$$
h(y) = \frac{\sqrt {n-1}\,C_n}{r_1r_2}\,y\,(1 - \gamma_y^2)^{(n-3)/2}.
$$
Putting it all together, the density of $\|U\|$ is
$$
h(y) = \frac{\Gamma(n/2)}{r_1r_2\sqrt\pi\,\Gamma((n-1)/2)}
\,y\,\left({1 - \left({
  \frac{y^2 - r_1^2 - r_2^2}{2r_1r_2}
}\right)^2}\right)^{(n-3)/2}
$$
for $r_1-r_2\le y\le r_1+r_2$, and $0$ otherwise.
EDIT:
Regarding the comment, when $n=3$, since $\Gamma(3/2)=\sqrt\pi\,/2$ and $\Gamma(1)=1$, this reduces to
$$
h(y) = \frac1{2r_1r_2}y.
$$
We then have
\begin{align}
\int_{r_1-r_2}^{r_1+r_2} h(y)\,dy
&= \frac1{4r_1r_2}y^2 \bigg|_{r_1-r_2}^{r_1+r_2}\\
&= \frac1{4r_1r_2}((r_1+r_2)^2 - (r_1-r_2)^2)\\
&= \frac1{4r_1r_2}(4r_1r_2) = 1.
\end{align}
