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$