Quantum harmonic oscillator

We kicked off with a new tutorial series on quantum optics. Our team regularly joins internal training sessions on the physics of quantum computers. With the core purpose to excel in quantum algorithm design we parallelly grow awareness of the fundamental limits imposed by physics. Understanding hardware operation details and acknowledging the latest technological breakthroughs gives us solid position to push the frontiers in quantum computing theory.

We share sample notes from an introductory session covering the fundamentals of the quantum harmonic oscillator model – a foundation of many quantum computer architectures.

1-D quantum harmonic oscillator

Below we discuss the quantum harmonic oscillator model. The Hamiltonian for the 1-D harmonic oscillator can be written as

\begin{equation} \hat{H}_{1D}=\frac{\hat{p}_{x}^{2}}{2\mu}+\frac{1}{2}\mu\omega^{2}\hat{x}^{2} \label{Harm:Ham1} \end{equation}

where \(\hat{p}_{x}\) denotes the \(x\)-component of the momentum operator and \(\hat{x}\) stands for the \(x\)-component of the position operator. \(\mu\) is the reduced mass of the system. We want to solve the stationary Schrödinger equation:

\begin{equation} \hat{H}_{1D}|\psi\rangle=E|\psi\rangle \label{harm:SE} \end{equation}

For wavefunctions belonging to \(L^{2}(\mathbb{R},dx)\), that is to the set of square-integrable functions, the Hamiltonian given in eq. \ref{Harm:Ham1} is Hermitian and has a discrete and positive spectrum. Eq. \ref{harm:SE} can be solved analytically using e.g. the Frobenius method. Here, we follow a more elegant algebraic approach to solving eq. \ref{harm:SE}.

First introduce a set of dimensionless operators by scaling the position and momentum operators. By doing so, one obtains convenient units of length, momentum and energy, which are natural for the harmonic oscillator. The new operators are defined as

\begin{equation} \label{Harm:alpha} \hat{x}=\alpha\hat{X}, \; \alpha \in R \qquad \hat{p}_{x}=\frac{\hslash}{\alpha}\hat{P} \end{equation}

where \(\alpha\) has the dimension of length. The Hamiltonian in the new dimensionless representation reads

\begin{equation} \hat{H}_{1D}=\frac{\hslash^{2}}{2\mu\alpha^{2}}\hat{P}_{x}^{2}+\frac{1}{2}\mu\omega^{2}\alpha^{2}\hat{X}^{2} \label{harm:ham2} \end{equation}

At this point one has the freedom to choose the value for the scaling constant \(\alpha\). Requiring the constants standing by \(\hat{P}_{x}^{2}\) and \(\hat{X}^{2}\) in eq. \ref{harm:ham2} to be equal, sets \(\alpha=\sqrt{\frac{\hslash}{\mu\omega}}\). Then the dimensionless position and momentum operators take the form

\begin{equation} \hat{X}=\sqrt{\frac{\mu\omega}{\hslash}}\hat{x}, \qquad \hat{P}=\sqrt{\frac{1}{\hslash\omega\mu}}\hat{p}_{x} \end{equation}

and the Hamiltonian simplifies to

\begin{equation} \hat{H}_{1D}=\frac{\hslash\omega}{2}\left(\hat{P}^{2}+\hat{X}^{2}\right) \label{harm:ham3} \end{equation}

The Hamiltonian given in eq. \ref{harm:ham3} can be cast in a more compact form with the following linear combination of the \(\hat{X}\) and \(\hat{P}\) operators :

\begin{equation} \label{creationdefinition} a=\frac{1}{\sqrt{2}}\left(\hat{X}+i\hat{P}\right), \qquad a^{\dagger}=\frac{1}{\sqrt{2}}\left(\hat{X}-i\hat{P}\right) \end{equation}

The new operators (non-unitary and non-Hermitian) are called the annihilation and creation operators, respectively. The creation and annihilation operators are each other's Hermitian conjugates and together they constitute ladder operators for the 1-D harmonic oscillator. The Hamiltonian operator expressed in these new operators is given by

