Thanks! I just checked, and the eigenfunction expansion matches the Gaussian kernel well locally around x=z=0, but then they go to zero pretty fast before I can come up with enough higher order Hermite terms. Sort of like Taylor expansion, locally matches very well, but requires much more terms for it to work further out.

stupid question, in (11), in 1D case, when x=z goes to infinity, left side of the equation remains one, but from (12) and (13), it looks like all eigenfunctions go to zero due to the exponential term? i.e., the 'diagonal' of the kernel should be the same value, regardless how large the x&z goes, but the eigenfunctions seem to diminish when x&z goes large, how come?