If $X,Y\sim N(0,1)$ with $E XY=\rho$ covariance, what is the expectation of their product after applying probabilist's Hermite polynomials $E_{X,Y} H_n(X) H_k(Y)$?
My initial guess is $\rho^{2k}\delta_{nk}$ where $\delta$ is dirac's delta. This would hold in the cases where $\rho=0,\rho=1$, ie, independent, or copy, as a direct result of orthogonality of Hermite's polynomials. However, I cannot prove or refute it in the general case.
The property can be deduced from Mehler's formula. The formula states that $$ E(x,y)=\frac{1}{\sqrt{1-\rho^2}}\exp\left(-\frac{\rho^2(x^2+y^2)-2xy\rho}{2(1-\rho^2)}\right) = \sum_{m=0}^\infty \frac{\rho^m}{m!}H_m(x)H_m(y) $$ Observe that $E(x,y)$ is equal to $p(x,y)/p(x)p(y)$, where $p(x,y)$ is the joint PDF of $(X,Y)$, and $p(x),p(y)$ are PDF of $X$ and $Y$ respectively. Therefore, we can take the expectation using the expansion $$ \mathbf{E}_{X,Y} [H_n(X) H_k(Y)]= \int H_n(x)H_k(y)p(x,y)dxdy \\ = \sum_{m=0}^\infty \frac{\rho^m}{m!}\int H_n(x)H_k(y) H_m(x)H_m(y) p(x)p(y) dx dy \\ = \sum_{m=0}^\infty \frac{\rho^m}{m!}\int H_n(x)H_m(x)dx\int H_k(y)H_m(y)dy \\ = (2\pi) n! \rho^n \delta_{nk} $$ where in the last line I used the orthogonality property $E_x H_k(x) H_n(x)=\sqrt{2\pi} n!\delta_{nk}$
This derivation has some constant factor $(2\pi)n!$ plus a different exponent $\rho^k$ than the one guessed in the question.