Problem proving $\int_{-\infty}^{\infty}x^2e^{-x^2}H_n(x)H_m(x)dx$ I am trying to show that the value of the following integral is:
$$\int_{-\infty}^{\infty}x^2e^{-x^2}H_n(x)H_m(x)dx=2^{n-1}\pi^{\frac{1}{2}}(2n+1)n!\delta_{mn}+2^n\pi^{\frac{1}{2}}(n+2)!\delta_{(n+2)m}+2^{n-2}\pi^{\frac{1}{2}}n!\delta_{(n-2)m},$$
where $H_n$ is the $n$-th Hermite polynomial.
My attempt is the following:
Suppose that the value of the integral is already known and that:
\begin{equation*}
\int_{-\infty}^{\infty}x^2e^{-x^2}H_n(x)H_m(x)dx=\lambda    
\end{equation*}
From the recurrence relation we have to:
\begin{equation*}
H_{n+1}(x)=2xH_{n}(x)-2nH_{n-1}(x)    
\end{equation*}
Therefore, since it is valid for any $ n $ then it will also be valid for any $ m $, that is:
\begin{equation*}
H_{m+1}(x)=2xH_m(x)-2mH_{m-1}(x)    
\end{equation*}
Solving for the terms $2xH_ {m} (x)$ and $2xH_ {n} (x)$ we have:
\begin{align}
2xH_ {n} (x) & = H_ {n + 1} (x) + 2nH_ {n-1} (x) \\
2xH_ {m} (x) & = H_ {m + 1} (x) + 2mH_ {m-1} (x)
\end{align}
And multiplying the \ textcolor {red} {\ textbf {equation 1}} by $ 2xH_m (x) $ and substituting the cleared value we have:
\begin{align*}
4x^2H_{n}(x)H_{m}(x)&=2xH_m(x)\left[H_{n+1}(x)+2nH_{n-1}(x)\right]\\
4x^2H_{n}(x)H_{m}(x)&=\left[H_{m+1}(x)+2mH_{m-1}(x)\right]\left[H_{n+1}(x)+2nH_{n-1}(x)\right]\\
4x^2H_{n}(x)H_{m}(x)&= H_{m+1}(x)H_{n+1}(x)+2nH_{m+1}(x)H_{n-1}(x)+2mH_{m-1}(x)H_{n+1}(x)+4mnH_{m-1}(x)H_{n-1}(x)\\
\end{align*}
And from what we have just found we can conclude that:
{\begin{align*}
4\lambda&=\int_{-\infty}^{\infty}e^{-x^2}\left[H_{m+1}(x)H_{n+1}(x)\right]dx+\int_{-\infty}^{\infty}e^{-x^2}2nH_{m+1}(x)H_{n-1}(x)dx+\int_{-\infty}^{\infty}2me^{-x^2}H_{m-1}(x)H_{n+1}(x)+\int_{-\infty}^{\infty}4mne^{-x^2}H_{m-1}(x)H_{n-1}(x)dx  
\end{align*}
So going back to the integral that we already had:
\begin{align*}
\cfrac{1}{4}\times4\lambda&=\cfrac{1}{4}\times \sqrt{\pi}\left\{2^{n+1}(n+1)!\left[\delta_{(m+1)(n+1)}+2m\delta_{(m-1)(n+1)}\right]+2^{n-1}(n-1)!\left[\delta_{(m+1)(n-1)}2n+4mn\delta_{(m-1)(n-1)}\right]\right\}\\
\lambda&=\cfrac{\sqrt{\pi}}{4}\left\{2^{n+1}(n+1)!\left[\delta_{(m+1)(n+1)}+2m\delta_{(m-1)(n+1)}\right]+2^{n-1}(n-1)!\left[\delta_{(m+1)(n-1)}2n+4mn\delta_{(m-1)(n-1)}\right]\right\}\\
\end{align*}
And as it was supposed from the beginning:
\begin{equation*}
\int_{-\infty}^{\infty}x^2e^{-x^2}H_n(x)H_m(x)dx=\lambda      
\end{equation*}
Then finally it will come to:
$\int_{-\infty}^{\infty}x^2e^{-x^2}H_n(x)H_m(x)dx=\cfrac{\sqrt{\pi}}{4}\left\{2^{n+1}(n+1)!\left[\delta_{(m+1)(n+1)}+2m\delta_{(m-1)(n+1)}\right]+2^{n-1}(n-1)!\left[\delta_{(m+1)(n-1)}2n+4mn\delta_{(m-1)(n-1)}\right]\right\}$
But I can't find the way to give the answer $2^{n-1}\pi^{\frac{1}{2}}(2n+1)n!\delta_{mn}+2^n\pi^{\frac{1}{2}}(n+2)!\delta_{(n+2)m}+2^{n-2}\pi^{\frac{1}{2}}n!\delta_{(n-2)m}$

I am going to work with the statistician's Hermite polynomials, which can be related to the physicists' (which appear in this question) via
$$
He_{\alpha}(x) = 2^{-\frac{\alpha}{2}}H_\alpha\left(\frac{x}{\sqrt{2}}\right)
$$
and the desired integral becomes
$$
I_{nm} = \int_{-\infty}^\infty x^2e^{-x^2}H_n(x)H_m(x)dx \\
= 2^{\frac{n+m}{2}-1}\sqrt{\pi}\int_{-\infty}^\infty x^2He_n(x)He_m(x)e^{-\frac{x^2}{2}}\frac{dx}{\sqrt{2\pi}}.
$$
We can use the identity $x^2 = He_2(x) + He_0(x)$ to split the integral into two parts
$$
I_{nm}=2^{\frac{n+m}{2}-1}\sqrt{\pi}\left(\int_{-\infty}^\infty He_n(x)He_m(x)e^{-\frac{x^2}{2}}\frac{dx}{\sqrt{2\pi}}+\int_{-\infty}^\infty He_2(x)He_n(x)He_m(x)e^{-\frac{x^2}{2}}\frac{dx}{\sqrt{2\pi}}\right)
$$
the first integral in the parentheses is equal to $n!\delta_{n,m}$ by orthogonality (of the statistician's Hermite polynomials)
$$
\int_{-\infty}^\infty He_{\alpha}(x)He_m(\beta)e^{-\frac{x^2}{2}}\frac{dx}{\sqrt{2\pi}}=n!\delta_{\alpha,\beta}
$$
and we will focus on the second
$$
J_{nm}=\int_{-\infty}^\infty He_2(x)He_n(x)He_m(x)e^{-\frac{x^2}{2}}\frac{dx}{\sqrt{2\pi}}.
$$
To calculate $J_{nm}$ we will make use of the Hermite linearization formula
$$
He_\alpha(x)He_\beta(x)=\sum_{k=0}^{\min(\alpha,\beta)}{\alpha \choose k}{\beta \choose k}k!He_{\alpha+\beta-2k}(x)
$$
which is proven in this paper (and elsewhere).
Using the linearization formula to write the product $He_2(x)He_n(x)$ in terms of a sum of individual Hermite polynomials gives
$$
J_{nm}=\sum_{k=0}^{\min(n,2)}{n \choose k}{2 \choose k}k!\int_{-\infty}^\infty He_{2+n-2k}(x)He_m(x)e^{-\frac{x^2}{2}}\frac{dx}{\sqrt{2\pi}}.
$$
There are at maximum three terms in this summation, therefore we can write out each term
$$
J_{nm}={n \choose 0}{2 \choose 0}0!\int_{-\infty}^\infty He_{n+2}(x)He_m(x)e^{-\frac{x^2}{2}}\frac{dx}{\sqrt{2\pi}}\\+{n \choose 1}{2 \choose 1}1!\int_{-\infty}^\infty He_{n}(x)He_m(x)e^{-\frac{x^2}{2}}\frac{dx}{\sqrt{2\pi}}\\+{n \choose 2}{2 \choose 2}2!\int_{-\infty}^\infty He_{n-2}(x)He_m(x)e^{-\frac{x^2}{2}}\frac{dx}{\sqrt{2\pi}}.
$$
Only the first term occurs when $n=0$, the first two when $n=1$ and all three terms otherwise.  These integrals can all be done by orthogonality with result
$$
J_{nm}=(n+2)!\delta_{n+2,m} + 2nn!\delta_{n,m} + n(n-1)(n-2)!\delta_{n-2,m}.
$$
Putting this together with the first term above leads to your sought after answer
$$
I_{nm}=\sqrt{\pi}2^{n}(n+2)!\delta_{n+2,m} + \sqrt{\pi}2^{n-1}(2n+1)n!\delta_{n,m} + \sqrt{\pi}2^{n-2}n!\delta_{n-2,m}.
$$
Note that the third term is not present if $n<2$
