Nilsson Orbital using diagonalization method

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 [521]1/2 state. $\kappa = 0.05, \mu(N=3) = 0.35, \mu(N=4) = 0.625, \mu(N=5) = 0.63$

Rotational Band

For deformed nuclei, it can be rotated in various angular momentum in Laboratory frame. Assume rigid body rotation, the energy is

$\displaystyle E_J = \frac{1}{2}I\omega^2 = \frac{1}{2I}J^2$

In QM, that becomes

$\displaystyle H = \sum_{i=1}^{3} \frac{\hbar^2}{2I_i} J_i^2$

For axial symmetry, $I_1 = I_2 = I$

$\displaystyle H = \frac{\hbar^2}{2I} (J^2 - J_3^2) + \frac{\hbar^2}{2I_3}J_3^2$

Remember, in deformed nuclei, the projection of $J$ along the symmetry axis in the body-frame is $K$ . The expected value of the Hamiltonian with state $|Nn_z m_l K \rangle$ in the body-frame is proportional to $J(J+1)$ for $J^2$ and $K$ for $J_3$. i.e.

$\displaystyle E_J = \frac{\hbar^2}{2I} J(J+1) + E_K$

From body-frame to Lab-frame, we should apply the Wigner D-Matrix to the intrinsic wave function. ( I am not sure the following equation is correct, but the idea is rotating the body-frame wavefunction with Wigner D-Matrix to get the Lab-frame wave function. In Lab frame the total angular momentum must be a good Quantum number as rotational symmetry restored, so as $J_z = M$. The problem of the following equation is that the J is not a good Q-number in Nilsson wavefunction )

$\displaystyle |JMK\rangle = \sum_{M} D_{MK}^{J} |Nn_zm_lK\rangle$

However, the Wigner D-Matrix does not conserve parity transform:

$\displaystyle D_{MK}^J \rightarrow (-1)^{J+K} D_{M-K}^{J}$

In order to restored the parity, we need to include $\pm K$ in the Lab-frame wave function.

$\displaystyle |JMK\rangle = \sum_{M} \left( D_{MK}^J \pm (-1)^{J+K} D_{M-K}^J \right) |Nn_zm_lK\rangle$

where + for positive parity, – for negative parity.

From the above equation, for $K^\pi = 0^+$ ($0^-$), $J$ must be even (odd). For $K > 0$, $J = K, K+1, K+2, ...$.

We can see for $K = 1/2$, the $J = 5/2, 9/2, 11/2$ are lower to the main sequence. This was explained by adding an extra term in the rotation Hamiltonian that connect $\Delta K = 1$.

$\displaystyle \langle JMK | H'(\Delta K = 1) |JMK \rangle$

$\displaystyle \rightarrow \langle D_{MK}^J | H' | D_{MK}^J \rangle+ \langle D_{M-K}^J |H'| D_{MK}^J \rangle + \langle D_{MK}^J | H' | D_{M-K}^J \rangle+ \langle D_{M-K}^J |H'| D_{M-K}^J \rangle$

The term with $\Delta K = 0$ vanished. And since $\latex K = 1$, the only non-zero case is $K = 1/2$.

A possible form of the $H' (\Delta K = 1) = \frac{1}{2} \omega (J_+ + J_-)$. These are the ladder operator to rise or lower the m-component by 1. In 19F case, we can think it is a single proton on top of 18O core.  A rotation core affect the proton with an additional force, similar to Coriolis force on earth.

Axial Harmonic Oscillator – Nilsson Orbit (III)

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).

Axial Harmonic Oscillator – Nilsson Orbit (II)

This time we will project the deformed orbit into spherical orbit. The wavefunctions are stated in here again.

$\displaystyle \Psi^D_{N n_z m}(z, \rho, \phi) = |Nn_z m\rangle_D \\ = \sqrt{\frac{1}{\alpha_z \alpha^2}}\sqrt{\frac{ n_\rho !}{2^n_z n_z! (m + n_\rho)!\sqrt{\pi^3}}} H_{n_z}\left(\frac{z}{\alpha_z}\right) \\ \exp\left(- \frac{1}{2}\left(\frac{z^2}{\alpha_z^2}+\frac{\rho^2}{\alpha^2}\right)\right) \left(\frac{\rho}{\alpha}\right)^{m} L_{n_\rho}^{m} \left(\frac{\rho^2}{\alpha^2}\right) \exp(i m \phi)$

$\displaystyle \Psi^S_{Nlm}(r, \theta, \phi) = |N l m_l\rangle_S\\ =\sqrt{ \frac{1}{\sqrt{\pi}\alpha^{2l+3}} \frac{(\frac{N-l}{2})! (\frac{N+l}{2})! 2^{N+l+2}}{(N+l+1)!}} r^l \exp\left(-\frac{r^2}{2\alpha^2}\right) L_{k}^{l+\frac{1}{2}}\left( \frac{r^2}{\alpha^2} \right) Y_{lm}(\theta, \phi)$

