Are polynomials dense in Gaussian Sobolev space? Let $\mu$ be standard Gaussian measure on $\mathbb{R}^n$, i.e. $d\mu = (2\pi)^{-n/2} e^{-|x|^2/2} dx$, and define the Gaussian Sobolev space $H^1(\mu)$ to be the completion of $C_c^\infty(\mathbb{R}^n)$ under the inner product
$$\langle f,g \rangle_{H^1(\mu)} := \int f g\, d\mu + \int \nabla f \cdot \nabla g\, d\mu.$$
It is easy to see that polynomials are in $H^1(\mu)$.  Do they form a dense set?
I am quite sure the answer must be yes, but can't find or construct a proof in general.  I do have a proof for $n=1$, which I can post if anyone wants.  It may be useful to know that the polynomials are dense in $L^2(\mu)$.
Edit: Here is a proof for $n=1$.
It is sufficient to show that any $f \in C^\infty_c(\mathbb{R})$ can be approximated by polynomials.  We know polynomials are dense in $L^2(\mu)$, so choose a sequence of polynomials $q_n \to f'$ in $L^2(\mu)$.  Set $p_n(x) = \int_0^x q_n(y)\,dy + f(0)$; $p_n$ is also a polynomial.  By construction we have $p_n' \to f'$ in $L^2(\mu)$; it remains to show $p_n \to f$ in $L^2(\mu)$.  Now we have
$$ \begin{align*} \int_0^\infty |p_n(x) - f(x)|^2 e^{-x^2/2} dx &= \int_0^\infty \left(\int_0^x (q_n(y) - f'(y)) dy \right)^2 e^{-x^2/2} dx \\
&\le \int_0^\infty \int_0^x (q_n(y) - f'(y))^2\,dy \,x e^{-x^2/2} dx \\
&= \int_0^\infty (q_n(x) - f'(x))^2 e^{-x^2/2} dx \to 0 \end{align*}$$
where we used Cauchy-Schwarz in the second line and integration by parts in the third.  The $\int_{-\infty}^0$ term can be handled the same with appropriate minus signs.
The problem with $n > 1$ is I don't see how to use the fundamental theorem of calculus in the same way.
 A: I think I have a different proof. Let $\gamma$ be the Gaussian measure, that is, $\gamma$ is given by the Radon-Nikodym density,
$$\mathrm{d}\gamma(x) = \frac{\mathrm{e}^{-x^2}}{\sqrt{\pi}} \mathrm{d}x.$$
Also, consider the Ornstein-Uhlenbeck operator given as,
$$L := -\frac12 \Delta + x \cdot \nabla.$$
We can verify that for $u$ and $v$ in $C_{\mathrm{c}}^\infty(\mathbf R^d)$ we have the symmetricity
$$\int_{\mathbf{R}^d} u L v \, \gamma(\mathrm{d}x) = \int_{\mathbf{R}^d} \nabla u \cdot \nabla v \, \gamma(\textrm{d}x).$$
The nice thing about this is that the Hermite polynomials are an orthogonal basis for the Gaussian Hilbert space, that is, $L^2(\gamma)$. These can be define using the Rodrigues' formula, that is,
$$H_n(x) = (-1)^n \mathrm{e}^{x^2} \partial_x^n \mathrm{e}^{-x^2}.$$
Furthermore, we have,
$$L H_n = n H_n.$$
Also, we have,
$$H_n' = 2n H_{n - 1}.$$
Proving that the polynomials form a basis for $L^2(\gamma)$ is not hard, just consider the entire function
$$F(z) = \int_{-\infty}^\infty \mathrm e^{zx - x^2} \, \frac{\mathrm{d}x}{\sqrt\pi}.$$
Hence, so are the Hermite polynomials as they are linear combinations of the monomials $(x \mapsto x^n)_n$. Higher-order Hermite polynomials are the canonical tensor extensions. So, given an $f$ in $L^2$ we have that $f$ is given in the form
$$f = \lim_N f_N = \lim_N \sum_{n = 0}^N a_n \frac{a_n}{\sqrt{n! 2^n}}.$$
And $f_N$ gives the limit of functions that converges to $f$.
Also, recall the orthogonality (after scaling),
$$\int_{\mathbf{R}^d} h_n h_m \, \gamma(\mathrm{d}x) = \delta_{nm}.$$
where
$$h_n = \frac{H_n}{\sqrt{n! 2^n}}.$$
We can rewrite the inner product on the Sobolev space as
$$\langle u, v \rangle_{H^1} = \langle u, (L + 1) v \rangle_{L^2(\gamma)}.$$
We only have to care about the bilinear form
$$\mathcal E(u, v) = \int_{\mathbf{R}^d} u L v \, \gamma(\mathrm{d}x).$$
Or so, picking $u = v = f_N - f$ we have after noting that
$$g_N := f - f_N = \sum_{n = N + 1}^\infty a_n \frac{H_n}{\sqrt{n! 2^n}},$$
$$
\begin{align}
\mathcal E(u, v) &= \int_{\mathbf{R}^d} f_N L f_N \, \gamma(\mathrm{d}x)\\
&\quad + \int_{\mathbf{R}^d} f L f \, \gamma(\mathrm{d}x)\\
&\quad - 2\int_{\mathbf{R}^d} f_N L f \, \gamma(\mathrm{d}x)\\
&= \sum_{n = 0}^\infty n \frac{|a_n|^2}{n! 2^n} - \sum_{n = 0}^N n \frac{|a_n|^2}{n! 2^n}\\
&= \sum_{n = N + 1}^\infty n \frac{|a_n|^2}{n! 2^n}.
\end{align}.
$$
And as $\sum |a_n|^2$ converges due to the $L^2$ density, so should this.
I hope I did not make any mistakes, just occurred to me while I was biking...
A: Nate, I once needed this result, so I proved it in  Dirichlet forms with polynomial domain (Math. Japonica 37 (1992) 1015-1024). There may be better proofs out there, but you could start with this paper. 
A: Byron's paper, which he linked in his (accepted) answer, has a proof in a more general setting (where the Gaussian measure can be replaced by any measure with exponentially decaying tails).  Here is a specialization of it to the Gaussian case, which I wrote up to include in some lecture notes.  I guess I was on the right track, the trick was to differentiate more times.
$\newcommand{\R}{\mathbb{R}}$
  For continuous $\psi : \R^k \to \R$, let
  \begin{equation*}
    I_i \psi(x_1, \dots, x_k) = \int_0^{x_i} \psi(x_1, \dots, x_{i-1},
    y, x_{i+1}, \dots, x_k)dy.
  \end{equation*}
  By Fubini's theorem, all operators $I_i, 1 \le i \le k$ commute.  If
  $\psi \in L^2(\mu)$ is continuous, then $I_i \psi$ is also
  continuous, and $\partial_i I_i \psi = \psi$.  Moreover,
  $$\begin{align*}
    \int_{0}^\infty |I_i \psi (x_1, \dots, x_k)|^2 e^{-x_i^2/2} dx_i &=
    \int_{0}^\infty \big\lvert\int_0^{x_i} \psi(\dots, y,\dots)\,dy\big\rvert^2
    e^{-x_i^2/2} dx_i \\
    &\le \int_0^\infty \int_0^{x_i} |\psi(\dots, y, \dots)|^2 \,dy x_i
    e^{-x_i^2/2}\,dx_i && \text{Cauchy--Schwarz} \\
    &= \int_0^\infty |\psi(\dots, x_i, \dots)|^2 e^{-x_i^2/2} dx_i
  \end{align*}$$
  where in the last line we integrated by parts.  We can make the same
  argument for the integral from $-\infty$ to $0$, adjusting signs as
  needed, so we have
  \begin{equation*}
    \int_\R |I_i \psi(x)|^2 e^{-x_i^2/2} dx_i \le \int_\R |\psi_i(x)|^2
    e^{-x_i^2/2} dx_i.
  \end{equation*}
  Integrating out the remaining $x_j$ with respect to $e^{-x_j^2/2}$
  shows
  $$
    ||{I_i \psi}||_{L^2(\mu)}^2 \le ||{\psi}||_{L^2(\mu)}^2,
  $$
  i.e. $I_i$ is a contraction on $L^2(\mu)$.
Now for $\phi \in C_c^\infty(\R^k)$, we can approximate $\partial_1
  \dots \partial_k \phi$ in $L^2(\mu)$ norm by polynomials $q_n$.
  If we let $p_n = I_1 \dots I_k q_n$, then $p_n$ is again a
  polynomial, and $p_n \to I_1 \dots I_k \partial_1 \dots \partial_k
  \phi = \phi$ in $L^2(\mu)$.  Moreover,
  $\partial_i p_n = I_1 \dots I_{i-1} I_{i+1} \dots I_k q_n \to I_1
  \dots I_{i-1} I_{i+1} \dots I_k \partial_1 \dots \partial_k \phi =
  \partial_i \phi$ in $L^2(\mu)$ also.