\begin{equation}\label{harm:ladderHamiltonian} \hat{H}_{1D}=\hslash\omega\left(a^{\dagger}a+\frac{1}{2}\right) \end{equation}

The form of the above representation to the 1-D harmonic oscillator Hamiltonian indicates that it is a sum of a single operator \(\hat{N}:=a^{\dagger}a\) and a constant term \(\frac{1}{2}\hslash\omega\). The former is called the number operator and the latter is the zero-point energy of the 1-D harmonic oscillator. It is evident that if function \(f\) is an eigenfunction of \(\hat{N}\), so it is an eigenfunction of \(\hat{H}_{1D}\). Assume that \(\hat{N}f=\lambda f\). Then

\begin{equation} \lambda=\langle f|\lambda f\rangle=\langle f|\hat{N}f\rangle=\langle\hat{N}\rangle_{f}=\langle af|af\rangle=\|af\|^{2}\ge 0 \Rightarrow \lambda \ge 0 \end{equation}

We can see that \(\hat{N}\) as well as \(\hat{H}_{1D}\) have a non-negative spectrum. Does there exist an eigenfunction which corresponds to the lowest possible \(\lambda=0\)? Assuming that \(\lambda=0\) exists we can write

\begin{equation} \hat{N}f_{0}=0 \end{equation}

for some \(f_0\). We then find that \(\langle f_{0}|\hat{N}f_{0}\rangle=0 \Rightarrow \langle f_{0}|\hat{N}f_{0}\rangle=\|af_{0}\|^{2}=0 \Rightarrow \|af_{0}\|=0\). As a consequence, in order to obtain the generic solution to the Schrödinger equation for \(\lambda=0\) we need to solve the equation \(af_{0}=0\), which can be readily written as

\begin{equation} \left(X+\frac{d}{dX}\right)f_{0}\left(X\right)=0 \end{equation}

with standard normalised textbook solution:

\begin{equation} f_{0}\left(X\right)=\pi^{-\frac{1}{4}}e^{-\frac{1}{2}X^{2}} \end{equation}

Our goal now is to derive a general formula for the \(n\)-th eigenfunction of \(\hat{N}\), that is for \(f_{n}\), utilizing only the creation and annihilation operators. This can be done with the aid of the following commutation relations:

\begin{equation} \left[a,a^{\dagger}\right]=\frac{1}{2}\left[\hat{X}+i\hat{P},\hat{X}-i\hat{P}\right]=\mathrm{1} \end{equation}

where we made use of the canonical commutation relations in the dimensionless representation \(\left[\hat{X},\hat{P}\right]=i\). The commutation relation between the creation and annihilation operator is helpful in calculating the following commutator

\begin{equation} \label{harm:laddercommutation} \left[a,\left(a^{\dagger}\right)^{n}\right]=a^{\dagger}\left[a,\left(a^{\dagger}\right)^{n-1}\right]+\left[a,a^{\dagger}\right]\left(a^{\dagger}\right)^{n-1}=...= n\left(a^{\dagger}\right)^{n-1} \end{equation}

We see that for \(n=2\): \(\left[a,\left(a^{\dagger}\right)^{n}\right]=2a^{\dagger}\). By induction, assuming that for some arbitrary \(n\) eq.\ref{harm:laddercommutation} holds, we find that:

\begin{equation} \left[a,\left(a^{\dagger}\right)^{n+1}\right]=a^{\dagger}\left[a,\left(a^{\dagger}\right)^{n}\right]+\left[a,a^{\dagger}\right]\left(a^{\dagger}\right)^{n}=n\left(a^{\dagger}\right)^{n-1}+ \left(a^{\dagger}\right)^{n}=\left(n+1\right)\left(a^{\dagger}\right)^{n} \end{equation}

But at the same time we have

