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?

1$\begingroup$ isn't this what is called subindependence ? $\endgroup$ – Carlo Beenakker Mar 13 at 16:11

1$\begingroup$ See this question on MSE. $\endgroup$ – Timothy Budd Mar 13 at 16:40
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.
Let the cumulative distribution function (cdf) $F_{X,Y}$ of the random pair $(X,Y)$ be the halfandhalf 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*}

$\begingroup$ Can you explain me please why $E X^mY^n=1/2EU^mEV^n+1/2EV^mEU^n$? @Iosif Pinelis $\endgroup$ – user157895564 Mar 13 at 18:10

$\begingroup$ @user157895564 : One way to do this is to use (1) to write the following: $EX^mY^m=\int_{\infty}^\infty\int_{\infty}^\infty x^my^n dF_{X,Y}(x,y) =\tfrac12\,\int_{\infty}^\infty x^m dF(x)\int_{\infty}^\infty y^n dG(y)+\tfrac12\,\int_{\infty}^\infty x^m dG(x)\int_{\infty}^\infty y^n dF(y)=\tfrac12\,EU^m\, EV^n+\tfrac12\,EV^m\,EU^n$. $\endgroup$ – Iosif Pinelis Mar 13 at 19:03

$\begingroup$ Why second equality is correct? @Iosif Pinelis $\endgroup$ – user157895564 Mar 13 at 19:46

$\begingroup$ @user157895564 : This equality immediately follows by (1) and the definitions of the LebesgueStieltjes integrals with respect to $dF_{X,Y}$, $dF$, and $dG$. To simplify this, we may use the fact that in the cited example of the distributions of r.v.'s $U$ and $V$, these r.v.'s actually have pdf's. I have added this detail to the answer. $\endgroup$ – Iosif Pinelis Mar 13 at 21:08

$\begingroup$ But it doesn't work, isn't it? $(F_{X,Y})'_x=f_{X,Y}$ only if $G=g,F=f$, but $f$ is always summable, but $F$ is not. Or I'm missed something? $\endgroup$ – user157895564 Mar 14 at 4:45