Finally I found how to do it. I post it, if someone is interested.
\begin{align}
D_{ln}(\varkappa) &= \frac{1}{\sqrt{2^nn!}\sqrt{2^ll!}}\frac{1}{\sqrt{\pi}}\int_{-\infty}^{+\infty}H_n(\tilde{x})e^{-\tilde{x}^2+\varkappa \tilde{x}}H_l(\tilde{x})\;\mathrm{d}\tilde{x} \notag\\
&= \frac{1}{\sqrt{2^nn!}\sqrt{2^ll!}}\frac{1}{\sqrt{\pi}}\int_{-\infty}^{+\infty}H_n(\tilde{x})e^{-\tilde{x}^2+\varkappa \tilde{x}-\varkappa^2/4}e^{\varkappa^2/4}H_l(\tilde{x})\;\mathrm{d}\tilde{x} \notag\\
&= \frac{1}{\sqrt{2^nn!}\sqrt{2^ll!}}\frac{1}{\sqrt{\pi}}e^{\varkappa^2/4}\underbrace{\int_{-\infty}^{+\infty}H_n(\tilde{x})e^{-(\tilde{x}-\varkappa/2)^2}H_l(\tilde{x})\;\mathrm{d}\tilde{x}}_I
\end{align}
If we pose $x = \tilde{x}-\frac{\varkappa}{2}$ in this expression, the integral $I$ becomes
\begin{equation*}
I = \int_{-\infty}^{+\infty}H_n(x+\varkappa/2)e^{-x^2}H_l(x+\varkappa/2)\;\mathrm{d}x
\end{equation*}
We know that
\begin{equation*}
H_n(x+a) = \sum_{p=0}^n \frac{n!}{(n-p)!p!}(2a)^{n-p}H_p(x)
\end{equation*}
Hence, the integral $I$ becomes
\begin{align*}
I &= \int_{-\infty}^{+\infty} \sum_{p=0}^n \frac{n!}{(n-p)!p!}\varkappa^{n-p}H_p(x) e^{-x^2} \sum_{q=0}^l \frac{l!}{(l-q)!q!}\varkappa^{l-q}H_q(x)\;\mathrm{d}x \\
&= \sum_{p=0}^n\sum_{q=0}^l \frac{n!}{(n-p)!p!}\varkappa^{n-p}\frac{l!}{(l-q)!q!}\varkappa^{l-q}\int_{-\infty}^{+\infty}H_p(x) e^{-x^2}H_q(x)\;\mathrm{d}x \\
\end{align*}
The Hermite polynomials are orthogonal in the range $(-\infty,\infty)$ with respect to the weighting function $e^{-x^2}$ and satisfy
\begin{alignat*}{2}
&&&\int_{-\infty}^{+\infty}H_p(x) e^{-x^2}H_q(x)\;\mathrm{d}x = \sqrt{\pi}2^pp!\;\delta_{pq} \\
&\Rightarrow\quad&& I = \sum_{p=0}^n\sum_{q=0}^l \frac{n!}{(n-p)!p!}\frac{l!}{(l-q)!q!}\varkappa^{n+l-p-q}\cdot \sqrt{\pi}2^pp!\;\delta_{pq}
\end{alignat*}
As this integral is nil if we have not $p=q$, we can replace the two sums by only one that goes from 0 to $\min(n,l)$. Let us say that $n<l$. Hence, the full expression for the $D$-matrix is
\begin{align}
D_{ln}(\varkappa) &= \frac{1}{\sqrt{2^nn!}\sqrt{2^ll!}}\frac{1}{\sqrt{\pi}}e^{\varkappa^2/4} \sum_{p=0}^n \frac{n!}{(n-p)!p!}\frac{l!}{(l-p)!p!}2^pp!\sqrt{\pi}\;\varkappa^{n+l-2p} \notag\\
&= \frac{\varkappa^{n+l}}{\sqrt{2^nn!}\sqrt{2^ll!}}e^{\varkappa^2/4} \sum_{p=0}^n \frac{n!}{(n-p)!p!}\frac{l!}{(l-p)!}2^p\;\varkappa^{-2p} \notag\\
&= \sqrt{\frac{n!}{l!}}\left(\frac{\varkappa}{\sqrt{2}}\right)^{n+l}e^{\varkappa^2/4} \sum_{p=0}^n \frac{l!}{(n-p)!(l-p)!p!}\left(\frac{\varkappa^2}{2}\right)^{-p} \notag\\
&= \sqrt{\frac{n!}{l!}}\left(\frac{\varkappa}{\sqrt{2}}\right)^{l-n}e^{\varkappa^2/4} \sum_{p=0}^n \frac{l!}{(n-p)!(l-p)!p!}\left(\frac{\varkappa^2}{2}\right)^{n-p}
\end{align}
Associated Laguerre polynomials $L_a^b(x)$ are given by
\begin{equation*}
L_a^b(x) = \sum_{k=0}^{a}(-1)^k \frac{(a+b)!}{(a-k)!(b+k)!k!}x^k
\end{equation*}
It suggests us to transform the expression of the $D$-matrix by posing $k=n-p$. Hence, we have
\begin{align}
D_{ln}(\varkappa) &= \sqrt{\frac{n!}{l!}}\left(\frac{\varkappa}{\sqrt{2}}\right)^{l-n}e^{\varkappa^2/4} \sum_{k=n}^0 \frac{(l)!}{(n-(n-k))!(l-(n-k))!(n-k)!}\left(\frac{\varkappa^2}{2}\right)^{n-(n-k)} \notag\\
&= \sqrt{\frac{n!}{l!}}\left(\frac{\varkappa}{\sqrt{2}}\right)^{l-n}e^{\varkappa^2/4} \sum_{k=0}^n \frac{l!}{k!(l-n+k)!(n-k)!}\left(\frac{\varkappa^2}{2}\right)^{k} \notag\\
&= \sqrt{\frac{n!}{l!}}\left(\frac{\varkappa}{\sqrt{2}}\right)^{l-n}e^{\varkappa^2/4} \sum_{k=0}^n (-1)^k\frac{([l-n]+n)!}{(n-k)!([l-n]+k)!k!}\left(-\frac{\varkappa^2}{2}\right)^{k} \notag\\
&= \sqrt{\frac{n!}{l!}}\left(\frac{\varkappa}{\sqrt{2}}\right)^{l-n}e^{\varkappa^2/4}L_n^{l-n}\left(-\frac{\varkappa^2}{2}\right)
\end{align}
It should be remembered that we had supposed that $n<l$. But that could be otherwise. In order to be general, $n_<$ and $n_>$ will be defined as $n_<=\min{(n,l)}$ and $n_>=\max{(n,l)}$ and $l-n=|l-n|$. We then have
\begin{equation}
D_{ln}(\varkappa) = \sqrt{\frac{n_<!}{n_>!}}\left(\frac{\varkappa}{\sqrt{2}}\right)^{|l-n|}L_{n_<}^{|l-n|}\left(-\frac{\varkappa^2}{2}\right)e^{\varkappa^2/4}
\end{equation}
- Method 1: recurrence relations and symmetries of the Hermite polynomials
Hint: I would use the relations
$$H_{n+1}(x)=xH_n(x)-H'_n(x), $$
$$H'_n(x)=nH_{n-1}(x) $$
and an induction argument, as follows. Let us suppose the steps $k=1,\dots,n-1$ are true; we want to show that
$$\int e^{-x^2}H_n^2(x)x^2dx=2^nn!\sqrt{\pi}(n+\frac{1}{2}). $$
Using the above relations the l.h.s is
$$\int e^{-x^2}H_n^2(x)x^2dx=\int e^{-x^2}x^2(xH_{n-1}(x)-H'_{n-1}(x))^2dx=
\int e^{-x^2}x^4H^2_{n-1}(x)dx \\
-2\int e^{-x^2}x^3H_{n-1}(x)H^{'}_{n-1}(x)dx
+\int e^{-x^2}x^2(H^{'}_{n-1}(x))^2dx=\\
\int e^{-x^2}x^4H^2_{n-1}(x)dx -2(n-1)\int e^{-x^2}x^3H_{n-1}(x)H_{n-2}(x)dx
+\\(n-1)^2\int e^{-x^2}x^2(H_{n-2}(x))^2dx. ~(*)
$$
Let us discuss $(*)$: the first 2 terms can be reduced using integration by parts on the products
$$e^{-x^2}x^4H^2_{n-1}(x)=-(-\frac{1}{2}(2x)e^{-x^2})x^3H_{n-1}(x),$$
$$e^{-x^2}x^3H_{n-1}(x)H_{n-2}(x)=-(-\frac{1}{2}(2x)e^{-x^2})x^2H_{n-1}(x)H_{n-2}(x), $$
while the last term can be evaluated by the induction hypothesis.
All we need is to remember that
$$\int e^{-x^2}x^{2q+1}H^2_{r}(x)dx=0,$$
$$\int e^{-x^2}x^{2q}H_{r}(x)H_{r-1}(x)dx=0,$$
for all $q\in\mathbb N$, $r\geq 1$ for symmetry. For the second equality we used
$H_r$ odd/even $\Leftrightarrow$ $H_{r-1}$ even/odd, which follows directly from the definition of the Hermite polynomials.
- Method 2: generating function
Using the generating function suggested by the OP, one needs to integrate w.r.t $x$ on the whole real axis both sides of
$$x^2e^{-x^2}e^{-\frac{t}{1+t}x^2}\frac{1}{\sqrt{1-t^2}}= \sum_{n=1}^\infty x^2e^{-x^2}\frac{H_n^2(x)}{2^nn!}t^n,$$
and expanding w.r.t. $t$. The integral on the l.h.s. is Gaussian and can be evaluated with standard techniques reducing to the case
$$\int_{-\infty}^\infty x^2e^{-\alpha(t)x^2}dx=\frac{1}{2}\sqrt{\frac{\pi}{\alpha(t)^3}}, ~(1)$$
for $\alpha(t)>0$ and $|t|<1$.
- Generating function with $\alpha(t)=\frac{1-t}{1+t}$
In the case under exam, multiplying both sides of the generating function identity w.r.t the product $x^2 e^{-x^2}$ ($e^{-x^2}$ is the weight for the integral identities involving the Hermite polynomials) and integrating, we arrive at (the exchange between summation and integration is allowed by the properties of Hermite poly.)
$$\frac{1}{\sqrt{1-t^2}} \int_{-\infty}^\infty x^2 e^{-x^2} e^{-\frac{t}{1+t}x^2}dx= \sum_{n=1}^\infty \int_{-\infty}^\infty x^2e^{-x^2}\frac{H_n^2(x)}{2^nn!}t^n; $$
to arrive at the thesis we write (using (1))
$$g(t):=\frac{1}{\sqrt{1-t^2}}\int_{-\infty}^\infty x^2e^{-\frac{1-t}{1+t}x^2}dx=\frac{1}{\sqrt{1-t^2}}\frac{1}{2}\sqrt{\frac{\pi}{\left(\frac{1-t}{1+t}\right)^3}}=\frac{\sqrt{\pi}}{2}\frac{1+t}{(1-t)^2}$$
for $|t|<1$ and we expand w.r.t $t$ at $t=0$ the function $g(t)$, obtaining
$$g(t)=\frac{\sqrt{\pi}}{2}\left(1+3t+\frac{1}{2!}10t^2+O(t^3)\right)$$
(the higher terms are easily computed, as well). This gives the thesis. I hope it helps.
Best Answer
While setting $t=1$ in the second-to-last equation gives a valid identity, it's not a very useful one. What you should instead do is expand the LHS of that equation in powers of $t$ just as you did already for the generating function. If you match this term-by-term to the RHS (i.e. require that the coefficient of $t^n$ be the same on both side) then you immediately obtain the desired result.