In Euclidean spaces of dimension n with time as 0-th coordinate, consider the representations of the Clifford algebra $$\mathcal {Cl}(1,n)$$ of the basis unit vectors $e_0, e_k, k=1,\dots n$.
The anticommutator algebra in $\mathcal {Cl}(1,n)$ for Cartesian metrics $\mathrm{diag}(1,-1,\dots,-1) $is
$$e_i \circ e_k + e_k \circ e_i = 2 g_{ik}$$
The fundamental geometric differential operator in Hilbert spaces over differentiable spin manifolds with positive metrics is the linear elliptic Dirac-Athiya operator
$$D=\sum_{i=1}^n e_i \nabla_i$$ i.e. the linear combination of the tangent vectors and th covariant derivative. In flat euclidean n-space, by the anticommutator algebra, we have
$$D^2= - \nabla^2 =-\Delta$$ a positive (second order elliptc) operator, that is invertible, such that from start values at $t=0$ in a suitable space of functions the values for later times can be calculated by linear integral operators: Fourier transform, solve for k.
Applying the Dirac operator twice we see, that all components $\psi_r$ of any solution $\psi$ of the free Dirac equation with fixed frequency $\omega$
$$ e_0 \partial_{ct} \ \psi + D \psi + \frac{mc}{\hbar} \psi = \frac{1}{c} \ i \omega \psi_\omega(x) + D \psi_\omega(x) + \frac{mc}{\hbar} \psi_\omega(x) =0 $$
are solutions of the free Klein-Gordon equation
$$\left( \frac{\omega^2}{c^2} -k^2\right)\hat{ \psi}_{\omega,r}(k) = \frac{mc^2}{\hbar^2} \hat{ \psi}_{\omega,r}(k)$$
The Dirac operator is a local mean value differential operator that is invariant with respect to coordinate transformations.
We do not try to write the Dirac operator in spherical coordinates, but for the Laplacian its easy.
The metric of the n-dimensional spherical coordinates line element
$$X=r \left( \cos \theta_{n-1}, \sin_{n-1} \cos \theta_{n-2}, \dots, \prod_k \sin \theta_k \sin \phi, \prod_k \sin \theta_k \cos \phi\right), \quad X^2=r^2, \quad e_i = \partial_{\xi_i} X$$
$$ds^2 = dr^2 + r^2 \sum_k \prod_l \sin \theta_{l<k}\ d\theta_k + \prod_l \sin \theta_l d\phi^2$$
All $ \theta_k $ vary in $(0,\pi)$, only $\phi$ is periodic in $(0,2\pi)$.
Instead of decribing the representation theory of the angular momentum tree for all dimensions, the features can already be seen for the K-G-operator for dimension 4:
$$ds^2 = dct^2 - dr^2 - r^2 (d \theta_1^2 + \sin\theta_1^2 d \theta_2^2 + \sin\theta_1^2 \sin\theta_2^2 d\phi^2 $$ such that the Laplacian decomposes by a tensor product of functions in Hilbert spaces for each single variable. The metric for the covariant derivatives is the inverse metric
$$\Delta = p_r^2 + \frac{1}{r^2} \ \left( p_\theta^2 + \frac{1}{\sin^2 \theta_1} \left( p_{\theta_1}^2 + \frac{1}{\sin^2 \theta_2} \ p_\phi^2\right) \right)$$
In contrast to the classical $\sum_k \frac{1}{g_{kk}} p_k^2 $ for covariant vectors the quantum version includes a first order term, but that does not alter the procedure.
With a product $\psi =R(r) \Theta_i(\theta_i) \dots e^{i\, m \, \phi}$
we see that each of the single variable operator produces its own discrete spectrum, depending on the eigenvalues of higher $p_{\theta_i}$.
The only periodic and cyclic variable has $p_\phi: \quad m = \pm n+\frac{1}{2}$ such that the double cover of the unit circle alang the equator produces the sign change by a rotation by $2\pi$, resulting in the Fermi-Dirac statistics of the many particles.
This fact yields the splitting of all the other eigenvalues of the diverse hypergeometric equations for $z_i = (\cos \frac{\theta}{2} \in (0,1))$, because only the eigenvalue $m^2$ of $p_{\phi}$ enters the eigenvalues for the diverse $p_\theta ^2$ operators. These are second order orbital angular momentum operators along meridians from pole to pole with polynomial solutions in $z_i$ and have no problems by rotations greater by angles outside $(0,\pi)$
For the radial equation only the total angular momentum enters, so the solution set is charcterized by by radial and total angular momentum squared.
All states of the same energy and total angular momentum are degenrated wrt to the deeper nested quantum numbers.
Bottom line: For the solution of the Dirac equation the eigensolutions of the Klein-Gordon operator can be grouped by energy, total angular momentum (half integer j) in pairs of the eigenvalues of the spin operator $$ \frac{1}{2}\ \sigma_\phi \psi = \pm \frac {1}{2} \psi$$.
Together with the pairs $l\mp+1$ eigenvalues of the orbital angular momentum operator $p_\phi$ states can be grouped by pairs $((l+1,-\frac{1}{2}),(l,\frac{1}{2}))$ in identical manner for each of the other quantum numbers from the lower dimensional sub-spheres.
This fundamental spin-orbit degeneration can only be broken by a vector field $A$ gauging the angualar momentum derivatives $p_k\to p_k-e\ A_k $ with an $A_\phi$ component producing a magnetic 2-form field in the main equatorial plane spanned by $e_r,e_\phi$.
Sorry for this lengthy explanation, it's someting like physicists Hilbert space complement of the algebraic deduction of the representation theory of $SO(n)$ as matrices .