\begin{equation} \left[\hat{N},\left(a^{\dagger}\right)^{n}\right]=n\left(a^{\dagger}\right)^{n} \end{equation}

which gives

\begin{equation} \left[\hat{N},\left(a^{\dagger}\right)^{n}\right]f_{0}=n\left(a^{\dagger}\right)^{n}f_{0} \end{equation}

To summarize, the eigenvalue problem for \(\hat{N}\) takes the form

\begin{equation} \hat{N}\left(a^{\dagger}\right)^{n}f_{0}=n\left(a^{\dagger}\right)^{n}f_{0} \end{equation}

which means that \(\left(a^{\dagger}\right)^{n}f_{0}\) is the eigenfunction of \(\hat{N}\) with eigenvalue \(n\), hence in fact must be proportional to \(f_{n}\). For this reason \(\hat{N}\) is called number operator. Normalization of \(f_{n}\) is straightforward:

\begin{equation} \begin{split} \|\left(a^{\dagger}\right)^{n}f_{0}\|^{2}=\langle \left(a^{\dagger}\right)^{n}f_{0}|\left(a^{\dagger}\right)^{n}f_{0}\rangle = \\ =\langle \left(a^{\dagger}\right)^{n-1}f_{0}| a\left(a^{\dagger}\right)^{n}f_{0}\rangle=n\langle \left(a^{\dagger}\right)^{n-1}f_{0}|\left(a^{\dagger}\right)^{n-1}f_{0}\rangle =...= n! \end{split} \end{equation}

where we used eq.~\ref{harm:laddercommutation} and the fact that \(af_{0}=0\). Finally, the normalized eigenfunctions are

\begin{equation}\label{harm:fn} f_{n}=\frac{1}{\sqrt{n!}}\left(a^{\dagger}\right)^{n}f_{0} \end{equation}

and the Schrödinger equation can be written as:

\begin{equation} \hat{H}_{1D}f_{n}=\hslash\omega\left(\hat{N}+\frac{1}{2}\right)f_{0}=\hslash\omega\left(n+\frac{1}{2}\right)f_{n}, \quad n=0,1,2,... \end{equation}

Functions \ref{harm:fn} form an orthonormal set in \(L^{2}\left(\mathbb{R},dx\right)\) as they are eigenfunctions of a Hermitian operator. For a more explicit form of the 1-D harmonic oscillator wavefunctions let us inspect eq. \ref{harm:fn}:

\begin{equation}\label{harm:explicit} f_{n}=\frac{1}{\sqrt{n!}}\left(a^{\dagger}\right)^{n}f_{0}=\frac{1}{\sqrt{2^{n}\pi^{\frac{1}{2}}n!}}\left(X-\frac{d}{dX}\right)^{n}e^{-\frac{1}{2}X^{2}} \end{equation}

According to the identity

\begin{equation}\label{harm:explicit2} -\frac{d}{dX}\left(g(X)e^{-\frac{1}{2}X^{2}}\right)=e^{-\frac{1}{2}X^{2}}\left(X-\frac{d}{dX}\right)g(X) \end{equation}

for any differentiable function \(g(X)\) we can write

\begin{equation}\label{harm:explicit3} (-1)^{n}\frac{d^{n}}{dX^{n}}\left(g(X)e^{-\frac{1}{2}X^{2}}\right)=e^{-\frac{1}{2}X^{2}}\left(X-\frac{d}{dX}\right)^{n}g(X) \end{equation}

By pasting \(g(X)=e^{-\frac{1}{2}X^{2}}\) we finally get

\begin{equation}\label{harm:explicit4} (-1)^{n}e^{\frac{1}{2}X^{2}}\frac{d^{n}}{dX^{n}}\left(e^{-X^{2}}\right)=\left(X-\frac{d}{dX}\right)^{n}e^{-\frac{1}{2}X^{2}} \end{equation}

Hence,

