Why do generalized Fourier series seem to converge to functions so far outside of the Hilbert space? My understanding is that one of the most important properties of a separable Hilbert space is that we can decompose its elements into generalized Fourier series. For example, consider a set of basis functions $f_n(x)$ that spans the Hilbert space $L^2(\mathbb{R})$. For any function $g(x) \in L^2(\mathbb{R})$, the generalized Fourier series
$$\sum_{n = 0}^\infty \langle g, f_n \rangle\, f_n(x), \qquad \qquad \langle g, f \rangle := \int_{-\infty}^\infty g(x)\, f(x)\, dx$$
converges to $g(x)$ in the $L^2$ norm. (Sometimes we can make even stronger statements; for example, for the Hilbert space $L^2([0, 2\pi])$ with the standard basis functions $e^{i n \theta}$, Carleson's theorem guarantees that the (standard) Fourier series above converges to $g(x)$ almost everywhere. As shown here, this is not true in general.)
However, the generalized Fourier series above seems to work much better than we might expect. That is, if we take a function $g(x)$ that does not lie in the Hilbert space and naively consider the above series, it appears to still converge to $g$. For example, the functions $g(x) \equiv 1$ and $g(x) = x^2$ obviously don't lie in $L^2(\mathbb{R})$. Nevertheless, when I try plotting the first few terms of the above generalized Fourier series with the Hermite basis functions, the series does indeed appear (by eye) to approach $g(x)$ more and more closely.
The above generalized Fourier series makes formal sense for any function $g(x)$ such that the inner product $\langle g, f \rangle$ is defined. For basis functions (like the Hermite functions) that lie in Schwartz space, this space of functions $g$ is much larger than $L^2(\mathbb{R})$ and even incorporates some functions that diverge exponentially at infinity. We can't apply the $L^2$ norm to functions $g(x)$ that don't lie in $L^2(\mathbb{R})$, but is there some other sense in which the above generalized Fourier series converges to $g$? If so, why?
 A: Given a regular Sturm Liouville equation
$$\frac{d}{dx} \left[ p(x) y' \right] +q(x)y+\lambda w(x)y=0$$
with $w(0)>0$, the eigenfunctions can be chosen to be an orthogonal basis with respect to the inner product
$$\langle f, g \rangle_w = \int_a^b f(x) g(x) w(x) dx$$
in the Hilbert space $H_w= \{ f | \langle f, f  \rangle_w < \infty \}$. 
The Hermite polynomials, $H_n$ are the eigenfunctions of the Hermite equation
$$\frac{d}{dx} \left[ e^{-x^2}y' \right] + \lambda e^{-x^2}y = 0$$
In this case $w= e^{-x^2}$ and the inner profuct is
$$\langle f, g \rangle_w = \int_{-\infty}^\infty f(x) g(x) e^{-x^2} dx$$
Therefore, by the above, for all $f \in H_w$ we have 
$$f=\sum_n \langle f, H'_n \rangle_w\, H'_n  \hspace{1cm}(*)$$
where $H_n'$ is the normalized Hermite Polynomial.
Now, note that 
$$\phi_n =e^{-\frac{x^2}{2}} H_n'$$
Now, let $g(x)=f(x) e^{-\frac{x^2}{2}}$. Then, from $(*)$ we get
$$g(x)= \sum_n \langle f, H'_n \rangle_w\, \psi_n$$
and 
$$\langle f, H_n' \rangle_w = \int_{-\infty}^\infty f(x) H_n' e^{-x^2} dx = \int_{-\infty}^\infty g(x) \psi_n(x) dx$$
which is the standard inner product $\langle g, \phi_n \rangle$.
So the theory of Sturm Liouville equations, applied to the Hermite equation, says that if you pick some $g(x)$, if $f(x) = e^{\frac{x^2}{2}}g(x) \in H_w$ then we have
$$g(x)= \sum_n \langle g, \psi_n \rangle \psi_n $$
and this is just the standard Fourier series convergence in the Hilbert space $\left( H_w, \langle \, \rangle_w \right)$. 
I am not sure if the integrability condition actually holds in this case, if not probably the reason is more subtle. But this seems to fit closely the general theory of eigenfuncton expansion for Sturm-Liouville equations.
P.S. I think that if one writes directly the Liouville form of the Hermite equation, then $\phi_n$ are the eigenfunctions and one gets directly 
$$$g(x)= \sum_n \langle g, \psi_n \rangle \psi_n $$
without having to go through the more general space $H_w$.
A: One possible way of thinking about this phenomenon is the following: The Hermite functions diagonalise the Fourier transform, that is $\mathcal{F}\phi_n = i^n \phi_n$ (if you choose the normalisations of the Fourier transform and the Hermite functions correctly). Therefore you get a "orthogonal decomposition" of anything has the Fourier transform as an endomorphism. This inclused $\mathcal{S}$, $L^2$ and $\mathcal{S}'$.
More precisely: Let $u\in\mathcal{S}'$ be any tempered distribution and define its Fourier coefficients as $c_n:=\langle u,\phi_n\rangle$. Using standard estimates, it follows that $c_n$ grows at most polynomially with growth rate depending on $u$ and that $\sum_n c_n\phi_n$ converges to $u$ in the $\mathcal{S}'$ sense.
Looking a bit closer one can indeed show that fast decay of $c_n$ translates to smoothness of $u$: If $u$ is in the Sobolev space $H^k(\mathbb{R})$, then $n^{k/2} c_n\in\ell^2$. This follows from partial integration and recursive formulas for $\psi_n'$. Conversely, if the coefficients decay that fast, the same argument backwards shows that $\sum c_n \phi_n$ converges to $u$ in some Sobolev norm. If $k\geq 0$, this shows in particular that it converges also in the $L^2$-norm.
Observe how this works nicely with the standard results for Fourier series for $L^2(S^1)$ where decay of the Fourier coefficients also translates to smoothness.
If one is very careful with the definitions, one can actually show that this holds for all $k\in\mathbb{R}$, not just for $k\in\mathbb{N}$. That is where the Fourier transform is useful, because the FT exchanges fractional Sobolev spaces $H^k$ with weighted $L^2$-spaces where the weight is something like $(1+|x|)^{k}$ for example.
In your example $1$ is an element of $H^{-1}$ and $x^2\in H^{-3}$ if I'm not mistaken. Therefore the series $\sum_n c_n \phi_n$ converges to $u$ in the $H^{-1}$ norm / $H^{-3}$ norm and in particular in the distributional sense.
