It's the 'accidental' degeneration of the spectrum of the Coulomb-Schrödinger operator.
For a general form of the spherical potential $V(r)$ one solves $$-R''(r) - \frac{2}{r} R' + \frac{l(l+1)}{r^2} \ R + V(r) R' = n_{r,l} \ R.$$ In the cases $V(r)=r^k, k=-1,0, 2$ there exists an additional constant of motion, the classical Runge-Lenz vector pointing along the constant perihelion vector and derived from $$p\times L = r p^2- (r\cdot p) \ p$$.
It follows, that the spectral sequences for knot lines on the unit sphere from Legendre polynomials and knots on the radius from the radial equation appear as a linear combination in the radial energy eigenvalue equation.
As a principle, the spectrum is ordered along the energy scale by the value of $n_r+l$ and the subspaces of degeneration are ordered by the eigenvalue of the Laplacian on the unit sphere $-\Delta_{S^2} = L^2$ with eigenvalues $l(l+1)$.
The fact that $l$ and not a linear form of the eigenvalue like $\sqrt{l(l+1)}$ appears in the main quantum number is again a fact of symmetry: All subspaces of states on the unit sphere with equal $L^2$, $l,m=-l\dots\l$$l,m=-l\dots l$ are degenerate and yield the same spherical energy contribution as the state with maximal $L_z=\pm l$ with no knots along meridians and kinettic energy in rotational form as opposed to the oscillating waves between the poles.
From the last sentence its clear whats the reason of the complicated representation of the symmetry group: The choice of the axis and the poles is arbitrary. A rotation of the sphere mixes all states of fixed quantum number $l$. The picture of static orbitals with fixed charge and current density has no absolute meaning; except if one breaks the symmetry by a weak magnetic field or considers the system in a fixed external lattice.