What actually causes the electromagnetic field to be quantized? My understanding is that quantization in QM is because the wavefunctions are analogous to classical standing waves, where the boundary conditions determine the specific quantization.  That's how it was explained when I took physics 3 a few years ago. And, more recently, in a class on semiconductor physics last semester, I saw the derivation of the wavefunction of  an electron in a hydrogen atom and how the discrete energy levels show up.
It makes intuitive sense that the electron being "bound" by its attraction to the proton could have such an effect.  But what about light?  What are the boundary conditions that result in photons having discrete energy levels, considering they aren't electrically charged?
To be clear, I'm not asking about how we discovered light is quantized, nor am I questioning that it is.  I'm also not looking for a mathematical derivation of the quantization of the EM field, as, without the proper understanding of the relevant concepts, such derivations tend to go in one ear and out the other, IME.  Rather, I want to understand, conceptually what the mechanism for the quantization of light/the EM field is in QM and/or QFT.
 A: You're confusing two things:

*

*The discreteness or "quantization" of available physical states, and

*the quantization of excitations in those states.

The first thing exists in classical physics and does indeed arise from boundary conditions.
The second thing is truly quantum.
It's not your fault that you're confused about this.
It's confusing because books, online videos, etc. all explain it terribly, if they bother to explain it at all.
Let's check my statement that the discreteness of physical states comes from boundary conditions and exists in classical physics.
Consider a violin string, described by a function $u(x, t)$ where $u$ is the displacement of the string.
The string is clamped on both ends (that's the boundary condition) so
$$u(x=0, t) = u(x=L, t) = 0 \, .$$
As you probably know, the displacement of the string can be expressed as a Fourier series
$$u(x, t) = \sum_{k=0}^\infty c_k(t) \sin(\pi k x / L) \, .$$
There are your discrete physical states right there: they have shape $\sin(\pi k x / L)$ and are enumerated by the integer $k$.
We've written the general state of the string $u(x, t)$ as a discrete sum over these other states, each with a time dependence $c_k(t)$.
This is exactly like expressing a general electron quantum state as
$$\left \lvert \psi(x, t) \right \rangle = \sum_{k} c_k(t) \left \lvert \psi_k(x) \right \rangle$$
where the set of states $\left \lvert \psi_k \right \rangle$ are something usually chosen to have simple time dependence, i.e. something such that $c_k(t)$ are simple.
Of course, the simplest choice of $\left \lvert \psi_k \right \rangle$ are those such that $H \left \lvert \psi_k \right \rangle = E_k \left \lvert \psi_k \right \rangle$ because then $c_k(t) = \exp(-i E_k t / \hbar)$.
We do exactly the same thing in classical physics: the Fourier series is convenient precisely because once you substitute it into the wave equation that describes the string
$$\frac{\partial^2 u}{\partial t^2} = v^2 \frac{\partial^2 u}{\partial x^2}$$
you get
\begin{align}
\sum_{k=0}^\infty \ddot{c}_k(t) \sin(\pi k x / L)
&= -v^2 \sum_{k=0}^\infty c_k(t) \left( \frac{\pi k}{L} \right)^2 \sin(\pi k x / L)\\
\ddot{c}_k(t) &= - \left( \frac{\pi v k}{L} \right)^2 c_k(t) \\
\rightarrow c_k(t) &= a_k \sin\left(\pi v k t / L \right) + b_k \cos\left(\pi v k t / L \right)
\end{align}
i.e. simple sinusoidal time dependence.
This is exactly the same procedure you use to solve a quantum particle in a well or the hydrogen atom (or any physical system described by linear time invariant differential equations).
These nice shapes of the string $\sin(\pi k x / L)$ which have simple time dependence are called "normal modes".
Now here's where classical and quantum mechanics differ:
In classical mechanics, the numbers $a_k$ and $b_k$, i.e. the things that say how excited each normal mode is, are continuous numbers, whereas in quantum mechanics their squares $n = a^\dagger a$ have a discrete spectrum.
That's the quantum part of quantum mechanics.
Now I invite you to read this other question which explains why you should not think of each excitation of a field (like an electron) as being a particle with its own identity, but should rather think of electrons as identity-less units of excitation like our classical $a_k$ and $b_k$.
A: For simplicity, this answer is for the scalar field. The argument for the electromagnetic field is similar.
The energy of a Scaler field looks something like:
$$\int d^3x (\pi^2+ |\vec{\nabla \phi}| ^2 + m^2\phi^2)$$
This mathematically looks like a bunch of coupled Harmonic oscillators. You can do a Fourier transform to make it look like a bunch of decoupled Harmonic oscillators:
$$\int d^3p \omega_{\vec{p}} a^{\dagger}_{\vec{p}} a_{\vec{p}}$$ .... (1)
where $\omega_{\vec{p}}=\sqrt{|\vec{p}|^2+m^2}$
The main point is that, just like how the position of an electron can be bound to the origin using the oscillator attractive force $-m\omega^2x$, you can interpret the field equation in such a way that that the field strength at a point is bound using an "attractive force" whose mathematical form is similar to an oscillator force.
I say "attractive force" in quotes because no one would use the word "force" to describe this phenomenon. What's really happening is that the field strength tries to bounce back to 0 if it deviates from there.  This is a consequence of the field equation. You can see that the Hamiltonian I wrote in the beginning has a quadratic term $m^2\phi^2$. This term is similar to the Harmonic oscillator potential $V=x^2$. Terms like this result in an attractive force according to $F=-\frac{dV}{dx}=-x$
Now, equation $(1)$ represents the energy of an infinite number of decoupled Harmonic oscillators of a continuous range of frequences. When you quantise this, the energy of each oscillator gets discretized according to $E=n\hbar \omega_{\vec{p}}$. The total energy of the field, which is integral of the energies of each oscillator, is not discretized because the oscilators have a continuous range of frequences.
