Asymptotic expansion of nonlinear Gaussian transformation in terms of covariance I'm reading this paper and on page 8 the authors state without proof an asymptotic expansion of a multivariate Gaussian integral in terms of the covariance obtained by applying what they call the "Lebesgue theorem" (my guess is  that they mean the Riemann-Lebesgue lemma).
Reformulated, the problem is as follows: Let $X_n, Y_n$ be standard Gaussians such that their covariance  $A(n)=\operatorname{E} (X_nY_n)$ tends to zero polynomially fast as $n \to \infty$. Furthermore, let $g$ be an even function with polynomial growth, such that $\operatorname{E}g(X_n)=0$. I would like to understand why 
$$\operatorname{E} \left( g(X_n) g(Y_n) \right) \sim
\operatorname{E} \left( g(X)g(Y) (X^2Y^2 - X^2 - Y^2 +1 ) \right) \frac{1}{2} A^2(n)$$
as $n \to \infty$, where $X, Y$ are two independent standard Gaussians.
Can somebody please hint me in the right direction?
Thanks a lot!
 A: Let $r:=r_n:=A(n)$, so that $r\to0$.  The right-hand side of the asymptotic relation in question is 
\begin{equation*}
 \text{rhs}=\frac{r^2}2\,Eg(X)(X^2-1)\;Eg(Y)(Y^2-1)\Big[=\frac{r^2}2\,(Eg(X)X^2)^2\Big], 
\end{equation*}
since $X,Y$ are iid [and $Eg(X)=0$]. 
The left-hand side of that asymptotic relation is 
\begin{equation*}
 \text{lhs}=\int_{-\infty}^\infty\int_{-\infty}^\infty g(x)g(y)h(r;x,y)\,dx\,dy, \tag{1}
\end{equation*}
where
\begin{equation}
 h(r;x,y):=\frac1{2 \pi  \sqrt{1-r^2}}\,\exp\Big\{-\frac{x^2-2 r x y+y^2}{2 \left(1-r^2\right)}\Big\}
\end{equation}
is the joint pdf of $(X_n,Y_n)$. Letting $f$ denote the standard normal pdf, we have 
\begin{equation*}
 h(0;x,y)=f(x)f(y),\quad h'_r(0;x,y)=xf(x)yf(y), 
\end{equation*}
\begin{equation*}
h''_{rr}(0;x,y)=(x^2-1)f(x)(y^2-1)f(y), 
\end{equation*}
\begin{equation*}
 h'''_{rrr}(r;x,y)=(1-r^2)^{-6}P(r,x,y)h(r;x,y), 
\end{equation*}
where $P(r,x,y)$ is a certain polynomial in $r,x,y$. 
Also, $Eg(X)=0$ and $Eg(X)X=0$ (since $g(X)X$ is odd in $X$, whereas $g$ grows slowly enough); here, as in the question, $X$ and $Y$ are iid standard normal random variables. 
So, using (1) with the Maclaurin expansion 
\begin{equation*}
 h(r;x,y)=h(0;x,y)+rh'_r(0;x,y)+\frac{r^2}2\,h''_{rr}(0;x,y)
 +\frac{r^3}2\,\int_0^1 h'''_{rrr}(\theta r;x,y)(1-\theta)^2\,d\theta, 
\end{equation*}
for $r\to0$ we have
\begin{align*}
 \text{lhs}  
 &=Eg(X)\,Eg(Y)+rEg(X)X\;Eg(Y)Y+\frac{r^2}2\,Eg(X)(X^2-1)\;Eg(Y)(Y^2-1)+O(r^3) \\ 
 &=\frac{r^2}2\,Eg(X)(X^2-1)\;Eg(Y)(Y^2-1)+O(r^3)\sim\text{rhs}
\end{align*}
-- assuming $Eg(X)X^2\ne0$. 