Since both functions span the entire space and are basis, thus, we can related them as

$\displaystyle|Nn_z m \rangle_D = \sum_{N' l' m'} C_{N n_z m}^{N' l' m'} |N' l' m' \rangle_S$

where

$C_{N n_z m}^{N' l' m'} = \langle (N'l'm')_S | (N n_z m)_D \rangle$

First thing we notice is that $m' = m$ . Because the $phi$ components are the same in both wave function. i.e. $\int \exp(- i m' \phi) \exp(i m \phi) d\phi = 0$ if $m' \neq m$. We can omit the $m, m'$, so that $C_{N n_z}^{N' l'}$

Second thing is the parity must be the same, thus when $N$ is even (or odd), $N'$ must be even (or odd).

For non-deformed $\delta = 0$, $N' = N$,, here are some results

$\displaystyle|0 0 0 \rangle_D = |0 0 0\rangle_S$

$\displaystyle|1 1 0 \rangle_D = |1 1 0 \rangle_S$

$\displaystyle|1 0 1 \rangle_D = |1 1 1 \rangle_S$

$\displaystyle|2 2 0 \rangle_D = \sqrt{\frac{2}{3}} |2 2 0 \rangle_S - \sqrt{\frac{1}{3}} |2 0 0 \rangle_S$

$\displaystyle|2 1 1 \rangle_D = |2 2 1 \rangle_S$

$\displaystyle|2 0 2 \rangle_D = |2 2 2 \rangle_S$

$\displaystyle|2 0 0 \rangle_D = \sqrt{\frac{1}{3}} |2 2 0 \rangle_S + \sqrt{\frac{2}{3}} |2 0 0 \rangle_S$

$\displaystyle|3 3 0 \rangle_D = \sqrt{\frac{2}{5}} |3 3 0 \rangle_S - \sqrt{\frac{3}{5}} |3 1 0 \rangle_S$

$\displaystyle|3 2 1 \rangle_D = -\sqrt{\frac{4}{5}} |3 3 1 \rangle_S + \sqrt{\frac{1}{5}} |3 1 1 \rangle_S$

$\displaystyle|3 1 2 \rangle_D = |3 3 2 \rangle_S$

$\displaystyle|3 0 3 \rangle_D = |3 3 3 \rangle_S$

$\displaystyle|3 1 0 \rangle_D = \sqrt{\frac{3}{5}} |3 3 0 \rangle_S + \sqrt{\frac{2}{5}} |3 1 0 \rangle_S$

$\displaystyle|3 0 1 \rangle_D = -\sqrt{\frac{2}{5}} |3 3 1 \rangle_S - \sqrt{\frac{4}{5}} |3 1 1 \rangle_S$

For $\delta = 0.3$

$\displaystyle|0 0 0 \rangle_D = 0.995|0 0 0 \rangle_S + 0.099|2 2 0 \rangle_S - 0.015|2 0 0 \rangle_S + 0.009|4 4 0 \rangle_S + ...$

$\displaystyle|1 1 0 \rangle_D = 0.987|1 1 0 \rangle_S + 0.132|3 1 0 \rangle_S - 0.014|5 3 0 \rangle_S + 0.014|5 5 0 \rangle_S + ...$

$\displaystyle|1 0 1 \rangle_D = 0.994|1 1 1 \rangle_S + 0.108|3 3 1 \rangle_S + 0.017|3 1 1 \rangle_S + 0.011|5 5 1 \rangle_S + ...$

$\displaystyle|2 2 0 \rangle_D = 0.790|2 2 0 \rangle_S - 0.564|2 0 0 \rangle_S - 0.153|4 2 0 \rangle_S + 0.139|4 4 0 \rangle_S + ...$

$\displaystyle|2 1 1 \rangle_D = 0.986|2 2 1 \rangle_S + 0.158|4 4 1 \rangle_S - 0.052|4 2 1 \rangle_S - 0.012|6 4 1 \rangle_S + ...$

For $\delta = 0.6$

$\displaystyle|0 0 0 \rangle_D = 0.958|0 0 0 \rangle_S + 0.258|2 2 0 \rangle_S - 0.084|2 0 0 \rangle_S + 0.061|4 4 0 \rangle_S + ...$

$\displaystyle|1 1 0 \rangle_D = 0.885|1 1 0 \rangle_S + 0.319|3 3 0 \rangle_S - 0.274|3 1 0 \rangle_S - 0.119|5 3 0 \rangle_S + ...$

$\displaystyle|1 0 1 \rangle_D = 0.955|1 1 1 \rangle_S + 0.281|3 3 1 \rangle_S + 0.078|5 5 1 \rangle_S - 0.043|5 3 1 \rangle_S + ...$

