First establish the basis. Integrating two generating functions
against \(e^{-sx^2}\) gives
\[\begin {aligned}
&\int _{\mathbb R}e^{-sx^2}
e^{\sqrt {2s}\,xz-z^2/2}
e^{\sqrt {2s}\,xt-t^2/2}\,dx\\
&\hspace {20mm}=\sqrt {\frac \pi s}\,e^{zt}.
\end {aligned}\]
Differentiation at \(z=t=0\) proves orthogonality with precisely
the displayed normalization. Every such differentiation is
justified by a polynomial times an integrable Gaussian, uniformly
for the complex parameters in a compact set.
For completeness, suppose \(f\in L^2\) is orthogonal to every
\(h_n^s\). Since \(\operatorname {He}_n\) has degree \(n\) and leading
coefficient one, all moments of
\(g(x)=f(x)e^{-sx^2/2}\) vanish. Cauchy–Schwarz shows that
\(g\in L^1\) and that
\[F(z)=\int _{\mathbb R}g(x)e^{zx}\,dx\]
is entire: on every compact set of complex \(z\), all derivative
integrands have an integrable majorant of the form
\(|f(x)|\,|x|^k e^{-sx^2/2+R|x|}\).
All Taylor coefficients of \(F\) vanish, so \(F=0\); in particular
the Fourier transform of \(g\) vanishes. To see the needed
uniqueness directly, convolve \(g\) with a normalized Gaussian.
The Gaussian Fourier integral and Fubini express this convolution
as an integral of its zero Fourier transform, hence it is zero.
As the Gaussian variance tends to zero, these convolutions tend
to \(g\) in \(L^1\), by continuity of translations in \(L^1\) and
dominated convergence. Thus \(g=0\) almost everywhere and \(f=0\).
The orthonormal family is complete.
Now diagonalize the kernel. With
\(a=w\cosh \omega /2\), the identities
\(s=w\sinh \omega \) and \(a+s/2=w/(2q)\) give, by completing
the square in \(y\),
\[\begin {aligned}
&\int _{\mathbb R}K(x,y)e^{-sy^2/2}
e^{\sqrt {2s}\,yz-z^2/2}\,dy\\
&\hspace {8mm}=\lambda _0e^{-sx^2/2}
e^{\sqrt {2s}\,xqz-q^2z^2/2}.
\end {aligned}\]
The constant is
\(\sqrt {\pi /(a+s/2)}=\sqrt {2\pi q/w}=\lambda _0\).
Differentiating \(n\) times at zero, again under a Gaussian
majorant, gives \(Th_n^s=\lambda _0q^nh_n^s\).
If \(f=\sum _n c_nh_n^s\), boundedness and completeness give
\[Tf=\sum _{n=0}^\infty \lambda _0q^nc_nh_n^s\]
in \(L^2\). Every eigenvalue is strictly positive, which proves
positivity and injectivity. Their maximum is \(\lambda _0\).
For a positive compact operator, being trace class means that
these eigenvalues have finite sum; here the sum is
\(\lambda _0/(1-q)\). This proves all operator assertions.
It also gives the Mehler kernel expansion
\[K(x,y)=\sum _{n=0}^\infty
\lambda _0q^nh_n^s(x)h_n^s(y)\]
in \(L^2(\mathbb R^2)\): the squared Hilbert–Schmidt norm of the
tail is \(\lambda _0^2\sum _{n>r}q^{2n}\) after truncation at \(r\).
This convergence statement suffices for the operator identity.