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.