Long time ago, I tried to tackle the Nilsson orbital by solving the Hamiltonian analytically. However, the Hamiltonian is without LS coupling. This times, I redo the calculation according to the reference B. E. Chi, Nuclear Phyiscs 83 (1966) 97-144.

The Hamiltonian is $\displaystyle H = \frac{P^2}{2m} + \frac{1}{2}m\left( \omega_\rho^2 (x^2+y^2) + \omega_z^2 z^2 \right) + C L\cdot S + D L\cdot L$

using $\displaystyle \omega_\rho = \omega_0 \left(\frac{3+2\delta}{3-4\delta}\right)^{1/6}$ $\displaystyle \omega_z = \omega_0 \left(\frac{3-4\delta}{3+2\delta}\right)^{1/3}$ $\displaystyle r^2 Y_{20}(\theta, \phi) = \frac{1}{4}\sqrt{\frac{5}{\pi}} (3z^2-r^2)$

The Hamiltonian becomes $\displaystyle H = -\frac{\hbar^2}{2m}\nabla^2 +\frac{1}{2} m \omega_0^2 r^2 - \frac{1}{2} m\omega_0^2 r^2 \beta Y_{20} + C L\cdot S + D L\cdot L$

Set $x_i^2 \rightarrow x_i^2 \frac{\hbar}{m \omega_0}$, and $r^2 \rightarrow \rho^2 \frac{\hbar}{ m \omega_0}$ $\displaystyle \frac{H}{m\omega_0} = \frac{1}{2}(-\nabla^2 + \rho^2) - \rho^2 \beta Y_{20} - 2 \kappa L\cdot S - \mu \kappa L\cdot L$

Set $\displaystyle \frac{H_0}{m\omega_0} = \frac{1}{2}(-\nabla^2 + \rho^2) - 2 \kappa L\cdot S - \mu \kappa L\cdot L$

and the perturbation is $\displaystyle \frac{H_p}{m\omega_0} = - \rho^2 Y_{20}$

The wavefunction for the spherical harmonic is $\displaystyle |Nljk\rangle = A r^l e^{-\frac{r^2}{2}} L_{\frac{N-l}{2}}^{l + \frac{1}{2}}(r^2) \sum_{m m_s} Y_{lm}(\theta, \phi) \chi_{\frac{1}{2} m_s} C_{lm\frac{1}{2} m_s}^{jk}$ $\displaystyle A = \sqrt{\frac{(\frac{N-l}{2})!(\frac{N+l}{2})! 2^{N+l+2}}{\sqrt{\pi} (N+l+1)!}}$

The diagonal elements are $\displaystyle \frac{1}{\hbar \omega_0 }\langle Nljk|H_0|Nljk\rangle = N + \frac{3}{2} - \kappa \langle L\cdot S \rangle - \mu \kappa l(l+1)$

where $\langle L \cdot S \rangle = \frac{1}{2} ( j(j+1) - l(l+1) - \frac{3}{4} )$

The off diagonal elements are $\displaystyle \frac{1}{\hbar \omega_0 }\langle Nljk|H_p|Nljk\rangle = - \langle Nljk| r^2 Y_{20}|Nljk\rangle$

( I will evaluate this integral in future )

The rest is diagonalization the Hamiltonian $\displaystyle H = H_0 + \beta H_p$

Here is the calculation for the 2nd harmonic for $\kappa = 0.05, \mu = 0$ The component of each orbital can be directly taken from the eigenvalue. Here is the 1/2 state. $\kappa = 0.05, \mu(N=3) = 0.35, \mu(N=4) = 0.625, \mu(N=5) = 0.63$ 