(I asked this question [on MSE](https://math.stackexchange.com/questions/4521591) 10 days ago, but I got no answer.)


Let $X$ and $Y$ be two independent identically distributed binomial random variables with parameters $n \in \mathbb{N}$ and $p \in (0,1)$. Let $Z := XY$ be their product.


Is it true or false that $\tilde{Z} := (Z - \mathbf{E}[Z]) / \sqrt{\mathbf{Var}[Z]}$ converges in distribution to a standard normal random variable (as $n \to \infty$) ?


At a first glance, I would be tempted to write $X = \sum\_{i=1}^n A\_i$ and $Y = \sum\_{i=1}^n B\_i$, where $A\_i$ and $B\_i$ are Bernoulli random variables, and then to apply the central limit theorem to $Z = \sum\_{i=1}^n \sum\_{j = 1} A\_i B\_j$... but $A\_i B\_j$ are *not* independent...


Thanks for any help


**P.S.1** For $p=1/2$, it is easy to check that $\mathbf{E}[Z] = n^2 / 4$ and $\mathbf{Var}[Z] = n^3 / 8 + n^2 / 16$. Moreover, expanding $(XY - n^2/4)^k$ with the binomial theorem and using the formula for the moments of the binomial distribution, I got that


$$\mathbf{E}\left[\tilde{Z}^k\right] = \frac1{(n^3 / 8 + n^2 / 16)^{k/2}} \sum\_{j=0}^k \binom{k}{j} \left(\sum\_{i=0}^j {j \brace i} (n)\_{i} (1/2)^i\right)^2 (-n^2/4)^{k-j}$$


where ${j \brace i}$ are Stirling numbers of second kind and $(n)\_{i}$ is a falling factorial. With this formula, I verified that the first 20 moments of $\tilde{Z}$ tend to the moments of a standard normal variable. However, I still do not know how to prove this for all moments.


**P.S.2** I think that one cannot prove the claim only using the fact that the (normalized) $X$ and $Y$ converge in distribution to normal variables. In fact, [it is known](https://mathworld.wolfram.com/NormalProductDistribution.html) that the product of two independent normal variables is not a normal variable.



Does the (normalized) product of two independent binomial variables converge in distribution to a normal variable?
$\newcommand{\R}{\mathbb R}\newcommand\ep\epsilon\newcommand\tsi{\tilde\sigma}$Yes, of course. This follows by the multivariate (here, bivariate) so-called delta method.


Indeed, we may assume that
\begin{equation\*}
 X=\sum\_{i=1}^n X\_i,\quad Y=\sum\_{i=1}^n Y\_i,
\end{equation\*}
where $X\_1,Y\_1,\dots,X\_n,Y\_n$ are independent identically distributed (iid) Bernoulli random variables (r.v.'s) with parameter $p\in(0,1)$.


For each $i\in[n]:=\{1,\dots,n\}$, let
\begin{equation\*}
 V\_i:=(X\_i-p,Y\_i-p),
\end{equation\*}
so that $V\_1,\dots,V\_n$ are iid zero-mean random vectors in $\R^2$. Then
\begin{equation\*}
 XY=n^2 f(\bar V)+n^2p^2, \tag{1}\label{1}
\end{equation\*}
where $\bar V:=\frac1n\,\sum\_{i=1}^n V\_i$ and for $(x,y)\in\R^2$
\begin{equation\*}
 f((x,y)):=f(x,y):=(x+p)(y+p)-p^2,
\end{equation\*}
so that $f(0,0)=0$ and for $L:=f'(0,0)$ we have
\begin{equation\*}
 L(x,y)=px+py
\end{equation\*}
and
\begin{equation\*}
 |f(x,y)-L(x,y)|=|xy|\le\tfrac12\,\|(x,y)\|^2,
\end{equation\*}
where $\|(x,y)\|:=\sqrt{x^2+y^2}$.
So, condition (2.1) of [this paper](https://projecteuclid.org/journals/electronic-journal-of-statistics/volume-10/issue-1/Optimal-order-bounds-on-the-rate-of-convergence-to-normality/10.1214/16-EJS1133.full) holds (for any real $\ep>0$ and $M\_\ep=1$).


Moreover, $v\_3:=E\|V\_1\|^3<\infty$ and
\begin{equation\*}
 \tsi:=\sqrt{EL(V\_1)^2}=\sqrt{p^2 E(X\_1-p+Y\_1-p)^2}=\sqrt{2p^3q}>0,
\end{equation\*}
where $q:=1-p$.
So, by \eqref{1} and Theorem 2.9 of the [same paper](https://projecteuclid.org/journals/electronic-journal-of-statistics/volume-10/issue-1/Optimal-order-bounds-on-the-rate-of-convergence-to-normality/10.1214/16-EJS1133.full),
\begin{equation\*}
 \frac{XY-n^2p^2}{n^2\sqrt{2p^3q/n}}
 =\frac{f(\bar V)}{\tsi/\sqrt n}\to Z\sim N(0,1) \tag{2}\label{2}
\end{equation\*}
(in distribution as $n\to\infty$).


Note also that $EXY=n^2p^2$ and
$$Var\,XY=2n^3p^3q+n^2p^2q^2=2n^3p^3q(1+O(1/n))\\ 
\sim2n^3p^3q=(n^2\sqrt{2p^3q/n})^2.$$
Thus, by \eqref{2},
\begin{equation\*}
 \frac{XY-EXY}{\sqrt{Var\,XY}}\to Z, 
\end{equation\*}
as desired.


Moreover, it follows from cited Theorem 2.9 that
\begin{equation\*}
 \sup\_{z\in\R}\Big|P\Big(\frac{XY-EXY}{\sqrt{Var\,XY}}\le z\Big)-P(Z\le z)\Big|\le 
 \frac{C\_p}{\sqrt n}
\end{equation\*}
for some real $C\_p>0$ depending only on $p\in(0,1)$ and all natural $n$.


