If I have 2 random variables $\xi, \eta$ and $\forall n,m \ \mathbb E\xi^n\eta^m=\mathbb E\xi^n \mathbb E\eta^m$, does this imply that $\xi,\eta$ are independent? How to show it?



Statement about independence of random variables
The answer is no. Indeed, let $X:=\xi$ and $Y:=\eta$. Let $U$ and $V$ be any independent random variables (r.v.'s) with different distributions but with the same finite moments of all orders: 
$$EU^m=EV^m=:\mu\_m$$
for all natural $m$. A standard example of the distributions of such r.v.'s $U$ and $V$ is given in [the answer by saz](https://math.stackexchange.com/questions/1166637/do-moments-define-distributions). 


Let the cumulative distribution function (cdf) $F\_{X,Y}$ of the random pair $(X,Y)$ be the half-and-half mixture of the cdf's $F\_{U,V}$ and $F\_{V,U}$, so that 
\begin{equation\*}
 F\_{X,Y}(x,y)=\frac{F(x)G(y)+G(x)F(y)}2 \tag{1}
\end{equation\*}
for all real $x,y$, where $F$ and $G$ are the cdf's of $U$ and $V$, respectively. Then for the cdf's $F\_X$ and $F\_Y$ one has $F\_X=F\_Y=\frac{F+G}2$ and hence 
\begin{equation\*}
4[F\_X(x)F\_Y(y)-F\_{X,Y}(x,y)]=[F(x)-G(x)][F(y)-G(y)]\ne0
\end{equation\*}
for some real $x,y$, so that $X$ and $Y$ are not independent. 


However, 
\begin{equation\*}
 EX^mY^n=\tfrac12\,EU^m\, EV^n+\tfrac12\,EV^m\,EU^n=\mu\_m\mu\_n=EX^m\,EY^n \tag{2}
\end{equation\*}
for all natural $m,n$. 




---


**Details on (2), in response to comments by the OP:** To simplify the matter, here we may use the fact that in the mentioned standard example of the distributions of r.v.'s $U$ and $V$, these r.v.'s actually have probability density functions (pdf's), say $f$ and $g$. Then one may rewrite (1) as
\begin{equation\*}
 f\_{X,Y}(x,y)=\frac{f(x)g(y)+g(x)f(y)}2 \tag{1a}
\end{equation\*}
for all real $x,y$, where $f\_{X,Y}$ is the joint pdf of $(X,Y)$. So, 
\begin{multline\*}
 EX^mY^m=\int\_{-\infty}^\infty\int\_{-\infty}^\infty x^my^n f\_{X,Y}(x,y)\,dx\,dy \\ 
=\frac12\,\int\_{-\infty}^\infty x^m f(x)\,dx\ \int\_{-\infty}^\infty y^n g(y)\,dy+\frac12\,\int\_{-\infty}^\infty x^m g(x)\,dx\ \int\_{-\infty}^\infty y^n f(y)\,dy \\ 
=\tfrac12\,EU^m\, EV^n+\tfrac12\,EV^m\,EU^n
=\mu\_m\mu\_n=EX^m\,EY^n. \tag{2a}
\end{multline\*}