$\displaystyle|2 2 0 \rangle_D = 0.598|2 2 0 \rangle_S - 0.45|2 0 0 \rangle_S - 0.379|4 2 0 \rangle_S + 0.299|4 4 0 \rangle_S \\ - 0.259|0 0 0 \rangle_S + 0.198|4 0 0 \rangle_S + 0.173|6 2 0 \rangle_S - 0.168|6 4 0 \rangle_S + ...$

$\displaystyle|2 1 1 \rangle_D = 0.882|2 2 1 \rangle_S + 0.381|4 4 1 \rangle_S - 0.191|4 2 1 \rangle_S - 0.129|6 6 1 \rangle_S \\ - 0.114|6 4 1 \rangle_S + ...$

We can see, more deform, more higher angular momentum states are involved.

Also, for a pure state when non-deform, the mixing is still small.

Axial Harmonic Oscillator – Nilsson Orbit

The Hamiltonian is

$\displaystyle H = -\frac{\hbar}{2m}\nabla^2 +\frac{m}{2}(\omega_\rho^2(x^2+y^2)+\omega_z z^2 )$

Use cylindrical coordinate, the Schrodinger equation is

$\displaystyle \left(-\frac{\hbar}{2m}\left(\frac{d^2}{dz^2} + \frac{1}{\rho}\frac{d}{d\rho}\left(\rho\frac{d}{d\rho}\right) + \frac{1}{\rho^2}\frac{d^2}{d\phi^2}\right) +\frac{m}{2}(\omega_\rho^2(\rho^2)+\omega_z z^2 ) \right) \Psi = E \Psi$

Where the energy is

$\displaystyle E = \hbar\omega_\rho(1+ n_x + n_y) + \hbar\omega_z\left(\frac{1}{2}+ n_z\right)$

Note that $n_x+n_y \neq n_\rho$. One of the reason is there are 2 degree of freedom, it cannot be solely expressed into 1 parameter.

Set $\Psi = Z P \Phi$,

$\displaystyle -\frac{\hbar}{2m}\left(\frac{d^2Z}{dz^2} P \Phi + \frac{1}{\rho}\frac{d}{d\rho}\left(\rho\frac{dP}{d\rho} \right) Z \Phi + \frac{1}{\rho^2}\frac{d^2\Phi}{d\phi^2} Z P\right) \\ + \frac{m}{2}\omega_\rho^2\rho^2 ZP\phi+\frac{m}{2}\omega_z z^2 ZP\Phi= E ZP\Phi$

$\displaystyle -\frac{\hbar}{2m}\left(\frac{d^2Z}{dz^2} \frac{1}{Z} + \frac{1}{\rho}\frac{d}{d\rho}\left(\rho\frac{dP}{d\rho} \right) \frac{1}{P} + \frac{1}{\rho^2}\frac{d^2\Phi}{d\phi^2} \frac{1}{\Phi}\right) \\ + \frac{m}{2}\omega_\rho^2\rho^2 +\frac{m}{2}\omega_z z^2 = E$

The angular part, we can set

$\displaystyle \frac{d^2\Phi}{d\phi^2} \frac{1}{\Phi} = -m_\phi^2$

The solution is $\Phi(\phi) = \exp(i m_\phi \phi)$

The z-part is usual 1D Harmonic oscillator

Thus the rest is

$\displaystyle -\frac{\hbar}{2m}\left( \frac{1}{\rho}\frac{d}{d\rho}\left(\rho\frac{dP}{d\rho} \right) \frac{1}{P} - \frac{m_\phi^2}{\rho^2}\right) + \frac{m}{2}\omega_\rho^2\rho^2 = \hbar\omega_\rho(1+ n_x + n_y)$

rearrange

$\displaystyle -\frac{\hbar}{2m}\left( \frac{1}{\rho}\frac{d}{d\rho}\left(\rho\frac{dP}{d\rho} \right) - \frac{m_\phi^2}{\rho^2}P\right) + \frac{m}{2}\omega_\rho^2\rho^2 P - \hbar\omega_\rho(1+ n_x + n_y) P = 0$

Set a dimensionless constant $\alpha^2 = \frac{\hbar}{m \omega_\rho}$, and $x = \rho/\alpha$

$\displaystyle \frac{1}{x}\frac{d}{dx}\left(x\frac{dP}{dx}\right) + \left( 2(1+n_x+n_y) - \frac{m_\phi^2}{x^2} - x^2 \right) P = 0$

Using the normalisation formula, $\int P^2 x dx = 1$, set $u = P \sqrt{x}$

$\displaystyle x^2\frac{d^2u}{dx^2}+ \left( 2(1+n_x+n_y) - \frac{m_\phi^2-\frac{1}{4}}{x^2} - x^2 \right) u = 0$