\begin{equation}\label{harm:harm:1Dsolution} f_{n}\left(X\right)=\frac{1}{\sqrt{2^{n}\pi^{\frac{1}{2}}n!}}(-1)^{n}e^{\frac{1}{2}X^{2}}\frac{d^{n}}{dX^{n}}e^{-X^{2}}:= \frac{1}{\sqrt{2^{n}\pi^{\frac{1}{2}}n!}}H_{n}\left(X\right)e^{-\frac{1}{2}X^{2}} \end{equation}

where \(H_{n}\left(X\right)\) is \(n-th\) Hermite polynomial defined as \(H_{n}\left(X\right):=(-1)^{n}e^{\frac{1}{2}X^{2}}\frac{d^{n}}{dX^{n}}e^{-\frac{1}{2}X^{2}}\). In the original representation of the position operator the 1-D harmonic oscillator wavefunction is given by

\begin{equation}\label{harm:explicit5} \psi_{n}\left(x\right)=\frac{1}{\sqrt{2^{n}\pi^{\frac{1}{2}}n!}}\left(\frac{\mu\omega}{\hslash}\right)^{\frac{1}{4}}H_{n}\left(\sqrt{\frac{\mu\omega}{\hslash}}x\right)e^{-\frac{1}{2}\frac{\mu\omega}{\hslash}x^{2}} \end{equation}

Lastly, let us derive the matrix elements of the creation and annihilation operators in the harmonic oscillator eigenbasis. The action of the annihilation operator \(a\) on the \(n-th\) eigenstate is given by

\begin{equation}\label{harm:explicit6} \begin{split} af_{n}=\frac{1}{\sqrt{n!}}a\left(a^{\dagger}\right)^{n}f_{0}=\frac{1}{\sqrt{n!}}\left[a,\left(a^{\dagger}\right)^{n}\right]f_{0}=\\ =\frac{n}{\sqrt{n!}}\left(a^{\dagger}\right)^{n-1}f_{0}=\sqrt{n} \frac{1}{\sqrt{(n-1)!}}\left(a^{\dagger}\right)^{n-1}f_{0}=\sqrt{n}f_{n-1} \end{split} \end{equation}

For the creation operator we can write

\begin{equation}\label{harm:explicit7} a^{\dagger}f_{n}=\frac{1}{\sqrt{n!}}\left(a^{\dagger}\right)^{n+1}f_{0}=\frac{\sqrt{n+1}}{\sqrt{(n+1)!}}\left(a^{\dagger}\right)^{n+1}f_{0}=\sqrt{n+1}f_{n+1} \end{equation}

The above relations can be used to calculate matrix elements of the position and momentum operator in the eigenbasis \(\{f_{n}\}_{n=0,1,2,...}\). In this basis the matrices for the creation and annihilation operators are given respectively as

\begin{equation} \left[a^{\dagger}\right]= \left(\begin{array}{cccccc} 0 & 0 & 0 & 0 & 0 &\ldots\\ \sqrt{1} & 0 & 0 & 0 & 0 & \ldots\\ 0 & \sqrt{2} & 0 & 0 & 0 &\ldots\\ 0 & 0 & \sqrt{3} & 0 & 0 & \ldots\\ \vdots \end{array}\right) \end{equation}
\begin{equation} \left[a\right]= \left(\begin{array}{cccccc} 0 &\sqrt{1}& 0 & 0 & 0 &\ldots\\ 0 & 0&\sqrt{2} & 0 & 0 & \ldots\\ 0 & 0 & 0 & \sqrt{3} & 0 &\ldots\\ 0 & 0 & 0 & 0 & \sqrt{4} & \ldots\\ \vdots \end{array}\right) \end{equation}

3-D isotropic harmonic oscillator in spherical coordinates

In the Cartesian coordinates the Hamiltonian for the isotropic three-dimensional harmonic oscillator has the form

\begin{equation}\label{harm:3Dcart} \hat{H}_{3D}(x,y,z)=-\frac{\hslash^{2}}{2\mu}\Delta+\frac{1}{2}\mu\omega^{2}\left(\hat{x}^{2}+\hat{y}^{2}+\hat{z}^{2}\right) \end{equation}

