Are the Hermite-Gauss functions linearly dense in $L^1(\mathbb{R})$?
Hermite functions are dense in $L^1[\mathbb{R}]$, but we cannot use coefficients computed as above. Every $L^1[\mathbb{R}]$ function can be approximated in $L^1[\mathbb{R}]$ by continuous bounded functions of compact support - so they are also in $L^2[\mathbb{R}]$ and can be approximated by linear combinations of Hermite functions.
The functions $g_m(t) = H_m(t)\,e^{-t^2/2} $ just give an orthogonal base of $L^2(\mathbb{R})$, since: $$ \int_{-\infty}^{+\infty} f_m(t)^2\,dt = 2^m m! \sqrt{\pi}. $$ An orthonormal base is given by: $$ f_m(x) = \frac{1}{\sqrt{2^n n! \sqrt{\pi}}}\,H_m(x)\, e^{-x^2/2}. $$ We may notice that if $m$ is odd then $\int_\mathbb{R}f_m(x)\,dx = 0$, while: $$ \int_{\mathbb{R}} f_{2n}(x)\,dx = \frac{(2n)!}{n!}\sqrt{2\pi}\frac{1}{\sqrt{4^n (2n)! \sqrt{\pi}}}=\sqrt{\frac{2\sqrt{\pi}}{4^n}\binom{2n}{n}}\approx \sqrt{2}\cdot n^{-1/4}.$$ Moreover, we have $\left|\,f_m(x)\,\right|\leq\frac{1}{2}$ for every $m$ and every $x$, and $f_m$ is essentially zero outside $\left[-\frac{\pi}{2}\sqrt{m},\frac{\pi}{2}\sqrt{m}\right]$. So we know the behaviour of our base with respect to $L^1,L^2,L^\infty$. By interpolation or other techniques, it is not difficult to check that linear combinations of $f_n$s cannot be dense in $L^1$.
For instance, we may consider the function $g(x)=\frac{1}{\sqrt{x}}\cdot\mathbb{1}_{(0,1)}(x)$ that belongs to $L^1\setminus L^2$.
Every $f_n$ has a rather large support, so if we compute: $$ a_n = \int_{0}^{1} g(x)\,f_n(x)\,dx $$ we have that: $$ g_N(x)=\sum_{n=0}^{N} a_n\,f_n(x) $$ is a not-so-bad approximation of $g(x)$ over $[0,1]$, but it has a long tail, such that $g_N(x)$ does not converge to $g(x)$ in $L^1(\mathbb{R})$. To make this argument visually appealing, this is the situation for $N=10$:
$\hspace0.5in$
This non-density argument is even more convincing (at least to me) if we notice that every $g_N$ is an entire function (of order 2) while $g$ is a compact-supported function with a branch-like singularity in a right neighbourhood of the origin. So the fact that $L^2$ and the Schwarz space are mapped into theirselves by the Fourier transform is really crucial for proving the density of the span of the "Hermite base".
The ordinary heat equation is $$ \frac{\partial F}{\partial t}=\frac{\partial^{2}F}{\partial x^{2}},\\ F(0,x)=f(x). $$ The initial heat distribution is $f$. The ordinary heat kernel is the Guassian, and the resulting time evolution operator $T(t)=e^{tL}$ is a constractive $C_0$ semigroup on every $L^{p}(\mathbb{R})$ for $1 \le p < \infty$; the fact that it is $C_{0}$ gives $$ \|e^{tL}f-f\|_{p}\rightarrow 0 \mbox{ as } t\downarrow 0,\;\; 1 \le p < \infty. $$ That gives a nice approximation. In fact, this approximation technique goes back to Weierstrass in his original proof of the Weierstrass Approximation Theorem.
The Hermite functions $h_{n}(x)=H_{n}(x)e^{-x^{2}/2}$ are the $L^{2}$ eigenfunctions of $$ Lf = -\frac{d^{2}f}{dx^{2}}+x^{2}f $$ with eigenvalues $\lambda = 2n+1$ for $n=0,1,2,3,\cdots$. The Hermite functions $\{ h_{n} \}_{n=0}^{\infty}$ form a complete orthonormal basis of $L^{2}(\mathbb{R})$. The heat equation associated with $L$ is $$ \frac{\partial F}{\partial t}=\frac{\partial^{2}F}{\partial^{2}x}-x^{2}F,\\ F(0,x)=f(x). $$ This heat equation is better behaved in many ways than the ordinary heat equation because $-x^{2}F$ pulls heat out of the system near $\pm\infty$ for positive $F$. The time evolution solution operator $T(t)=e^{tL}$ in this case is $$ T(t)f = \sum_{n=0}^{\infty}e^{-(2n+1)t}(f,h_n)h_n. $$ So the approximation problem by Hermite functions is closely related to the continuity properties of the heat solution at $t=0$. That's why one studies the Hermite kernel function $$ K(r,x,y)=\sum_{n=0}^{\infty}r^{n}h_{n}(x)h_{n}(y), $$ which has an explicit representation as a bivariate Gaussian: $$ K(r,x,y) = \frac{1}{\sqrt{\pi(1-r^{2})}} \exp\left\{-\frac{1}{4}\frac{1-r}{1+r}(x+y)^{2}-\frac{1}{4}\frac{1+r}{1-r}(x-y)^{2}\right\}. $$ So the approximation problem can be studied by looking at the question $$ f\;\; ? = ?\;\; \lim_{r\downarrow 0}\int_{-\infty}^{\infty}K(r,x,y)f(y)dy = \lim_{r\downarrow 0}\sum_{n=0}^{\infty}r^{n}(f,h_n)h_n(x). $$