The Hamiltonian for the Quantum Harmonic Oscillator is (disregarding constants) the Hermite operator $$ Hf = -f''+x^{2}f, $$ where $\mathcal{D}(H)$ consists of all twice absolutely continuous functions $f \in L^{2}(\mathbb{R})$ for which $Hf \in L^{2}(\mathbb{R})$.
Question: Without using properties of Hermite functions or Hermite polynomials, is there a direct method to show that
$f \in \mathcal{D}(H) \implies xf, f' \in L^{2}(\mathbb{R})$.
$f,g \in \mathcal{D}(H) \implies (Hf,g)= (f',g')+(xf,xg)=(f,Hg)$.
$H$ is selfadjoint. (2 implies the spectrum is non-negative.)
Background: Using these facts, the standard ladder argument used in Physics becomes a rigorous proof that $f$ is an $L^{2}(\mathbb{R})$ eigenfunction of $H$ with eigenvalue $\lambda$ iff
$\lambda=2n+1$ for some $n=0,1,2,3,\cdots$, and
$f$ is a constant multiple of the Hermite function $$ h_{n}(x) = (-1)^{n}e^{x^{2}/2}\frac{d^{n}}{dx^{n}}e^{-x^{2}}. $$
(Actually, only properties (1) and (2) are needed to fully justify the ladder argument.)
Summary
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