- The revolutionary insight here is precisely that the prediction from classical electrodynamics - which would suggest that the allowed energies for each mode are continuous - is wrong. If you instead restrict the allowed energies to discrete multiples of $h\nu$ (where $\nu$ is the frequency of the mode in question and $h$ is some constant), then the resulting prediction matches experiment. Einstein suggested that we take this apparently ad-hoc discretization seriously, which was one of the major stepping stones toward quantum mechanics.
- Standing waves and plane waves are not mutually exclusive phenomena. $\mathbf E = \sin(kx)\sin(\omega t) \hat y$ is a plane wave (the wavefronts are planes $x=const.$) and also a standing wave.
- The two quadratic degrees of freedom in the Hamiltonian are not quite $\mathbf E$ and $\mathbf B$. If you formulate electromagnetism in the Hamiltonian framework, then (choosing the Coulomb gauge as per knzhou's nice answer here) the Hamiltonian can be expressed in terms of $\mathbf E$ and the vector potential $\mathbf A$ as $H \sim \frac{1}{2}(E^2+\omega^2A^2)$. It turns out that $\mathbf E$ is the canonically conjugate "momentum" to $\mathbf A$, and it is in this sense that the electromagnetic field is analogous to a harmonic oscillator.
- This is answered by $(3)$.
The electromagnetic field Lagrangian density is given by
$$\mathcal L = \frac{B^2}{2\mu_0}-\frac{\epsilon_0E^2}{2}$$
Writing $\mathbf E = -\nabla \phi - \dot{\mathbf A}$ and $\mathbf B = \nabla \times \mathbf A$, we can express this in terms of the vector potential as
$$\mathcal L = \frac{1}{2\mu_0}(\nabla\times\mathbf A)^2-\frac{\epsilon_0}{2}\left(-\nabla \phi - \dot{\mathbf A}\right)^2 $$
Treating $\mathbf A$ like our generalized coordinates, we define the canonical momenta $\boldsymbol \Pi$ to be
$$\boldsymbol \Pi = \frac{\partial \mathcal L}{\partial \dot{\mathbf A}} = \epsilon_0 (-\nabla\phi - \dot{\mathbf A})$$
$$\iff \dot{\mathbf A} = -\frac{1}{\epsilon_0}\boldsymbol\Pi + \nabla \phi$$
Performing the requisite Legendre transformation to move to the Hamiltonian formalism yields the Hamiltonian density:
$$\mathcal H = \boldsymbol \Pi \cdot \dot{\mathbf A} - \mathcal L = \frac{1}{2\epsilon_0}\Pi^2 + \frac{1}{2\mu_0} (\nabla\times \mathbf A)^2 + \boldsymbol{\Pi} \cdot \nabla \phi$$
Once we integrate over space, we'll be able to integrate the last term by parts. If the fields go to zero at infinity, it simply vanishes in the absence of charges, yielding
$$\mathcal H = \frac{1}{2\epsilon_0}\Pi^2+ \frac{1}{2\mu_0}(\nabla \times \mathbf A)^2$$
The last thing is to rewrite $(\nabla \times \mathbf A)^2$. This is somewhat tedious in general, but the idea is that if we're looking at a single electromagnetic wave mode in a cavity, then $\mathbf A$ will be a plane wave of the form $\mathbf A_0 e^{i\mathbf k \cdot \mathbf r - \omega t}$, where $\omega = ck$ and $c$ is the speed of light. If that's the case, then $\nabla \times \mathbf A = i\mathbf k \times \mathbf A$. If we further work in the Coulomb gauge, then $\nabla \cdot \mathbf A =i \mathbf k\cdot \mathbf A = 0$; this allows us to write (after some algebra)
$$(\nabla \times \mathbf A) = k^2 A^2$$
and so
$$\mathcal H = \frac{1}{2\epsilon_0}\Pi^2 + \frac{k^2}{2\mu_0} A^2$$
Looking back, you should identify $\boldsymbol \Pi = \epsilon_0 \mathbf E$. Factoring out the $\epsilon_0$ and noting that $\frac{1}{\epsilon_0\mu_0}=c^2$, we have finally that
$$\mathcal H = \epsilon_0 \left(\frac{1}{2}E^2+ \frac{1}{2}\omega^2 A^2\right)$$
