I understand why the angular momentum along a particular axis ought to be quantized to values $m\hbar$. This is because the rotation unitary operator should return to the identity after each rotation by $2\pi$ and angular momentum (along that axis), being the generator of these rotations is also then quantized accordingly.
But why is the total orbital angular momentum $L^2$ quantized? More precisely why is it quantized to values $l(l+1)\hbar^2$.
The quantization of $L_z$ (or along any direction, really) can be easily understood by making an analogy to a wave that is trapped into moving along a circle - you need to have an integer number of half-wavelengths on the circle in order for the wave to meet itself to create a standing wave.
Other answers here on physics.stackexchange do a good job explaining why $J_z$ is bounded by $J^2$, naming the generalized area of study, and links to other resources, so I won't go into those details. Instead, I want to try to build intuition by building up a physical picture similar to the waves on a ring.
What's going on is instead of the wave being confined to a ring, it's a wave that lives on the surface of a sphere. That confinement, just like the waves on a ring, leads directly to quantization of the possible waves.
Try to imagine the ways that such a sphere could vibrate. The easiest one to imagine is the whole sphere growing and shrinking. This is known as a "breathing mode" and corresponds to the $\ell = 0$, no angular momentum vibration. Next, imagine the top and bottom half of the sphere oscillating in opposite directions (top expands, bottom contracts, and vice versa). In this case, the wave is described with a Legendre polynomial evaluated at the cosine of the polar angle ($P_1(\cos \theta)$, specifically). This corresponds to the $\ell = 1,\ m=0$ wave. For the $m=\pm 1$ waves, you have to imagine a wave propagating around the equator of the sphere in one direction or the other.
When you start moving on to even higher order vibrations with more nodal lines on the sphere, you need to start being systematic. The description of the polar direction oscillations form standing waves with amplitudes described by the associated Legendre polynomials and the $z$-component is a traveling wave that comes from multiplying those polynomials by $\mathrm{e}^{i m\phi}$. When that combination is normalized to unity, by taking the absolute-square and integrating over $4\pi$ steradians, the combination becomes the spherical harmonics. Wikipedia's article on them has a nice animation that can help you visualize what I tried to describe here with words.
