One of the problem (or difficulty) is the conversion between cylindrical quantum number $|Nn_zm\rangle$ and spherical quantum number $|Nlm\rangle$. In the limit of $\delta = 0$, all state with same $N$ has same energy and the conversion is not necessary.

With out LS coupling, the energy level can be found by the approximation

$\displaystyle E(N, n_z) = \hbar\omega \left( \frac{3}{2} + \left(1 + \frac{\delta}{3} \right) N - \delta n_z \right)$

We can see from the formula, for $\delta>0$, the largest $n_z$ has lowest energy. And the conversion to spherical harmonics can be done using projection method.

With LS coupling, the conversion becomes complicated. One issue is the coupling with the spin. Since the total angular momentum $J = L + S$ is not a conservative quantity due to deformation. In the body-fixed frame, the additional good quantum number is the z-projection of $J$, which denote as $K = |m_l \pm 1/2 |$, the quantum state becomes

$|Nn_z m_l K\rangle$

Note that each state will have only 2 degeneracy with negative $K$. On the spherical basis, the coupling with $J$ is a standard textbook content. the eigen state is

$|N l j mj \rangle = |N l j K\rangle$

The conversion get complicated, because there is no clear rule on how $| Nn_z m_l K\rangle \rightarrow |Nl j m_j \rangle$. When we want to calculate the energy level using cylindrical basis with the LS coupling, since the $J$ is not clear, even through we can write down the formula

$\displaystyle E( N, n_z, l, j) = \hbar \omega_z \left(\frac{1}{2} + n_z \right) + \hbar \omega_\rho \left(1 + N - n_z \right) \\ + a \frac{1}{2} ( j(j+1) - l(l+1) - s(s+1)) + b l^2$

It is difficult to know the $j , l$.

The go around this conversion, we can diagonal the Hamiltonian using spherical basis. The Hamiltonian is can be written as

$\displaystyle H = - \frac{\hbar^2}{2m} \nabla^2 + \frac{m\omega}{2}r^2 \left( 1 - \frac{2}{3}\sqrt{\frac{16\pi}{5}} Y_{20}(\theta, \phi) \delta \right) + a L\cdot S + b L^2$

The first 2 terms is the spherical harmonic oscillator. The energy is $\hbar\omega( 3/2 + N)$. The $Y_20$ restricts the coupling between $l, l \pm 2$ (see here.)

For example, when using the basis on $N = 1$, which are

$\displaystyle |N l j m_j \rangle = \left( |1 1 \frac{3}{2} \frac{3}{2}\rangle , |1 1 \frac{1}{2} \frac{1}{2}\rangle, |1 1 \frac{3}{2} \frac{1}{2}\rangle \right)$

at $\delta = 0.5$, the matrix is

$\begin{pmatrix} 2.633 & 0 & 0 \\ 0 & 2.233 & 0 \\ 0 & 0 & 2.233\end{pmatrix}$

The eigen values and eigen state are

$\displaystyle 2.633 , |1 1 \frac{3}{2} \frac{3}{2}\rangle$

$\displaystyle 2.233 , |1 1 \frac{1}{2} \frac{1}{2}\rangle , |1 1 \frac{3}{2} \frac{1}{2}\rangle$

The lower level has mixed $j$ !

Of course, a more complete diagonalization should be admixture from other major shell as well. When we disable the LS coupling, using 7 major shells, the Nilsson diagram look like this

This is consistence with $omega_0 = const.$, i.e. no volume conservation. With LS coupling switched on.

Although it is not so clear, we can see the state avoid “crossing”.  The s-state are unaffected by the LS coupling.

One deflect in the plot is that, the slope seem to be incorrect. And also, I tried to add volume conservation, but the result does not consistence….. (sad).