where \(\mu\) is the reduced mass of the system, \(\omega\) is the harmonic frequency and \(\Delta\) is the Laplacian. The form of the above Hamiltonian suggests that it is spherically symmetric. Therefore a reasonable choice for a coordinate system is spherical coordinates

\begin{equation} \hat{H}_{3D}(r,\theta,\phi)=-\frac{\hslash^{2}}{2\mu r^{2}}\frac{\partial}{\partial r}\left(r^{2}\frac{\partial}{\partial r}\right)+\frac{1}{2\mu r^{2}}\hat{L}^{2}+\frac{1}{2}\mu\omega^{2}r^{2} \end{equation}

where, \(\hat{L}\) stands for total angular momentum operator. Because \(\hat{L}\) depends only on spherical angles, it commutes with the quadratic potential and the radial kinetic energy operator \(-\frac{\hslash^{2}}{2\mu r^{2}}\frac{\partial}{\partial r}\left(r^{2}\frac{\partial}{\partial r}\right)\). This result can be proved using the following commutation relations:

\begin{equation*} \left[\hat{L}_{i},\hat{p}_{j}^{2}\right]=\epsilon_{ikl}\left[\hat{x}_{k}\hat{p}_{l},\hat{p}_{j}^{2}\right]= \epsilon_{ikl}\left[\hat{x}_{k}, \hat{p}_{j}^{2}\right]\hat{p}_{l}=2i\hslash\epsilon_{ikl}\delta_{kj}\hat{p}_{j}\hat{p}_{l}=0 \end{equation*}

since we sum in default over \(k,l\) and the expression is a product of totally symmetric \(\hat{p}_{j}\hat{p}_{l}\) and the totally antisymmetric species \(\epsilon_{ikl}\). For this reason there must exist a common eigenbasis for the total Hamiltonian and the angular part of the Hamiltonian. This leads to factorization of the total wavefunction for 3-D harmonic oscillator

\begin{equation} \psi(r,\theta,\phi)=R_{nl}(r)Y_{lm}(\theta,\phi) \end{equation}

After some manipulations the Schrödinger equation reads

\begin{equation} \left(-\frac{\hslash^{2}}{2\mu r^{2}R(r)_{nl}}\left(\frac{\partial}{\partial r}\left(r^{2}\frac{\partial}{\partial r}\right)+\frac{1}{2}\mu\omega^{2}r^{2}\right)\right)R_{nl}(r)+\frac{1}{2\mu r^{2}Y_{lm}(\theta,\phi)}\hat{L}^{2}Y_{lm}(\theta,\phi)=E \end{equation}

The angular part of the above equation is the 3-D homogeneous Laplace equation, with spherical harmonics as solutions:

\begin{equation} \frac{1}{2\mu r^{2}}\hat{L}^{2}Y_{lm}(\theta,\phi)=\frac{\hslash^{2}l(l+1)}{2\mu r^{2}}{Y_{lm}(\theta,\phi)}\equiv E_{\text{rot}}Y_{lm}(\theta,\phi) \end{equation}

and because the angular part depends also on \(r\) it must be incorporated in the radial equation, providing a coupling between the rotational and vibrational motion.

\begin{equation}\label{harm:radial} \left(-\frac{\hslash^{2}}{2\mu r^{2}}\frac{\partial}{\partial r}\left(r^{2}\frac{\partial}{\partial r}\right)+\frac{1}{2}\mu\omega^{2}r^{2}\right)R_{nl}(r)+\frac{\hslash^{2}l(l+1)}{2\mu r^{2}}R_{nl}(r)=E R_{nl}(r) \end{equation}

Analogically to the one-dimensional case we introduce dimensionless quantities: \(\alpha=\sqrt{\frac{\hslash}{\mu\omega}}\), which yields a representation of the dimensionless position operator:

\begin{equation} \hat{\rho}=\alpha^{-1}\hat{r} \end{equation}

