Assume an $\ell_n$-fold degeneracy at $E_n$: then choose an orthonormal set
$\psi^{(n)}_\nu$, $\nu = 1, 2, \ldots, \ell_n$.  The function $P_R
\psi^{(n)}_\nu$ is in the same subspace:
$$
P_R \psi^{(n)}_\nu = \sum_{\kappa = 1}^{\ell_n} \psi^{(n)}_\kappa
\Gamma^{(n)}_{\kappa\nu}(R)
$$
where $\Gamma^{(n)}$ is an {\it irreducible, unitary\/} representation of
the symmetry group $G$ of the system.  Each $n$ corresponds with another
energy level.  One can purely mathematical derive irreducible
representations of a symmetry group and label the energy levels with a
quantum number this way.  A fixed choice of $\Gamma^{(n)}(R)$ defines the
base functions $\psi^{(n)}_\nu$.  This way one can also label each separate
base function with a quantum number.
\bye