$\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 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, \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*} \sqrt n\,\sup_{z\in\R}\Big|P\Big(\frac{XY-EXY}{\sqrt{Var\,XY}}\le z\Big)-P(Z\le z)\Big|\le C_p \end{equation*} for some real $C_p>0$ depending only on $p\in(0,1)$ and all natural $n$.

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

$\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 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, \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/\sqrt 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*} \sqrt n\,\sup_{z\in\R}\Big|P\Big(\frac{XY-EXY}{\sqrt{Var\,XY}}\le z\Big)-P(Z\le z)\Big|\le C_p \end{equation*} for some real $C_p>0$ depending only on $p\in(0,1)$ and all natural $n$.

$\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 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, \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*} \sqrt n\,\sup_{z\in\R}\Big|P\Big(\frac{XY-EXY}{\sqrt{Var\,XY}}\le z\Big)-P(Z\le z)\Big|\le C_p \end{equation*} for some real $C_p>0$ depending only on $p\in(0,1)$ and all natural $n$.