so that eq. \ref{harm:radial} takes the form

\begin{equation}\label{harm:radial1} \frac{d^{2}R(\rho)}{d\rho^{2}}+\frac{2}{\rho}\frac{dR(\rho)}{d\rho}-\rho^{2}R(\rho)-\frac{l(l+1)}{\rho^{2}}R(\rho)+\frac{2E}{\hslash\omega}R(\rho)=0 \end{equation}

noting that \(R_{nl}(r)\) and \(R(\rho)\) are different functions. In order to predetermine an appropriate ansatz for the radial wavefunction it is a common practice to investigate the asymptotic behaviour of the solutions. Let us focus on eq. \ref{harm:radial1}. In the limit of large \(\rho\) the equation reduces to

\begin{equation} \frac{d^{2}R(\rho)}{d\rho^{2}}-\rho^{2}R(\rho)=0 \end{equation}

providing Gaussian-type solutions:

\begin{equation} R(\rho)=e^{-\frac{1}{2}\rho^{2}} \end{equation}

In the opposite limit of small \(\rho\) the equation

\begin{equation} \frac{d^{2}R(\rho)}{d\rho^{2}}+\frac{2}{\rho}\frac{dR(\rho)}{d\rho}-\frac{l(l+1)}{\rho^{2}}R(\rho)=0 \end{equation}

suggests the following ansatz \(R(\rho)=\rho^{t}\), which generates the algebraic condition

\begin{equation} t(t-1)+2t-l(l+1)=0 \end{equation}

setting value to the exponent at \(t=l\). Finally it is possible to propose a reasonable ansatz for the radial wavefunction in the following form

\begin{equation} R(\rho)=\rho^{l}\sum_{k=0}^{+\infty}a_{k}\rho^{k}e^{-\frac{1}{2}\rho^{2}} \end{equation}

First we compute derivatives

\begin{equation} \frac{dR(\rho)}{d\rho}=\sum_{k=0}^{+\infty}a_{k}\left[(l+k)\rho^{k+l-1}-\rho^{l+k+1}\right]e^{-\frac{1}{2}\rho^{2}} \end{equation}
\begin{equation} \frac{d^{2}R(\rho)}{d\rho^{2}}=\sum_{k=0}^{+\infty}a_{k}\left[(l+k)(l+k-1)\rho^{k+l-2}-(2l+2k+1)\rho^{l+k}+\rho^{l+k+2}\right]e^{-\frac{1}{2}\rho^{2}} \end{equation}

Next we can place the above expressions into \ref{harm:radial1}, divide by the exponent factor and obtain

\begin{equation} \begin{split} \sum_{k=0}^{+\infty}a_{k}\biggr[(l+k)(l+k-1)\rho^{k+l-2}-(2l+2k+1)\rho^{l+k}+\rho^{l+k+2}+2(l+k)\rho^{k+l-2}\\ -2\rho^{l+k}-\rho^{l+k+2}-l(l+1)\rho^{l+k-2}+\frac{2E}{\hslash\omega}\rho^{l+k}\biggr]=0 \end{split} \label{eq:3dhofrob} \end{equation}

Note that powers in above equation differ by two. After some rearrangement we get

\begin{equation} \sum_{k=0}^{+\infty}a_{k}\left[k(2l+k+1)\rho^{k+l-2}+\left(\frac{2E}{\hslash\omega}-(2l+2k+3)\right)\rho^{l+k}\right]=0 \end{equation}

First two terms for \(k=0,1\) in the above equation determine coefficients \(a_{0},a_{1}\). For \(k=0\) the expression standing by \(a_{0}\) vanishes so \(a_{0}\) is allowed in the solution and is undetermined (may be chosen arbitrarily). For \(k=1\):

\begin{equation} a_{1}\left[(2l+2)\rho^{l-1}+\left(\frac{2E}{\hslash\omega}-(2l+5)\right)\rho^{l+1}\right]=0 \end{equation}

