Random reflections unexpectedly produce banded distributions Start with $p_1$ a random point on the origin-centered unit circle $C$.
At step $i$, select a random point $q_i$ on $C$, and a random mirror line
$M_i$ through $q_i$, and reflect $p_i$ in $M_i$ to point $p_{i+1}$.
(Note that if every $M_i$ passed through the origin rather than $q_i \in C$,
all $p_i$ would lie on $C$.)
In the example below, $p_1$ on $C$ reflects to $p_2$ near the origin,
$p_2$ reflects out to $p_3$ near $(1,1.5)$, etc.

          


          

Green connects $p_1,\ldots,p_5$. Dark lines are mirrors $M_i$.


I expected that repeating this process would produce some smooth
distribution around the origin, as the points reflect further and further out,
when the mirrors effectively pass through the origin at that expanded scale.
But I find not infrequently bands of density, as depicted below.

   


   

Each distribution includes $n=5000$ reflected points. The unit circle is red.
Approx. radii: $107,67,47,148$.


Extending $n$ to larger values seems to continue banding effects.
Here is the 3rd example above, the most uniform of those four, extended.
Its radius increases from $47$ to $225$:

          


          

Continuation of 3rd example above, to $n=25000$. Approx. radius $225$.



Q. What explains this radial clustering/banding behavior?

 A: Let $p_i$ be at distance $r_i > 2$ from the origin, WLOG assume $p_i = (r_i, 0)$. Let $q_i = (\cos \alpha, \sin \alpha)$, and $v_i = (\cos \beta, \sin \beta)$ be a direction vector of $M_i$. By choice $\alpha$ and $\beta$ are independent and equidistributed on $[0, 2 \pi)$. Simple calculations show that $||p_{i + 1}||^2 =: r_{i + 1}^2 = r_i^2 + 4r_i \sin \beta \sin \gamma + 4 \sin^2 \gamma$, where $\gamma = \alpha - \beta$. It follows straightforwardly that $|r_{i + 1} - r_i| \leq 2$, further ${\mathrm E} r_{i + 1} = r_i$. Indeed, $\beta$ is independent of $\gamma$, and, conditioned on $\gamma$, $r_{i + 1}$ is distributed symmetrically on $[r_i - 2 \sin \gamma, r_i + 2 \sin \gamma]$. (In fact, $r_{i + 1} - r_i$ has density $p(x) = \frac{\pi - 2 \arcsin(x/2)}{8}$, but we won't need that) Note that when $r_i < 2$, $r_{i + 1}$ is larger on average since it can "bounce" off zero.  This turns out to be wrong! Empirically $\mathrm{E}(r_{i + 1} - r_i) = \frac{1}{2r_i} + O(r_i^{-2})$. What follows is not technically correct, but I'll leave it as is until I figure out a proper way to fix it.
Now, assuming $r_i$ is large enough, the sequence $r_i, r_{i + 1}, \ldots, r_{i + k}$ is practically a martingale as long as no $r_j < 2$. Since $|r_{i + 1} - r_i| \leq 2$, Azuma's inequality would yield $\mathrm{Prob}(|r_{i + k} - r_i|) > t) \leq 2\exp(\frac{-t^2}{8k})$, hence a lot of subsequent $r_{i + k}$ are likely to be concentrated around $r_i$ (for instance, a very stupid bound shows that out of subsequent $20000$ points at least $99.95\%$ will stay within $|r_{i + k} - r_i| \leq 1000$). This could be made into a more quantitative statement, but I feel it still explains the phenomenon quite well.
A: (Too long for a comment)
Let $n_i$ denote the unit normal vector to the mirror line $M_i$. Then we can provide a simple recursive formula $p_{i+1} = p_i - 2\langle p_i-q_i, n_i \rangle n_i$. Since $p_i - \langle p_i, n_i \rangle n_i$ is orthogonal to $n_i$, it follows that
\begin{align*}
\| p_{i+1} \|
&= \| p_i - 2\langle q_i, n_i \rangle n_i \|.
\end{align*}
From this, we can realize the law of $\|p_i\|$'s by the following construction:


*

*Let $(\theta_i)_{i=1}^{\infty}$ and $(\varphi_i)_{i=1}^{\infty}$ i.i.d. having the uniform distribution over $[0, 2\pi)$.

*Set $R_1 = 1$ and $R_{i+1} = \sqrt{R_i^2 - 4R_i \cos\varphi_i \cos\theta_i + 4\cos^2\theta_i}$.
Then $(R_i)_{i=1}^{\infty} \stackrel{d}= (\|p_i\|)_{i=1}^{\infty}$. Here are some easy observations:


*

*If $R_i \geq 2$, then $R_{i+1} \geq R_i - 2\cos\varphi_i\cos\theta_i $

*$R_i \gg 1$, then $R_{i+1} = R_i - 2\cos\varphi_i\cos\theta_i + \mathcal{O}(R_i^{-1})$
This means that $(R_i)$ almost behave like a random walk, thus the 'banding' may be explained by the language of random walk, such as local times, although I am no expert on this topic.

Figure. A sample path of $(\|p_i\|)$ for $i \in [1, 10^6]$.
A: As said by Sangchui Lee, we are interested in $r_i = |p_i|$ and we assume that the reflections induce enougth random rotation such that the "bands" of radius $r$ just reveal the number of time $r_i$ was closed to $r$. And we have
$$| p_{i+1}| =|p_i -2\langle q_i, n_i \rangle n_i|$$ 
We define then a random march $X_i$ defined by $X_0=0$ and $$X_{i+1}=X_i-2\langle q_i, n_i \rangle n_i$$
We immediately have that $|p_i|$ has the same law as $|X_i|$ and that for large $i$, $X_i$ behave like the brownian motion on $\mathbb{R}^2$.
One can calculate that the variance of $2\langle q_i, n_i \rangle n_i$ is $\sigma^2 = 2$, the points should be then in a disk of randius $\sim \sqrt{n \sigma^2 }=100$ for $n=5000$ and $224$ for $n=25000$ . 
Conclusion for large $i$, $|p_i|$ is given by the Bessel process. And as Sangchui Lee said we are interested in the local time. 
https://en.wikipedia.org/wiki/Local_time_(mathematics). 
The Ray-Knigh theorem should be usefull here even if it is state for the Brownian motion and not for the Bessel process.
(original paper of Knigh https://www.jstor.org/stable/pdf/1993647.pdf)