Because of the long-range and short-range behavior, set

$\displaystyle u(x) = f(x) x^{m_\phi+\frac{1}{2}} \exp\left(-\frac{x^2}{2}\right)$

$\displaystyle x\frac{d^2f}{dx^2} + (1+2m_\phi-2x^2)\frac{df}{dx} + 2(n_x+n_y-m_\phi)x f= 0$

Set $y = x^2$

$\displaystyle y \frac{d^2f}{dy^2} + (m_\phi+1-y)\frac{df}{dy} + \frac{n_x+n_y-m_\phi}{2} f= 0$

Define $n_\rho = \frac{n_x+n_y-m_\phi}{2} \rightarrow n_x+n_y = 2 n_\rho + m_\phi = N - n_z$.

$\displaystyle y \frac{d^2f}{dy^2} + (m_\phi+1-y)\frac{df}{dy} + n_\rho f= 0$

This is our friend again! The complete solution is

$\displaystyle \Psi_{n_z n_\rho m_\phi}(z, \rho, \phi) \\ = \sqrt{\frac{1}{\alpha_z \alpha^2}}\sqrt{\frac{ n_\rho !}{2^n_z n_z! (m_\phi + n_\rho)!\sqrt{\pi^3}}} H_{n_z}\left(\frac{z}{\alpha_z}\right) \\ \exp\left(- \frac{1}{2}\left(\frac{z^2}{\alpha_z^2}+\frac{\rho^2}{\alpha^2}\right)\right) \left(\frac{\rho}{\alpha}\right)^{m_\phi} L_{n_\rho}^{m_\phi} \left(\frac{\rho^2}{\alpha^2}\right) \exp(i m_\phi \phi)$

where $\alpha_z^2 = \hbar/m/\omega_z$

The energy is

$\displaystyle E = \hbar\omega_z \left(\frac{1}{2} + n_z \right) + \hbar\omega_\rho \left( 2n_\rho + m_\phi + 1 \right)$

The quantum number $m_\phi$ has same meaning as $m_l$ in spherical case.

The notation for the state is

$|Nn_z m_l \rangle$

with the spin, the only good quantum number is the z-component of the total angular momentum $J = L+S$  long the body axis, i.e. $K = m_\phi \pm 1/2$, thus, the state is

$|Nn_z m_l K \rangle$

Note that the total angular momentum $J^2$ is not a good quantum number in deformation, as the rotational symmetry is lost. However, the quantum number $K$ is linked with the angular momentum of the Nilsson single particle orbit. It is because when a particle has $m_j = K$, the angular momentum must at least $j \geq K$.

The above is a general solution for the harmonic oscillator in cylindrical coordinate. When $\omega_z = \omega_\rho$, it reduce to spherical case.

According to P. Ring & P. Schuck (2004) (The Nuclear Many-Body Problem, P.68), the Hamiltonian can be expressed as a quadruple deform field by setting

$\displaystyle \omega_\rho^2 = \omega^2 \left(1+\frac{2}{3} \delta\right) \\ \omega_z^2 = \omega^2 \left(1-\frac{4}{3} \delta\right)$

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

This Hamiltonian has similarity with the deformation

$\displaystyle R(\theta, \phi) = R_0 \left(1 + \alpha_00 + \sum_{\lambda=1}^{\infty}\sum_{\mu=-\lambda}^{\lambda} \alpha_{\lambda \mu} Y_{\lambda \mu}(\theta, \phi) \right)$

take the first quadruple term, and calculate $R^2(\theta, \phi)$

$R^2(\theta, \phi) = R_0^2 (1 + 2 \beta Y_{20}(\theta,\phi) )$

Compare, we have

$\displaystyle \beta = \frac{1}{3}\sqrt{\frac{16\pi}{5}}\delta \approx 1.05689 \delta$

The ratio

$\displaystyle \frac{\omega_z}{\omega_\rho} = \left(\frac{\alpha_\rho}{\alpha_z} \right) = \sqrt{\frac{1-\frac{4}{3}\delta}{1+\frac{2}{3}\delta}} = 1- \delta + \frac{1}{6}\delta^2 - \frac{5}{18}\delta^3...$

We need volume conservation

$\omega_x \omega_y \omega_z = \omega_0^3$

Thus,

$\displaystyle \omega = \omega_0 \left(\frac{1}{(1+2/3\delta)(1-4/3\delta)} \right)^{\frac{1}{6}} \approx \omega_0 \left(1+\frac{2}{9}\delta^2 + \frac{8}{81}\delta^3 \right)$

The energy is

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

Here is some plots with various $\delta = 0, 0.3, 0.5$

Next time, I will add the LS coupling and L-square term to recreate the Nilsson diagram. Also I will expand the solution from cylindrical coordinate into spherical coordinate. This unitary transform is the key to understand the single particle-ness.