we can see that the function in brackets do not vanish for any \(l\) value, hence \(a_{1}\) must be zero. By shifting the summation in eq. \ref{eq:3dhofrob} to obtain equal powers in the monomials we get

\begin{equation} \sum_{k=2}^{+\infty}\left[a_{k+2}(k+2)(2l+k+3)+a_{k}\left(\frac{2E}{\hslash\omega}-(2l+2k+3)\right)\right]\rho^{l+k}=0 \end{equation}

Forcing all coefficients standing by the powers of \(\rho\) to vanish results in the following recurrence relation:

\begin{equation} a_{k+2}=-\frac{\frac{2E}{\hslash\omega}-(2l+2k+3)}{(k+2)(2l+k+3)}a_{k} \label{eq:recurrence} \end{equation}

In the limit of large \(k\) this series behaves as a harmonic-like, i.e. \(a_{k+2}=\frac{2}{k}a_{k}\), hence is divergent. Therefore there must exist such \(s>0\) for which \(a_{s}=0\), which implies all subsequent coefficients to be also zero. From the above recurrence relation is it possible to read a condition on allowed energy levels for the 3-D harmonic oscillator

\begin{equation} \frac{2E}{\hslash\omega}-(2l+2s+3)=0 \end{equation}

and after some rearrangement we arrive at the familiar expression for eigenvalues of the 3-D harmonic oscillator:

\begin{equation} \label{eigenenergy} E=\hslash\omega \left( l+s+\frac{3}{2} \right) \end{equation}

Note that here the total energy depends explicitly on the \(l\)-angular momentum quantum number, meaning that there exists a rotational-vibrational coupling, which manifests itself in the dependence of the height of the potential barrier on the rotational state. Inserting the expression for the total energy into the recurrence relation given in eq. \ref{eq:recurrence} gives

\begin{equation} a_{k+1}=\frac{k-s}{(k+1)(l+k+\frac{3}{2})}a_{k} \end{equation}

and the final form of the solution is given by

\begin{equation} R_{nl}(\rho)=\sum_{k=0}^{n}a_{k}(s)\rho^{l+2k}e^{-\frac{1}{2}\rho^{2}} \end{equation}

The total energy can be written as

\begin{equation} E=\hslash\omega\left(l+2s+\frac{3}{2}\right)=\hslash\omega\left(n+\frac{3}{2}\right),\qquad n=0,1,2,... \end{equation}

The total wavefunction of the 3-D isotropic harmonic oscillator is given by

\begin{equation}\label{harm:ultimatesolution} \psi(r,\theta,\phi)=R_{nl}(r)Y_{lm}(\theta,\phi) \end{equation}

The radial part of this wavefunction can be expressed by associated Laguerre polynomials:

\begin{equation}\label{sphericaloscillator} R_{nl}(r)=N_{nl}r^{l}L^{l+\frac{1}{2}}_{n}\left(\frac{\mu\omega}{\hslash}r^2\right)e^{-\frac{1}{2}\frac{\mu\omega}{\hslash}r^{2}} \end{equation}

where the normalization constant is given by

\begin{equation*} N_{nl}=\left(\frac{\mu\omega}{\hslash}\right)^{\frac{1}{2}l+\frac{3}{4}}\sqrt{\frac{2^{\left(n+2l+\frac{3}{2}\right)}n!}{\pi(2n+2l+1)!!}} \end{equation*}

Functions from eq. \ref{sphericaloscillator} are called the spherical oscillator functions and are commonly used as basis functions in solving nuclear motion problems. The energy levels of the 3-D isotropic harmonic oscillator are \(\frac{1}{2}(n+1)(n+2)\)-degenerate. This result can be derived in a straightforward way when one considers a 3-D harmonic oscillator in the Cartesian representation and counts how many possible combinations of quantum numbers \((n_x,n_y,n_z)\) associated with the \(x\), \(y\) and \(z\) dimension, respectively, add up to the same value.