$\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](https://math.stackexchange.com/a/2599580/96609).


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



Normal distribution by successive approximation?
In the rewritten form it doesn't require any ingenuity at all. For brevity, introduce the notation
$$
(Th)(s)=\frac 1\pi\int\_0^s\frac{h(s-u)h(u)}{\sqrt{u(s-u)}}\,du=\frac 1\pi\int\_0^1\frac{h(s(1-v))h(sv)}{\sqrt{v(1-v)}}\,dv
$$
Suppose that $h\_0$ is any real valued Lipschitz function on $[0,+\infty)$ (the finiteness of the second moment guarantees Lipschitzness for your $h\_0$ and, if you really want to discuss it, I can show how this particular condition can be dropped) such that $\|h\_0\|\_\infty=h\_0(0)=1$ and there exists $h\_0'(0)=-c< 0$ (the latter two conditions are essential). Then the iterations $h\_{n+1}=Th\_n$ converge to $e^{-cx}$ uniformly on compact subsets of $[0,+\infty)$.


The proof consists of several easy observations:


1. $\|Th\|\_\infty\le\|h\|\_\infty^2$ and, thereby, $\|h\_n\|\_\infty\le 1$ for all $n$.
2. If $h$ is bounded and $L$-Lipschitz, then $Th$ is $L\|h\|\_\infty$-Lipschitz.


Indeed,
$$
Th(s)-Th(S)=\frac1\pi\int\_0^1\frac{dv}{\sqrt{v(1-v)}}[h(s(1-v))h(sv)-h(S(1-v))h(Sv)]
$$
and
$$
|h(s(1-v))h(sv)-h(S(1-v))h(Sv)|
\\
\le\|h\|\_\infty[|h(s(1-v))-h(S(1-v))|+|h(sv)-h(Sv)|]
\\
\le\|h\|\_\infty L[|s-S|(1-v+v)]=\|h\|\_\infty L|s-S|\,.
$$
Thus all $h\_n$ are Lipschitz with the same Lipschitz constant $L$ as $h\_0$.


3. If $A\in\mathbb R$ and $h(s)\ge e^{-As}$ on $[0,s\_0]$, then $(Th)(s)\ge e^{-As}$ on $[0,s\_0]$.
3') If $a\in\mathbb R$ and $0\le h(s)\le e^{-as}$ on $[0,s\_0]$, then $(Th)(s)\le e^{-as}$ on $[0,s\_0]$.


(all exponential functions are fixed points of $T$ and we have monotonicity as long as $h$ stays non-negative)


4. Now choose any $0<b<a<c<A<B$ and choose $s\_0>0$ so that $e^{-As}\le h\_0(s)\le e^{-as}$ on $[0,s\_0]$ (by the derivative at zero condition such $s\_0$ exists).


Consider the largest $S\_n$ such that $e^{-Bs}\le h\_n(s)\le e^{-bs}$ on $[0,S\_n]$. We have $S\_0\ge s\_0$ and $S\_{n+1}\ge S\_n$. We want to improve the latter trivial inequality to some quantitative advance $S\_{n+1}\ge S\_n+\delta(S\_n)$ where $\delta>0$ is separated from $0$ on any compact subinterval of $[s\_0,+\infty)$. That is quite easy:
$$
h\_{n+1}(S\_n)=(Th\_n)(S\_n)=\frac1\pi\int\_0^{S\_n}\frac{h\_n(S\_n-u)h\_n(u)}{\sqrt{u(S\_n-u)}}\,du\ge
\\
\frac1\pi\int\_0^{S\_n}\frac{e^{-B(S\_n-u)}e^{-Bu}}{\sqrt{u(S\_n-u)}}\,du
+\frac1\pi\int\_0^{s\_0}\frac{e^{-B(S\_n-u)}[e^{-Au}-e^{-Bu}]}{\sqrt{u(S\_n-u)}}\,du
\\
\ge e^{-BS\_n}+e^{-BS\_n}S\_n^{-1}\int\_0^{s\_0}[e^{-Au}-e^{-Bu}]\,du=e^{-BS\_n}+\Delta(S\_n)
$$
Then, by the Lipschitz property of both $h\_{n+1}$ and $e^{-Bs}$, the inequality $h\_{n+1}(s)\ge e^{-Bs}$ persists on $[S\_n,S\_n+\delta(S\_n)]$ with $\delta(S\_n)=\frac{\Delta(S\_n)}{B+L}$, say. The extension of the upper bound is similar.


The outcome is that the double inequality $e^{-Bs}\le h\_n(s)\le e^{-bs}$ propagates from $[0,s\_0]$ to the entire real line. Since $0<b<a<c<A<B$ were arbitrary, we conclude that $h\_n$ tend to $e^{-cs}$ uniformly on compact intervals, finishing the story.


