Essential Selfadjointness of Quantum Harmonic Oscillator Hamiltonian

Summary

1. The symmetry (and hence self-adjointness) of the maximal operator is a case of the so-called Sears Theorem; the presentation of the in the book of F. A. Berezin and M. Shubin [BerShu91] (Section 2.1, starting on page 50) should handle the boundary issues to your satisfaction. (In particular, they only work on $[-T, T]$ until they're ready to handle the limit $T \to \infty$, and use clever integrating factors roughly akin to $\left( 1 - \frac{|x|}{T}\right)$ to force the necessary cancellation at the endpoints, when needed.)
2. Once that is established, the compatibility with the quadratic form $\langle f', g' \rangle + \langle xf, xg \rangle$ arises by the theory of closed, densely defined, semibounded-below quadratic forms (see Reed and Simon [ReeSim72], Volume I, Chapter VIII [I think? I don't have it on me]; in the alternative, see T. Kato's text [Kat76] on the more general closed sectorial forms [Chapter VI]).
3. The general family of problems (in one dimension) is the study of limit-point and limit-circle operators; see D.B. Sears's paper [Sea50], and Volume II of Reed and Simon [ReeSim75], Chapter X, Appendix 1 for more information.

Self-adjointness -- Sears Theorem

More generally, for the (symmetric, densely-defined) operator (and using $\mathfrak{D}$ for domains)

\begin{align} \mathfrak{D}(H_0) & = C^{\infty}_{0}(\mathbb{R}) = \lbrace C^{\infty}(\mathbb{R})\text{ functions with compact support} \rbrace\\ H_0 y &= -y'' + v(x) y, \end{align} where $v(x)$ is real-valued, measurable, locally-bounded, and not too negative, a proof of the symmetry of $H_0^*$ (i.e., a proof of the essential self-adjointness of $H_0$) is given in Berezin and Shubin (Section 2.1, starting p. 50). They call the result the Sears Theorem, which I believe refers to D. B. Sears and his paper [Sea50], handling a similar problem on $(0, \infty)$. In the particular case $v(x) = x^2$, I believe that their $H_0^*$ is your $H$, so that you have the necessary self-adjointness ($H_0$ is essentially self-adjoint if and only if $\overline{H_0} = H_0^* = H_0^{**}$.)

Moreover, to avoid special notation, we note that if $v(x) \geq -1$ for all $x$, they show that for all $f \in \mathfrak{D}(H_0^*)$, $f'(x) \in L^2(\mathbb{R})$.

The proof is long enough that I don't think I can summarize it with justice here. For example, one of the steps you are effectively left to do for yourself is that for $f \in \mathfrak{D}(H_0^*)$, the following expression is finite (if $f$-dependent): $$ \sup_{T > 0} \frac{1}{T^2} \left( \int_{-T}^T |f(x)|^2 \, dx - 2 |f(0)| T \right). $$

Quadratic-Form Compatibility -- Lots of Theory

The work of Berezin and Shubin does not directly show $xf(x) \in L^2(\mathbb{R})$, or compatibility with the semi-bounded-below (hence sectorial), closed quadratic form \begin{align} \mathfrak{D}(\mathfrak{h}) &= \lbrace f \in L^2(\mathbb{R}): f'(x) \in L^2(\mathbb{R}), xf(x) \in L^2(\mathbb{R}) \rbrace\\ \mathfrak{h}(f, g) &= \langle f', g' \rangle + \langle xf, xg \rangle \end{align} (I leave the proof of closure to you unless requested.)

Fortunately, there is a workaround for this with the theory of quadratic forms and some symbol-pushing. Define $H_{\text{sum}}$ to be the literal operator sum of $-D^2$ and the multiplication-by-$x^2$ operator, \begin{align} \mathfrak{D}(H_{\text{sum}}) & = \lbrace f \in L^2(\mathbb{R}): f''(x) \in L^2(\mathbb{R}), x^2f(x) \in L^2(\mathbb{R}) \rbrace\\ H_{\text{sum}}y & = -y'' + x^2y \end{align} and define $H_{\text{quad}}$ to be the nonnegative self-adjoint operator coming by theory from the nonnegative, closed, quadratic form $\mathfrak{h}$ (see, e.g., [Kat76], Chapter VI, Section 2.1, Thm. 2.1, pp. 322-323), which in particular satisfies \begin{equation} \tag{$\dagger$} \label{eq:domincl} \mathfrak{D}(H_{\text{quad}}) \subseteq \mathfrak{D}(\mathfrak{h}), \end{equation} and for all $f \in \mathfrak{D}(H_{\text{quad}})$ and $g \in \mathfrak{D}(\mathfrak{h})$, \begin{equation} \tag{$\dagger \dagger$} \label{eq:compatible} \langle H_{\text{quad}}f, g \rangle = \mathfrak{h}(f, g). \end{equation}

We have, or can show, the inclusions $$H_0 \subseteq H_{\text{sum}} \subseteq H_{\text{quad}},$$ and hence by taking adjoints $$H_{\text{quad}}^* \subseteq H_{\text{sum}}^* \subseteq H_0^*.$$ Yet $H_{\text{quad}}$ is self-adjoint and closed, and $H_0$ is essentially self-adjoint, and so by combining the above two lines, $$H_0 \subseteq H_{\text{quad}} \subseteq H_0^* = \overline{H_0}$$,

and hence $H_{\text{quad}}$ is a closed extension of $H_0$ inside the closure $\overline{H_0}$, so $H_{\text{quad}} = H_0^* = \overline{H_0} = H$. Therefore, the properties \eqref{eq:domincl} and \eqref{eq:compatible} apply to $H$, so $xf(x) \in L^2(\mathbb{R})$ by the former, and the compatibility with the quadratic form by the latter.

(In fact, $H_{\text{sum}}$ is equal to $H$ as well, but I don't know how to prove that without involving the Hermite functions.)

References:

[BerShu91] F. A. Berezin and M. Shubin, The Schrödinger Equation, Mathematics and its Applications vol. 76, Springer, 1991. Springer link

[Kat76] T. Kato, Perturbation Theory for Linear Operators, Second Edition, Grundlehren der mathematischen Wissenschaftern 132, Springer-Verlag, 1976.

[ReeSim72] Michael Reed and Barry Simon, Methods of Mathematical Physics, vol I: Functional Analysis, Academic Press, 1972.

[ReeSim75] Michael Reed and Barry Simon, Methods of Mathematical Physics, vol II: Fourier Analysis, Self-Adjointness, Academic Press, 1975.

[Sea50] D. B. Sears, "Note on the uniqueness of Green's functions associated with certain differential equations," Canad. J. Math, 2, 314 -- 325. Scopus link