Proof of the Variational Theorem I have trouble understanding the proof of the Variational Theorem. I'll recall quickly the proof to show my problems (see also this post and the answer by Mateus Sampaio for a detailed proof).
Let $H$ be the Hamiltonian operator of a system. Since it is self-adjoint, we can find an orthonormal basis $\{\Psi_n\}$ with corresponding real eigenvalues $E_n$. Suppose that $\min_n E_n$ exists and denote it by $E_0$. We may then order the eigenvalues $E_0\leq E_1\leq\dots$. We can write any normalized state $|\Psi\rangle$ as
$$
|\Psi\rangle = \sum_n c_n |\Psi_n\rangle \;\; \text{with} \;\; \sum_n |c_n|^2 = 1.
$$
We can now compute for the expected value of the energy
$$
\begin{align*}
\langle \Psi | H | \Psi \rangle &= \langle \sum_n c_n \Psi_n | H | \sum_m c_m \Psi_m\rangle\\
&= \langle \sum_n c_n \Psi_n | \sum_m c_m H \Psi_m\rangle \\
&= \sum_{m,n} c_n^* c_m \langle \Psi_n | E_m \Psi_m \rangle \\
&= \sum_n |c_n|^2 E_n \geq \sum _n |c_n|^2 E_0 = E_0.
\end{align*}
$$
What I don't understand is the second step where we pull the Hamiltonian inside the second sum over $m$ to let it act on $\Psi_m$. The Hamiltonian is linear, so interchanging it with a finite sum is not a problem. But here we interchange it with an infinite sum, so in particular we do the following:
$$
H\left(\sum_{m=0}^\infty c_m \Psi_m\right) = H\left(\lim_{M\to\infty} \sum_{m=0}^Mc_m\Psi_m\right)
= \lim_{M\to\infty} H\left(\sum_{m=0}^M c_m \Psi_m\right).
$$
But the Hamiltonian is not continuous. If it was, it would be bounded because it is linear. But it is unbounded because it involves second order derivatives.
My question is: Why can we nevertheless interchange the Hamiltonian with the limit?
 A: Yes, in general, it is not justified to exchange a limit and an unbounded operator because an unbounded operator is not continuous. However, in your case, you can apply Parseval's identity: If $\{e_n\}_{n\in\mathbb{N}}$ is an orthonormal basis of a separable Hilbert space $\mathcal{H}$, then, for any $x\in \mathcal{H}$,
$$x = \sum_{n\in \mathbb{N}}\langle e_n, x\rangle e_n.$$
Thus, if $H$ is a bounded from below self-adjoint operator (defined on $D(H) \subset \mathcal{H}$) with discrete spectrum $\sigma(H) = \{ E_0, E_1, \dots \}$, $E_0 = \min \sigma(H)$, and $\Psi_1, \Psi_2, \dots \in D(H)$ are the corresponding eigenvectors that form an orthonormal basis of $\mathcal{H}$, then, for $\Psi = \sum_{n\in\mathbb{N}} c_n \Psi_n \in D(H)$,
$$\begin{align} H\Psi &= \sum_{n\in \mathbb{N}} \langle\Psi_n, H\Psi\rangle \Psi_n \\ 
&= \sum_{n\in \mathbb{N}} \langle H\Psi_n, \Psi\rangle \Psi_n \\ 
&= \sum_{n\in \mathbb{N}} E_n \langle\Psi_n, \Psi\rangle \Psi_n \\
&= \sum_{n\in \mathbb{N}} E_n c_n \Psi_n, \end{align}$$
which is the identity you required.
As a remark, let me mention that the given assumptions on $H$ (discrete spectrum and eigenvectors form an orthonormal basis) are pretty strong and in most situations of physical interest not satisfied. A more general proof of the variational theorem employs the spectral theorem. According to the spectral theorem, for every self-adjoint operator $H$, a projection-valued measure $E$ exists such that
$$H = \int_{\sigma(H)} \lambda \ \mathrm{d}E(\lambda).$$
Hence, for every $\Psi \in D(H)$:
$$\begin{align}\langle \Psi, H \Psi\rangle = \int_{\sigma(H)} \lambda \ \mathrm{d}\langle \Psi, E(\lambda) \Psi\rangle = \int_{\sigma(H)} \lambda \ \mathrm{d}\|E(\lambda) \Psi\|^2 \\
\geq \min \sigma(H) \int_{\sigma(H)} \mathrm{d}\|E(\lambda) \Psi\|^2 = \min \sigma(H) \|\Psi\|^2. \end{align}$$
