The problem is related to the paper 1509.02093 by Oh and Thomann, where the authors considered the 2D Wick ordered NLS.
Let $g=a+ib$ be a complex number. Then it is claimed (see (2.7) in the paper and the text therein below) that the following identity (slightly different from the original one up to prefactors) holds for any $m\in\mathbb{N}$:
$$ \sum_{\ell=0}^m \begin{pmatrix}m\\ \ell\end{pmatrix}H_{2\ell}(a) H_{2(m-\ell)}(b)=P_m(a^2+b^2),$$
where $H_k$ is the $k$-th Hermite polynomial and $P_m$ is some polynomial of $m$-th order. This identity is simply stated without a detailed proof, but the authors claimed that it can be proved by using the ratational invariance of the complex Gaussian distributions. In fact, it was originally assumed that the numbers $a$ and $b$ were i.i.d. Gaussian distributions. But as one can observe from the examples computed on page 14, the identity clearly holds for all real numbers $a$ and $b$. As also pointed out by the authors, it is not quite obvious to derive the identity without using the probabilistic tools. However, I can neither find a probabilistic nor a deterministic proof. Some of my attempts would be to replace the numbers $a$ and $b$ by Gaussians with prefactors and then take expectation, but this was without success.
It would be nice that some hints with more details can be given.
Update 1: The original formula was wrong and a binomial factor was missing. Now corrected.
Update 2: I think that Lemma 2.1 in the paper has actually proved the claim. But the proof does not really make use of the rotational invariance of the complex Gaussian (in my opinion). It would still be nice that one can reveal the connection of the rotational invariance of the Gaussian and the claim.
Best Answer
$$\begin{eqnarray*} \textrm{LHS} &=& \exp\left(\frac{a^2+b^2}2\right) \left(\frac{\partial^2}{\partial a^2} + \frac{\partial^2}{\partial b^2}\right)^m \exp\left(-\frac{a^2+b^2}2\right) \\ &=& \exp\left(\frac{a^2+b^2}2\right) \Delta^m \exp\left(-\frac{a^2+b^2}2\right) \end{eqnarray*}$$ where $\Delta = \nabla^2$ is the Laplacian operator.
In polar coordinates $a = r \cos \theta$, $b = r \sin \theta$, $\Delta = \frac{\partial^2}{\partial r^2} + \frac 1r \frac{\partial}{\partial r} + \frac 1{r^2}\frac {\partial^2}{\partial \theta^2}$, but because the original sum $\exp(\tfrac{r^2}2) \Delta^m \exp(-\tfrac{r^2}2)$ is independent of $\theta$, the $\frac 1{r^2}\frac {\partial^2}{\partial \theta^2}$ component is zero. I think this explains the comment on rotation-invariance.
We get $$\exp\left(\frac{r^2}2\right) \left( \frac{\partial^2}{\partial r^2} + \frac 1r \frac{\partial}{\partial r} \right)^m \exp\left(-\frac{r^2}2\right)$$ and note that $$\begin{eqnarray*} && \left( \frac{\partial^2}{\partial r^2} + \frac 1r \frac{\partial}{\partial r} \right) Q(r^2) \exp\left(-\frac{r^2}2\right) \\ &=& \left( \frac{\partial}{\partial r} + \frac 1r \right) \frac{\partial}{\partial r} Q(r^2) \exp\left(-\frac{r^2}2\right) \\ &=& \left( \frac{\partial}{\partial r} + \frac 1r \right) \left( 2r Q'(r^2) - r Q(r^2) \right) \exp\left(-\frac{r^2}2\right) \\ % &=& \left((r^2-2) Q(r^2) + 4(1-r^2) Q'(r^2) + 4r^2 Q''(r^2) \right) \exp\left(-\frac{r^2}2\right) \\ % \end{eqnarray*}$$ so we derive $$P_{m+1}(x) = (x-2) P_m(x) + 4(1-x) P_m'(x) + 4x P_m''(x)$$ which might perhaps (I'm not sure, but I find it suggestive) be better regrouped as $$P_{m+1}(x) = (x-1) (P_m(x) - 4 P_m'(x)) - (P_m(x) - 4x P_m''(x))$$