$\newcommand\R{\mathbb R}\newcommand\la\lambda$It is well known and easy to see that the rotationally invariant 
product of two probability measures on $\R$ has to be a Gaussian (or Dirac) measure; see e.g. [this answer][1]. 

This appears to make the conjecture below somewhat plausible. 

Let $\mu$ be any probability measure on $\R^2$ with a finite nonzero covariance matrix. Let $\mu_1:=\mu$. For each natural $n$, consider the following three-step procedure: 

1. let $$\la_n:=\mu_n^{(1)}\otimes\mu_n^{(2)},$$
where $\mu_n^{(1)}$ and $\mu_n^{(2)}$ are the marginals of $\mu_n$; 
2. let $$\nu_n:=\frac1{2\pi}\int_0^{2\pi}\la_n R_t\,dt,$$
where  $\la_n R_t$ is the pushforward measure obtained from $\la_n$ by the rotation about the origin through angle $-t$;
3. let $\mu_{n+1}$ be obtained by rescaling the probability measure $\nu_n$ so that the covariance matrix of $\mu_{n+1}$ be the unit matrix. 

**Conjecture:** $\mu_n$ converges weakly (as $n\to\infty$) to the standard Gaussian measure on $\R^2$. 

Is this conjecture true? 

*Comment:* Perhaps Step 3, the rescaling, is not essential. Of course, if we have the convergence to a (nondegenerate) Gaussian measure without Step 3, then we have such a convergence with Step 3 as well. 

--- 

We can restate the problem (without rescaling) analytically as follows. Let $f_n$ denote the characteristic function of $\nu_n$, so that 
$$f_n(u,v)=\int_{\R^2}\nu_n(dx\times dy)e^{i(ux+uy)} \\
=\int_{\R^2}\nu_n(dx\times dy)\cos(ux+uy)\quad \text{(by symmetry)}$$
for all real $u$ and $v$. Then 
$$f_n(u,v)=g_n\big(\sqrt{u^2+v^2}\big)$$
for some function $g_n\colon[0,\infty)\to\mathbb R$ and all real $u$ and $v$. Then for all natural $n$ and all real $r\ge0$ 
$$g_{n+1}(r)=\frac2\pi\int_0^{\pi/2}
dt\,g_n(r\cos t)g_n(r\sin t). \tag{1}$$
We want to show that 
$g_n(r)\to e^{-c^2 r^2/2}$ for some $c\in(0,\infty)$ and all real $r\ge0$. 

So, analytically, the problem may be viewed as one of stability of a (nonlinear) integral equation or as one of solving such an integral equation by iterations. 

---

Using substitutions $g_n(r)=h_n(r^2)$ and $r^2=s$, we can rewrite (1) as 
$$h_{n+1}(s)=\frac1\pi\int_0^s
du\,\frac{h_n(s-u)}{\sqrt{s-u}}\frac{h_n(u)}{\sqrt{u}}$$
and then as 
$$\pi H_{n+1}(s)\sqrt s=[(H_n*H_n)(s)=]\int_0^s
du\,H_n(s-u)H_n(u)$$
for all natural $n$ and all real $s\ge0$, where $H_n(u):=h_n(u)/\sqrt u$. 
We want to show that $H_n(u)\to e^{-c^2 u/2}/\sqrt u$ for some $c\in(0,\infty)$ and all real $u>0$. 



  [1]: https://math.stackexchange.com/a/2599580/96609