[proofplan]
We pass to the spectral representation of the self-adjoint operator $\mathcal H$. In that representation, $\mathcal H$ is multiplication by a real measurable function and the unitary evolution $e^{-it\mathcal H}$ is multiplication by a complex phase of modulus $1$. The modulus-one phase preserves the weighted $L^2$ condition defining the operator domain, and it cancels from the quadratic form representing the energy.
[/proofplan]
[step:Represent the Hamiltonian and the propagator as multiplication operators]
By the [spectral theorem for self-adjoint operators](/theorems/6911) (citing a result not yet in the wiki: Spectral Theorem for [Self-Adjoint Operators](/page/Self-Adjoint%20Operators)), there exist a [measure space](/page/Measure%20Space) $(X,\mathcal A,\mu)$, a real-valued measurable function $h:X\to\mathbb R$, and a unitary map
\begin{align*}
U:H\to L^2(X,\mathcal A,\mu)
\end{align*}
such that $\mathcal H=U^{-1}M_hU$, where the multiplication operator
\begin{align*}
M_h:D(M_h)\subset L^2(X,\mathcal A,\mu)\to L^2(X,\mathcal A,\mu)
\end{align*}
has domain
\begin{align*}
D(M_h)=\{f\in L^2(X,\mathcal A,\mu):hf\in L^2(X,\mathcal A,\mu)\}
\end{align*}
and satisfies $M_h f=hf$ for $f\in D(M_h)$.
The same functional calculus gives, for each $t\in\mathbb R$,
\begin{align*}
Ue^{-it\mathcal H}U^{-1}=M_{\theta_t},
\end{align*}
where $\theta_t:X\to\mathbb C$ is the measurable function $\theta_t(x)=e^{-ith(x)}$. Since $h$ is real-valued, $|\theta_t(x)|=1$ for every $x\in X$.
[/step]
[step:Use the modulus-one phase to preserve the operator domain]
Fix $t\in\mathbb R$. Define $f_0\in L^2(X,\mathcal A,\mu)$ by $f_0=U\psi_0$. Since $\psi_0\in D(\mathcal H)$ and $\mathcal H=U^{-1}M_hU$, we have $f_0\in D(M_h)$, equivalently $hf_0\in L^2(X,\mathcal A,\mu)$.
Now
\begin{align*}
U\psi(t)=Ue^{-it\mathcal H}\psi_0=M_{\theta_t}f_0=\theta_t f_0.
\end{align*}
Because $|\theta_t|=1$, we have
\begin{align*}
|h\theta_t f_0|^2=|hf_0|^2
\end{align*}
pointwise on $X$. Therefore $h\theta_t f_0\in L^2(X,\mathcal A,\mu)$, so $\theta_t f_0\in D(M_h)$. Since $D(\mathcal H)=U^{-1}D(M_h)$, it follows that $\psi(t)\in D(\mathcal H)$.
[guided]
Fix a time $t\in\mathbb R$. The domain question becomes transparent after applying the unitary spectral transform. Define
\begin{align*}
f_0=U\psi_0.
\end{align*}
This is an element of $L^2(X,\mathcal A,\mu)$ because $U:H\to L^2(X,\mathcal A,\mu)$ is unitary. More is true: since $\psi_0\in D(\mathcal H)$ and $\mathcal H=U^{-1}M_hU$, the vector $f_0$ belongs to the domain of $M_h$. By the definition of $D(M_h)$, this means exactly that
\begin{align*}
hf_0\in L^2(X,\mathcal A,\mu).
\end{align*}
The evolved vector is obtained by multiplying $f_0$ by the phase $\theta_t$. Indeed, the functional calculus identity $Ue^{-it\mathcal H}U^{-1}=M_{\theta_t}$ gives
\begin{align*}
U\psi(t)=Ue^{-it\mathcal H}\psi_0=M_{\theta_t}f_0=\theta_t f_0.
\end{align*}
To prove that $\psi(t)\in D(\mathcal H)$, it is enough to prove that $U\psi(t)\in D(M_h)$, because $D(\mathcal H)=U^{-1}D(M_h)$. Thus we must check that $h\theta_t f_0\in L^2(X,\mathcal A,\mu)$. Since $h$ is real-valued, the complex exponential $\theta_t(x)=e^{-ith(x)}$ has modulus $1$ for every $x\in X$. Hence
\begin{align*}
|h(x)\theta_t(x)f_0(x)|^2=|h(x)f_0(x)|^2
\end{align*}
for every $x\in X$ where the functions are represented. Therefore
\begin{align*}
\int_X |h\theta_t f_0|^2\,d\mu(x)=\int_X |hf_0|^2\,d\mu(x)<\infty.
\end{align*}
So $\theta_t f_0\in D(M_h)$, and consequently $\psi(t)\in D(\mathcal H)$.
[/guided]
[/step]
[step:Compute the energy in the spectral representation and cancel the phase]
Since $\psi(t)\in D(\mathcal H)$, the energy [inner product](/page/Inner%20Product) is well-defined. Using $\mathcal H=U^{-1}M_hU$, unitarity of $U$, and $U\psi(t)=\theta_t f_0$, we obtain
\begin{align*}
(\mathcal H\psi(t),\psi(t))_H=(M_h\theta_t f_0,\theta_t f_0)_{L^2(X,\mathcal A,\mu)}.
\end{align*}
By the definition of the $L^2$ inner product, linear in the first argument,
\begin{align*}
(M_h\theta_t f_0,\theta_t f_0)_{L^2(X,\mathcal A,\mu)}=\int_X h\theta_t f_0\,\overline{\theta_t f_0}\,d\mu(x).
\end{align*}
Since $|\theta_t|^2=1$, this becomes
\begin{align*}
\int_X h\theta_t f_0\,\overline{\theta_t f_0}\,d\mu(x)=\int_X h|f_0|^2\,d\mu(x).
\end{align*}
Again using unitarity of $U$ and $\mathcal H=U^{-1}M_hU$,
\begin{align*}
\int_X h|f_0|^2\,d\mu(x)=(M_hf_0,f_0)_{L^2(X,\mathcal A,\mu)}=(\mathcal H\psi_0,\psi_0)_H.
\end{align*}
Thus
\begin{align*}
(\mathcal H\psi(t),\psi(t))_H=(\mathcal H\psi_0,\psi_0)_H.
\end{align*}
Since $t\in\mathbb R$ was arbitrary, the equality holds for every $t\in\mathbb R$.
[/step]