This is the random walk on a sphere Roberts & Ursell (1960).
Suppose that $X_i$ are rotations by an angle $\alpha$ on the sphere in a random direction
with probability $p(\mu)d\mu$ where $\mu=\cos\alpha$. The
probability density on the sphere after the action of $X_n$ can be expressed with the probability density of the preceding step as
\begin{align}
\rho_n(\vec{r}) =
\frac{1}{2\pi}\int_{-1}^1 d \mu\int_{S_2} d \vec{r}' p(\mu)
\delta(\vec{r} \cdot \vec{r}'- \mu)\,\rho_{n-1}(\vec{r}')\,.
\end{align}
This equation is linear in $\rho$ and can be solved using the
eigenbasis of the corresponding linear operator.
Next we show that if all steps have the same angular value $\alpha$, the eigenfunctions are
the spherical harmonics $Y_{\ell m}$ with eigenvalue $P_{\ell}(\mu)$.
We will use the following identities of the Legendre polynomials,
\begin{align}
&\delta(\cos \gamma - \mu)= \sum_{\ell=0}^{\infty}
\frac{2\ell+1}{2}P_{\ell}(\cos\gamma)P_{\ell}(\mu)\,,\\
&P_{\ell}(\cos\gamma)
=
\frac{4\pi}{2\ell+1}\sum_{m=-\ell}^{\ell}
Y_{\ell m}(\vec{r})Y_{\ell m}^*(\vec{r}')\,.
\end{align}
where $Y_{\ell m}(\vec{r})$ are orthonormal spherical harmonics. Substituting
\begin{align}
\rho_n(\vec{r}) = \sum_{\ell=0}^{\infty}\sum_{m=-\ell}^{\ell}
P_{\ell}(\mu)\, Y_{\ell m}(\vec{r}) \int_{S_2} d\vec{r}'\,Y_{\ell
m}^*(\vec{r}')\rho_{n-1}(\vec{r}').
\end{align}
Expanding the initial distribution in this basis as
\begin{align}
\rho_0(\vec{r}) = \sum_{\ell, m}a_{\ell m,0}Y_{\ell m}(\vec{r})
\end{align}
the distribution after a single step is
\begin{align}
\rho_1(\vec{r}) = \sum_{\ell,m}
\langle P_{\ell}(\mu) \rangle\, a_{\ell m,0} Y_{\ell m}(\vec{r}),
\end{align}
where $\langle P_{\ell}(\cos \alpha) \rangle = \int_{-1}^1 P_{\ell}(\mu) p(\mu) d \mu$.
After the $n^{\rm th}$ step we have
\begin{align}
\rho_{n}(\vec{r}) = \sum_{\ell, m} \langle P_{\ell}(\mu) \rangle^{n} a_{\ell m,0}Y_{\ell m}(\vec{r}).
\end{align}
Since $|\langle P_\ell(\mu)\rangle| \leq 1$ for $\ell>0$,
each $a_{\ell m}$ multipole moment with $\ell>0$ decays exponentially in the
number of steps from their initial value $a_{\ell m,0}$
as $\exp[n\ln|\langle P_\ell(\mu)\rangle|]$; the
system ``isotropizes'' with a decay time of $\Delta t /\ln |\langle
P_{\ell}(\mu) \rangle|$ where $\Delta t$ is the timestep.
The Green's function corresponding to an initial density
$\rho_0$ that is concentrated at the $\theta=0$ pole corresponds to
$a_{\ell m} =\sqrt{(2\ell + 1) /(4\pi)}\delta_{m,0}$. Thus the
probability for the angle $\theta$
between the initial and final position after $n$ steps is given by
\begin{align}
\rho_n(\theta_n) d\cos\theta_n = \sum_{\ell=0}^\infty
\frac{2\ell +1}{2}\, \langle P_{\ell}(\mu)\rangle^n\,
P_{\ell}(\cos\theta_n)d\cos\theta_n\,,
\end{align}
which implies that
\begin{align}
\langle P_{\ell}(\cos\theta_n)\rangle \equiv \int_{-1}^{1} P_{\ell}(\cos\theta_n)\,\rho_n(\theta_n)d \cos\theta_n =
\langle P_{\ell}(\mu)\rangle^n\,.
\end{align}
In the limiting case of Brownian motion the angular step $\alpha$ and
the timestep $\Delta t$ both approach zero with $\alpha^2\sim \Delta
t$. In this limit
$P_{\ell}(\mu) \approx 1-\frac14 \ell (\ell+1) \alpha^2$, and so
\begin{align}
\rho_{n}(\vec{r}) = \sum_{\ell,m} a_{\ell m} Y_{\ell m}(\vec{r}) e^{-\frac14 \ell (\ell+1) v}
\end{align}
where $v = n \langle \alpha^2\rangle = \langle\alpha^2\rangle t
/\Delta t $ is the variance of the corresponding planar motion.
Thus
\begin{align}
\langle P_{\ell}(\cos\theta_n) \rangle = e^{-\frac14 \ell(\ell+1) n \langle\alpha^2\rangle}=e^{-\frac14 \ell(\ell+1) v}.
\end{align}
Up to this point we have only calculated the expectation value of the action of $X_n\dots X_2 X_1$. The distribution function may be similarly expressed with the distribution of multipole moments $x_{\ell} = P_{\ell}(\cos \theta_n)$ as
\begin{align}
p(x_{\ell}) = \prod_k \int_{-1}^{1} d \mu_k P_{\ell}(\mu_k)\delta\left( x_{\ell} - \prod_i P_{\ell}(\mu_i) \right)\,.
\end{align}
