$\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.

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:

- let $$\la_n:=\mu_n^{(1)}\otimes\mu_n^{(2)},$$ where $\mu_n^{(1)}$ and $\mu_n^{(2)}$ are the marginals of $\mu_n$;
- 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$;
- 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$.

2more comments