Central limit theorem for 3d rotations Let $X_i$ (where $i=1, \dots , n$) be independent and identically distributed 3d rotations. What is the distribution of $X_1X_2\dotsb X_n$ in the limit of large $n$?
I'm especially interested in the special case where $X_i$ have mean 0 and are highly concentrated i.e. they are very small rotations. Just to make this clear, an example of this in the 2d case would be the rotations $(x,y)\rightarrow(x\cos \epsilon -y\sin \epsilon, x\sin \epsilon +y\cos \epsilon)$ where $\epsilon$ is a random variable with mean 0 and small variance.
 A: About your 2D example: for $\epsilon$ of fixed variance, as $n\to\infty$, the composition will converge in distribution to the uniform (= Haar) measure on the orthogonal group. In less fancy words, the angle of the rotation will converge to the uniform measure on the circle. Note that this will hold whether $\epsilon$ has mean $0$ or not, just some non-degeneracy (typically, the additive group generated by the support $\epsilon$ should not be a discrete sub-group of the circle ...) and positive variance are enough.
In the more general case, typically (meaning under reasonable assumptions about the distribution) if you compose rotations around the origin with fixed variance, you converge to a uniform distribution on the orthogonal group (the Haar measure).
There will be two more interesting setups: the affine case (i.e. the rotation center is not always the same) and the "CLT" case when you choose $\epsilon_n$ together with $n$, something like $1/\sqrt n$. The first one I am not sure about; for the second one, you will end up with Brownian motion on the orthogonal group and the limiting distribution should be the exponential of a random, Gaussian, anti-symmetric matrix (i.e. a Gaussian on the corresponding Lie algebra).
A: 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}
