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.