$\ell$, i.e. $ H_{n \ell,n^\prime \ell^\prime} = \delta_{\ell, \ell^\prime} H_{n,n^\prime}$, and the eigenvalue problem [*for a given $\ell$ only*]{} needs be solved. While a diagonalization is now required for each $\ell$ separately, this greatly reduces the computational cost. Evaluating the last two inner products of Eq. (\[hmatrix2\]), we obtain all the matrix elements for the Hamiltonian, $$H_{n,n^\prime}=\delta_{n,n^\prime}\left[\frac{\hbar^2\beta_{n,\ell}^2}{2ma^2}+\frac{\hbar\omega_c\ell}{ {2}}\right]+\frac{1}{8}m\omega_c^2a^2\rho^2_{n,n^\prime}.
\label{hmatrix3}$$ Using energy units of $\hbar^2/(2ma^2)$, we can write this in dimensionless form $$\frac{H_{n,n^\prime}}{\hbar^2/(2ma^2)} \equiv h_{n,n^\prime} = \delta_{n,n^\prime}\left[\beta_{n,\ell}^2+2{\gamma} \ell \right]+
    {\gamma^2}\rho^2_{n,n^\prime}
    \label{hmatrix4}$$ with $$\rho^2_{n,n^\prime}=\frac{\bra{n}r^2\ket{n^\prime}}{a^2}=\frac{2\int_{0}^1J_\ell(\beta_{n^\prime,\ell}\rho)\rho^2J_\ell(\beta_{n,\ell}\rho)\rho d\rho}{J_{\ell+1}(\beta_{n^\prime,\ell})J_{\ell+1}(\beta_{n,\ell})}.
    \label{rho_defn}$$\
Note that we have suppressed the dependency on $\ell$ in the matrix labels, but it is carried as a parameter in all of these matrices. Recall that the constant $G$ is related to the ratio of flux to the quantum of flux (given by $\gamma$), through $$\gamma = \pi G = \frac{B\pi a^2}{h/e}=\frac{\Phi}{\Phi_0}.
    \label{gamma_defn}$$ As was argued earlier, we expect $\gamma$ to roughly determine the degeneracy in this system. The integral in Eq. (\[hmatrix\]) can be evaluated analytically; we get $$\rho^2_{n,n^\prime} = \delta_{n,n^\prime} \Bigl[ \frac{\beta^2_{n,\ell} + 2\ell^2 -2} {3 \beta^2_{n,\ell}} \Bigr] + (1 - \delta_{n,n^\prime}) \Bigl[\frac{8\beta_{n^\prime,\ell}\beta_{n,\ell}}{(\beta_{n^\prime,\ell}^2-\beta_{n,\ell}^2)^2} \frac{J_{\ell-1}(\beta_{n^\prime,\ell})J_{\ell-1}(\beta_{n,\ell})}{J_{\ell+1}(\beta_{n^\prime,\ell})J_{\ell+1}(\beta_{n,\ell})} \Bigr]
\label{bess_int1}$$ We now have all the matrix elements analytically, to insert into Eq. (\[matrix1\]) for diagonalization. The big 