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.