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

Firstly I just wanted to clarify what you meant by "twice absolutely continuous", my assumption is that means "twice absolutely continuously differentiable". If that is the case, then the fact that the question is posed on $\Bbb R$ is important, of course $f\in L^2(\Bbb R)$ is a strong assumption.

Assume $Hf=-f''+x^2f\in L^2(\Bbb R)$, this tells us there is some $h\in L^2$ such that $Hf=h$, and so we obtain $$(Hf,f)=\int_\Bbb R-f''f+x^2f^2=\int_\Bbb R(f')^2+x^2f^2$$ $$=(h,f)\le\|h\|_{L^2(\Bbb R)}\|f\|_{L^2(\Bbb R)}<\infty,$$ note that by the comment of @mgn $f$ and $f'$ vanish at infinity, so the integration by parts is justified. This gives us $(1)$.

Take $f,g\in\mathcal{D}(H)$, and we have $$(Hf,g)=\int_\Bbb R-f''g+x^2fg=\int_\Bbb Rf'g'+x^2fg=(f',g')+(xf,xg)$$ $$=\int_\Bbb R-fg''+x^2fg=(f,Hg)$$ again the integration by parts is justified by the decay of $f,f',g,g'$ at infinity. So we have $(2)$.

If we posed the question on $(-1,1)$ without imposing conditions on the end points we will have problems, since we cannot justify the integration by parts, we could achieve the same "problem" by setting $\mathcal{D}(H)$ to be the set of functions $f\in L^2_{loc}(\Bbb R)$, such that $Hf\in L^2(\Bbb R)$ since then we cannot justify the integration by parts.

Conversly, the problem may be well posed on a finite interval if we set zero boundary conditions at the two end points, then I imagine that the twice absolutely continuous assumption means that $f',g'$ are at least bounded, so we can justify the integration by parts. Seeking $f\in L^2(\Bbb R)$ is a similar analogy to this, as we are setting a condition at the "endpoints" $\pm\infty$ (really it's the limiting case).