# Isotropic harmonic oscillator in 3-d – use of spherical harmonics

Shankar, R. (1994), Principles of Quantum Mechanics, Plenum Press. Chapter 12, Exercise 12.6.11.

[If some equations are too small to read easily, use your browser’s magnifying option (Ctrl + on Chrome, probably something similar on other browsers).]

We’ve solved the 3-d isotropic harmonic oscillator before, so we’ve already solved most of Shankar’s exercise 12.6.11. We can quote the results here. The solution has the form

$\displaystyle \psi_{Elm}=\frac{U_{El}\left(r\right)}{r}Y_{l}^{m}\left(\theta,\phi\right) \ \ \ \ \ (1)$

The earlier solution uses notation from Griffiths’s book, but as the end result is the same, it’s not worth going through the derivation again using Shankar’s notation.

The potential is

$\displaystyle V(r)=\frac{1}{2}m\omega^{2}r^{2} \ \ \ \ \ (2)$

The radial equation to be solved is

$\displaystyle \frac{d^{2}u}{d\rho^{2}}=\left(-1+\frac{l(l+1)}{\rho^{2}}+\rho_{0}^{2}\rho^{2}\right)u \ \ \ \ \ (3)$

If we define

 $\displaystyle \kappa^{2}$ $\displaystyle \equiv$ $\displaystyle \frac{2\mu E}{\hbar^{2}}\ \ \ \ \ (4)$ $\displaystyle \rho$ $\displaystyle \equiv$ $\displaystyle \kappa r\ \ \ \ \ (5)$ $\displaystyle \rho_{0}$ $\displaystyle \equiv$ $\displaystyle \frac{\mu\omega}{\hbar\kappa^{2}}=\frac{\hbar\omega}{2E} \ \ \ \ \ (6)$

Taking the asymptotic behaviour of the radial function for small and large ${r}$ into account leads us to a solution of form

$\displaystyle u(\rho)=\rho^{l+1}e^{-\rho_{0}\rho^{2}/2}v(\rho) \ \ \ \ \ (7)$

Note that Griffiths’s ${v}$ is not the same as Shankar’s ${v}$, the latter of which is defined by Shankar’s equation 12.6.49.

This gives a differential equation for Griffiths’s ${v}$

$\displaystyle \rho\frac{d^{2}v}{d\rho^{2}}+2\left(l+1-\rho_{0}\rho^{2}\right)\frac{dv}{d\rho}+\rho(1-\rho_{0}(2l+3))v=0 \ \ \ \ \ (8)$

The function ${v}$ can be solved as a power series, giving

$\displaystyle v(\rho)=\sum c_{j}\rho^{j} \ \ \ \ \ (9)$

Substituting into 8 leads to the recursion relation

$\displaystyle c_{q+2}=\frac{\rho_{0}(2q+2l+3)-1}{(q+2)(q+2l+3)}c_{q} \ \ \ \ \ (10)$

with ${c_{1}=0}$, so that ${c_{q}=0}$ for all odd ${q}$. The requirement that the series terminates at some finite value of ${j}$ leads to the quantization condition on ${E}$:

 $\displaystyle E$ $\displaystyle =$ $\displaystyle \hbar\omega\left(q_{max}+l+\frac{3}{2}\right) \ \ \ \ \ (11)$

or, defining ${n=q_{max}+l}$,

$\displaystyle E_{n}=\hbar\omega\left(n+\frac{3}{2}\right) \ \ \ \ \ (12)$

We worked out the degeneracies in the earlier post as well, so that the degeneracy of ${E_{n}}$ is

$\displaystyle d\left(n\right)=\frac{1}{2}(n+1)(n+2) \ \ \ \ \ (13)$

To complete Shankar’s exercise, we need to work out the eigenfunctions for ${n=0}$ and ${n=1}$. For ${n=0}$, ${q_{max}=l=0}$, so only ${c_{0}\ne0}$ and we have

 $\displaystyle v\left(\rho\right)$ $\displaystyle =$ $\displaystyle c_{0}\ \ \ \ \ (14)$ $\displaystyle u\left(\rho\right)$ $\displaystyle =$ $\displaystyle c_{0}\rho e^{-\rho_{0}\rho^{2}/2}\ \ \ \ \ (15)$ $\displaystyle \psi_{000}$ $\displaystyle =$ $\displaystyle \frac{u}{r}Y_{0}^{0}\ \ \ \ \ (16)$ $\displaystyle$ $\displaystyle =$ $\displaystyle c_{0}\kappa e^{-\rho_{0}\rho^{2}/2}Y_{0}^{0}\ \ \ \ \ (17)$ $\displaystyle$ $\displaystyle =$ $\displaystyle c_{0}\sqrt{\frac{2\mu3\omega}{4\pi\hbar}}e^{-\mu\omega r^{2}/2\hbar} \ \ \ \ \ (18)$

where in the fourth line we used

 $\displaystyle \kappa$ $\displaystyle =$ $\displaystyle \frac{\sqrt{2\mu E}}{\hbar}=\frac{\sqrt{2\mu\frac{3}{2}\hbar\omega}}{\hbar}=\sqrt{\frac{3\mu\omega}{\hbar}}\ \ \ \ \ (19)$ $\displaystyle Y_{0}^{0}$ $\displaystyle =$ $\displaystyle \frac{1}{\sqrt{4\pi}} \ \ \ \ \ (20)$

Normalizing this requires that

 $\displaystyle \int_{0}^{2\pi}\int_{0}^{\pi}\int_{0}^{\infty}\psi_{000}^{2}r^{2}\sin\theta dr\;d\theta\;d\phi$ $\displaystyle =$ $\displaystyle c_{0}^{2}\frac{6\mu\omega}{\hbar}\int_{0}^{\infty}e^{-\mu\omega r^{2}/\hbar}r^{2}dr\ \ \ \ \ (21)$ $\displaystyle$ $\displaystyle =$ $\displaystyle 1 \ \ \ \ \ (22)$

This is a standard Gaussian integral and can be done using software or tables so we get

$\displaystyle c_{0}=\frac{\sqrt{6}}{3}\left(\frac{\mu\omega}{\pi\hbar}\right)^{1/4} \ \ \ \ \ (23)$

This gives a wave function of

$\displaystyle \psi_{000}=\left(\frac{\mu\omega}{\pi\hbar}\right)^{3/4}e^{-\mu\omega r^{2}/2\hbar} \ \ \ \ \ (24)$

which agrees with the earlier result.

For ${n=1}$, the degeneracy is, from 13

$\displaystyle d\left(1\right)=3 \ \ \ \ \ (25)$

The three possibilities are ${m=0,\pm1}$ which are reflected in the three spherical harmonics ${Y_{1}^{0,\pm1}}$. The radial function is the same in all cases, and is obtained from ${q_{max}=0}$, ${l=1}$. From 7, this gives

 $\displaystyle v\left(\rho\right)$ $\displaystyle =$ $\displaystyle c_{0}\ \ \ \ \ (26)$ $\displaystyle u\left(\rho\right)$ $\displaystyle =$ $\displaystyle c_{0}\rho^{2}e^{-\rho_{0}\rho^{2}/2}\ \ \ \ \ (27)$ $\displaystyle \psi_{11m}$ $\displaystyle =$ $\displaystyle \frac{u}{r}Y_{1}^{m}\ \ \ \ \ (28)$ $\displaystyle$ $\displaystyle =$ $\displaystyle c_{0}\kappa^{2}re^{-\rho_{0}\rho^{2}/2}Y_{1}^{m}\ \ \ \ \ (29)$ $\displaystyle$ $\displaystyle =$ $\displaystyle c_{0}\frac{5\mu\omega}{\hbar}re^{-\mu\omega r^{2}/2\hbar}Y_{1}^{m} \ \ \ \ \ (30)$

Again, we work out ${c_{0}}$ by imposing normalization. For example

 $\displaystyle \psi_{111}$ $\displaystyle =$ $\displaystyle c_{0}\frac{5\mu\omega}{\hbar}re^{-\mu\omega r^{2}/2\hbar}Y_{1}^{1}\ \ \ \ \ (31)$ $\displaystyle$ $\displaystyle =$ $\displaystyle -c_{0}\frac{5\mu\omega}{\hbar}re^{-\mu\omega r^{2}/2\hbar}\sqrt{\frac{3}{8\pi}}\sin\theta e^{i\phi} \ \ \ \ \ (32)$

The normalization integral is

 $\displaystyle \int_{0}^{2\pi}\int_{0}^{\pi}\int_{0}^{\infty}\psi_{111}^{2}r^{2}\sin\theta dr\;d\theta\;d\phi$ $\displaystyle =$ $\displaystyle c_{0}^{2}\left(\frac{5\mu\omega}{\hbar}\right)^{2}\frac{3}{8\pi}2\pi\int_{0}^{\pi}\int_{0}^{\infty}e^{-\mu\omega r^{2}/\hbar}r^{4}\sin^{3}\theta dr\;d\theta\ \ \ \ \ (33)$ $\displaystyle$ $\displaystyle =$ $\displaystyle c_{0}^{2}\frac{75}{8}\sqrt{\frac{\pi\hbar}{\mu\omega}}=1\ \ \ \ \ (34)$ $\displaystyle c_{0}$ $\displaystyle =$ $\displaystyle \frac{2\sqrt{6}}{15}\left(\frac{\mu\omega}{\pi\hbar}\right)^{1/4} \ \ \ \ \ (35)$

I used Maple to do the integrals. This gives a wave function of

$\displaystyle \psi_{111}=-\sqrt{\frac{\mu\omega}{\hbar}}\left(\frac{\mu\omega}{\pi\hbar}\right)^{3/4}re^{-\mu\omega r^{2}/2\hbar}\sin\theta e^{i\phi} \ \ \ \ \ (36)$

We can work out the other two wave functions the same way (I used Maple, so I won’t go into the details):

 $\displaystyle \psi_{11-1}$ $\displaystyle =$ $\displaystyle \sqrt{\frac{\mu\omega}{\hbar}}\left(\frac{\mu\omega}{\pi\hbar}\right)^{3/4}re^{-\mu\omega r^{2}/2\hbar}\sin\theta e^{-i\phi}\ \ \ \ \ (37)$ $\displaystyle \psi_{110}$ $\displaystyle =$ $\displaystyle \sqrt{\frac{2\mu\omega}{\hbar}}\left(\frac{\mu\omega}{\pi\hbar}\right)^{3/4}re^{-\mu\omega r^{2}/2\hbar}\cos\theta \ \ \ \ \ (38)$

The ${\psi_{110}}$ here is the same as ${\psi_{001}}$ in our rectangular solution set. The other two are linear combinations of ${\psi_{100}}$ and ${\psi_{010}}$ from our rectangular set, which were (the suffixes in these 2 equations stand for ${x}$, ${y}$ and ${z}$, and not ${n}$, ${l}$ and ${m}$):

 $\displaystyle \psi_{100}$ $\displaystyle =$ $\displaystyle \sqrt{\frac{2m\omega}{\hbar}}\left(\frac{m\omega}{\pi\hbar}\right)^{3/4}e^{-m\omega r^{2}/2\hbar}r\sin\theta\cos\phi\ \ \ \ \ (39)$ $\displaystyle \psi_{010}$ $\displaystyle =$ $\displaystyle \sqrt{\frac{2m\omega}{\hbar}}\left(\frac{m\omega}{\pi\hbar}\right)^{3/4}e^{-m\omega r^{2}/2\hbar}r\sin\theta\sin\phi \ \ \ \ \ (40)$

We have

 $\displaystyle \psi_{111}$ $\displaystyle =$ $\displaystyle \frac{1}{\sqrt{2}}\left(\psi_{100}+i\psi_{010}\right)\ \ \ \ \ (41)$ $\displaystyle \psi_{11-1}$ $\displaystyle =$ $\displaystyle \frac{1}{\sqrt{2}}\left(\psi_{100}-i\psi_{010}\right) \ \ \ \ \ (42)$

# Harmonic oscillator in 2-d and 3-d, and in polar and spherical coordinates

Shankar, R. (1994), Principles of Quantum Mechanics, Plenum Press. Chapter 10, Exercises 10.2.2 – 10.2.3.

We’ve seen that the 3-d isotropic harmonic oscillator can be solved in rectangular coordinates using separation of variables. The Hamiltonian is

$\displaystyle H=\frac{p_{x}^{2}+p_{y}^{2}+p_{z}^{2}}{2m}+\frac{m\omega^{2}}{2}\left(x^{2}+y^{2}+z^{2}\right) \ \ \ \ \ (1)$

The solution to the Schrödinger equation is just the product of three one-dimensional oscillator eigenfunctions, one for each coordinate. That is

$\displaystyle \psi_{n}\left(x,y,z\right)=\psi_{n_{x}}\left(x\right)\psi_{n_{y}}\left(y\right)\psi_{n_{z}}\left(z\right) \ \ \ \ \ (2)$

Each one-dimensional eigenfunction can be expressed in terms of Hermite polynomials as

$\displaystyle \psi_{n_{x}}\left(x\right)=\left(\frac{m\omega}{\pi\hbar}\right)^{1/4}\frac{1}{\sqrt{2^{n_{x}}n_{x}!}}H_{n_{x}}\left(\sqrt{\frac{m\omega}{\hbar}}x\right)e^{-m\omega x^{2}/2\hbar} \ \ \ \ \ (3)$

with the functions for ${y}$ and ${z}$ obtained by replacing ${x}$ by ${y}$ or ${z}$ and ${n_{x}}$ by ${n_{y}}$ or ${n_{z}}$. We also saw earlier that in the 3-d oscillator, the total energy for state ${\psi_{n}\left(x,y,z\right)}$ is given in terms of the quantum numbers of the three 1-d oscillators as

$\displaystyle E_{n}=\hbar\omega\left(n+\frac{3}{2}\right)=\hbar\omega\left(n_{x}+n_{y}+n_{z}+\frac{3}{2}\right) \ \ \ \ \ (4)$

and that the degeneracy of level ${n}$ is ${\frac{1}{2}\left(n+1\right)\left(n+2\right)}$.

Since the Hermite polynomial ${H_{n_{x}}}$ has parity ${\left(-1\right)^{n_{x}}}$ (that is, odd (even) polynomials are odd (even) functions), the 3-d wave function ${\psi_{n}}$ has parity ${\left(-1\right)^{n_{x}}\left(-1\right)^{n_{y}}\left(-1\right)^{n_{z}}=\left(-1\right)^{n}}$.

We can write the one ${n=0}$ state and three ${n=1}$ states in spherical coordinates using the standard transformation

 $\displaystyle x$ $\displaystyle =$ $\displaystyle r\sin\theta\cos\phi\ \ \ \ \ (5)$ $\displaystyle y$ $\displaystyle =$ $\displaystyle r\sin\theta\sin\phi\ \ \ \ \ (6)$ $\displaystyle z$ $\displaystyle =$ $\displaystyle r\cos\phi \ \ \ \ \ (7)$

Using the notation ${\psi_{n}=\psi_{n_{x}n_{y}n_{z}}=\psi_{n_{x}}\psi_{n_{y}}\psi_{n_{z}}}$, we have, using ${H_{0}\left(y\right)=1}$ and ${H_{1}\left(y\right)=2y}$:

 $\displaystyle \psi_{000}$ $\displaystyle =$ $\displaystyle \left(\frac{m\omega}{\pi\hbar}\right)^{3/4}e^{-m\omega r^{2}/2\hbar}\ \ \ \ \ (8)$ $\displaystyle \psi_{100}$ $\displaystyle =$ $\displaystyle \sqrt{\frac{2m\omega}{\hbar}}\left(\frac{m\omega}{\pi\hbar}\right)^{3/4}e^{-m\omega r^{2}/2\hbar}r\sin\theta\cos\phi\ \ \ \ \ (9)$ $\displaystyle \psi_{010}$ $\displaystyle =$ $\displaystyle \sqrt{\frac{2m\omega}{\hbar}}\left(\frac{m\omega}{\pi\hbar}\right)^{3/4}e^{-m\omega r^{2}/2\hbar}r\sin\theta\sin\phi\ \ \ \ \ (10)$ $\displaystyle \psi_{001}$ $\displaystyle =$ $\displaystyle \sqrt{\frac{2m\omega}{\hbar}}\left(\frac{m\omega}{\pi\hbar}\right)^{3/4}e^{-m\omega r^{2}/2\hbar}r\cos\theta \ \ \ \ \ (11)$

We can check that these are the correct spherical versions of the eigenfunctions by using the Schrödinger equation in spherical coordinates, which is

$\displaystyle H\psi=\left[-\frac{\hbar^{2}\nabla^{2}}{2m}+\frac{m\omega^{2}}{2}r^{2}\right]\psi=E\psi \ \ \ \ \ (12)$

The spherical laplacian operator is

$\displaystyle \nabla^{2}\psi=\frac{1}{r^{2}}\frac{\partial}{\partial r}\left(r^{2}\frac{\partial\psi}{\partial r}\right)+\frac{1}{r^{2}\sin\theta}\frac{\partial}{\partial\theta}\left(\sin\theta\frac{\partial\psi}{\partial\theta}\right)+\frac{1}{r^{2}\sin^{2}\theta}\frac{\partial^{2}\psi}{\partial\phi^{2}} \ \ \ \ \ (13)$

You can grind through the derivatives by hand if you like, but I just used Maple to do it, giving the results

 $\displaystyle H\psi_{000}$ $\displaystyle =$ $\displaystyle \frac{3}{2}\hbar\omega\psi_{000}\ \ \ \ \ (14)$ $\displaystyle H\psi_{100}$ $\displaystyle =$ $\displaystyle \frac{5}{2}\hbar\omega\psi_{100}\ \ \ \ \ (15)$ $\displaystyle H\psi_{010}$ $\displaystyle =$ $\displaystyle \frac{5}{2}\hbar\omega\psi_{010}\ \ \ \ \ (16)$ $\displaystyle H\psi_{001}$ $\displaystyle =$ $\displaystyle \frac{5}{2}\hbar\omega\psi_{001} \ \ \ \ \ (17)$

In two dimensions, the analysis is pretty much the same. In the more general case where the masses are equal, but ${\omega_{x}\ne\omega_{y}}$, the Hamiltonian is

$\displaystyle H=\frac{p_{x}^{2}+p_{y}^{2}}{2m}+\frac{m}{2}\left(\omega_{x}^{2}x^{2}+\omega_{y}^{2}y^{2}\right) \ \ \ \ \ (18)$

A solution by separation of variables still works, with the result

$\displaystyle \psi_{n}\left(x,y\right)=\psi_{n_{x}}\left(x\right)\psi_{n_{y}}\left(y\right) \ \ \ \ \ (19)$

The total energy is

$\displaystyle E_{n}=E_{n_{x}}+E_{n_{y}}=\hbar\omega\left(n_{x}+\frac{1}{2}+n_{y}+\frac{1}{2}\right)=\hbar\omega\left(n+1\right) \ \ \ \ \ (20)$

For a given energy level ${n=n_{x}+n_{y}}$, there are ${n+1}$ ways of forming ${n}$ out of a sum of 2 non-negative integers, so the degeneracy of level ${n}$ is ${n+1}$.

The one ${n=0}$ state and two ${n=1}$ states are

 $\displaystyle \psi_{00}$ $\displaystyle =$ $\displaystyle \left(\frac{m\omega}{\pi\hbar}\right)^{1/2}e^{-m\omega\left(x^{2}+y^{2}\right)/2\hbar}\ \ \ \ \ (21)$ $\displaystyle \psi_{10}$ $\displaystyle =$ $\displaystyle \sqrt{\frac{2m\omega}{\hbar}}\left(\frac{m\omega}{\pi\hbar}\right)^{1/2}e^{-m\omega\left(x^{2}+y^{2}\right)/2\hbar}x\ \ \ \ \ (22)$ $\displaystyle \psi_{01}$ $\displaystyle =$ $\displaystyle \sqrt{\frac{2m\omega}{\hbar}}\left(\frac{m\omega}{\pi\hbar}\right)^{1/2}e^{-m\omega\left(x^{2}+y^{2}\right)/2\hbar}y \ \ \ \ \ (23)$

To translate to polar coordinates, we use the transformations

 $\displaystyle x$ $\displaystyle =$ $\displaystyle \rho\cos\phi\ \ \ \ \ (24)$ $\displaystyle y$ $\displaystyle =$ $\displaystyle \rho\sin\phi \ \ \ \ \ (25)$

so we have

 $\displaystyle \psi_{00}$ $\displaystyle =$ $\displaystyle \left(\frac{m\omega}{\pi\hbar}\right)^{1/2}e^{-m\omega\rho^{2}/2\hbar}\ \ \ \ \ (26)$ $\displaystyle \psi_{10}$ $\displaystyle =$ $\displaystyle \sqrt{\frac{2m\omega}{\hbar}}\left(\frac{m\omega}{\pi\hbar}\right)^{1/2}e^{-m\omega\rho^{2}/2\hbar}\rho\cos\phi\ \ \ \ \ (27)$ $\displaystyle \psi_{01}$ $\displaystyle =$ $\displaystyle \sqrt{\frac{2m\omega}{\hbar}}\left(\frac{m\omega}{\pi\hbar}\right)^{1/2}e^{-m\omega\rho^{2}/2\hbar}\rho\sin\phi \ \ \ \ \ (28)$

Again, we can check this by plugging these polar formulas into the polar Schrödinger equation, where the 2-d Laplacian is

$\displaystyle \nabla^{2}=\frac{\partial^{2}}{\partial\rho^{2}}+\frac{1}{\rho}\frac{\partial}{\partial\rho}+\frac{1}{\rho^{2}}\frac{\partial^{2}}{\partial\phi^{2}} \ \ \ \ \ (29)$

The results are

 $\displaystyle H\psi_{00}$ $\displaystyle =$ $\displaystyle \hbar\omega\psi_{00}\ \ \ \ \ (30)$ $\displaystyle H\psi_{10}$ $\displaystyle =$ $\displaystyle 2\hbar\omega\psi_{10}\ \ \ \ \ (31)$ $\displaystyle H\psi_{01}$ $\displaystyle =$ $\displaystyle 2\hbar\omega\psi_{01} \ \ \ \ \ (32)$

# Creation and annihilation operators for the 3-d harmonic oscillator

Reference: Tom Lancaster and Stephen J. Blundell, Quantum Field Theory for the Gifted Amateur, (Oxford University Press, 2014), Problem 3.3.

For the 3-d harmonic oscillator, we can write the hamiltonian as

$\displaystyle \hat{H}=\frac{1}{2m}\left(\hat{p}_{1}^{2}+\hat{p}_{2}^{2}+\hat{p}_{3}^{2}\right)+\frac{m\omega^{2}}{2}\left(\hat{x}_{1}^{2}+\hat{x}_{2}^{2}+\hat{x}_{3}^{2}\right) \ \ \ \ \ (1)$

This is effectively three independent oscillators, one in each of the three coordinate directions, so using the 1-d form of the hamiltonian in terms of creation and annihilation operators, we have

$\displaystyle \hat{H}=\hbar\omega\sum_{i=1}^{3}\left(\frac{1}{2}+\hat{a}_{i}^{\dagger}\hat{a}_{i}\right) \ \ \ \ \ (2)$

Using the position and momentum operators expressed as

 $\displaystyle \hat{x}_{i}$ $\displaystyle =$ $\displaystyle \sqrt{\frac{\hbar}{2m\omega}}\left(\hat{a}_{i}^{\dagger}+\hat{a}_{i}\right)\ \ \ \ \ (3)$ $\displaystyle \hat{p}_{i}$ $\displaystyle =$ $\displaystyle i\sqrt{\frac{\hbar m\omega}{2}}\left(\hat{a}_{i}^{\dagger}-\hat{a}_{i}\right) \ \ \ \ \ (4)$

we can write the components of angular momentum in terms of creation and annihilation operators. For example, the ${z}$ component can be written as

 $\displaystyle L^{3}$ $\displaystyle =$ $\displaystyle \hat{x}_{1}\hat{p}_{2}-\hat{x}_{2}\hat{p}_{1}\ \ \ \ \ (5)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{i\hbar}{2}\left[\left(\hat{a}_{1}^{\dagger}+\hat{a}_{1}\right)\left(\hat{a}_{2}^{\dagger}-\hat{a}_{2}\right)-\left(\hat{a}_{2}^{\dagger}+\hat{a}_{2}\right)\left(\hat{a}_{1}^{\dagger}-\hat{a}_{1}\right)\right]\ \ \ \ \ (6)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{i\hbar}{2}\left[-\hat{a}_{1}^{\dagger}\hat{a}_{2}+\hat{a}_{1}\hat{a}_{2}^{\dagger}-\hat{a}_{2}\hat{a}_{1}^{\dagger}+\hat{a}_{2}^{\dagger}\hat{a}_{1}\right]\ \ \ \ \ (7)$ $\displaystyle$ $\displaystyle =$ $\displaystyle -i\hbar\left(\hat{a}_{1}^{\dagger}\hat{a}_{2}-\hat{a}_{2}^{\dagger}\hat{a}_{1}\right) \ \ \ \ \ (8)$

where we’ve used the fact that all operators of one coordinate commute with all operators of the other two coordinates. If we work out the other two coordinates we get the general formula

$\displaystyle L^{i}=-i\hbar\epsilon^{ijk}\hat{a}_{j}^{\dagger}\hat{a}_{k} \ \ \ \ \ (9)$

where ${\epsilon^{ijk}}$ is the Levi-Civita symbol, which is +1 if ${ijk}$ is a even permutation of 1,2,3, ${-1}$ if ${ijk}$ is an odd permutation of 1,2,3 and 0 if any two of ${ijk}$ are equal, and repeated indices are summed.

We can define some new creation and annihilation operators as follows:

 $\displaystyle \hat{b}_{1}^{\dagger}$ $\displaystyle \equiv$ $\displaystyle -\frac{1}{\sqrt{2}}\left(\hat{a}_{1}^{\dagger}+i\hat{a}_{2}^{\dagger}\right)\ \ \ \ \ (10)$ $\displaystyle \hat{b}_{0}^{\dagger}$ $\displaystyle \equiv$ $\displaystyle \hat{a}_{3}^{\dagger}\ \ \ \ \ (11)$ $\displaystyle \hat{b}_{-1}^{\dagger}$ $\displaystyle \equiv$ $\displaystyle \frac{1}{\sqrt{2}}\left(\hat{a}_{1}^{\dagger}-i\hat{a}_{2}^{\dagger}\right) \ \ \ \ \ (12)$
 $\displaystyle \hat{b}_{1}$ $\displaystyle \equiv$ $\displaystyle -\frac{1}{\sqrt{2}}\left(\hat{a}_{1}-i\hat{a}_{2}\right)\ \ \ \ \ (13)$ $\displaystyle \hat{b}_{0}$ $\displaystyle \equiv$ $\displaystyle \hat{a}_{3}\ \ \ \ \ (14)$ $\displaystyle \hat{b}_{-1}$ $\displaystyle \equiv$ $\displaystyle \frac{1}{\sqrt{2}}\left(\hat{a}_{1}+i\hat{a}_{2}\right) \ \ \ \ \ (15)$

From the commutation relations for ${\hat{a}_{i}^{\dagger}}$ and ${\hat{a}_{i}}$:

$\displaystyle \left[\hat{a}_{i},\hat{a}_{j}^{\dagger}\right]=\delta_{ij} \ \ \ \ \ (16)$

we get

 $\displaystyle \left[\hat{b}_{0},\hat{b}_{j}^{\dagger}\right]$ $\displaystyle =$ $\displaystyle \delta_{0,j}\ \ \ \ \ (17)$ $\displaystyle \left[\hat{b}_{1},\hat{b}_{1}^{\dagger}\right]$ $\displaystyle =$ $\displaystyle \frac{1}{2}\left(\left[\hat{a}_{1},\hat{a}_{1}^{\dagger}\right]+\left(-i\right)i\left[\hat{a}_{2},\hat{a}_{2}^{\dagger}\right]\right)\ \ \ \ \ (18)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{1}{2}\left(1+1\right)\ \ \ \ \ (19)$ $\displaystyle$ $\displaystyle =$ $\displaystyle 1\ \ \ \ \ (20)$ $\displaystyle \left[\hat{b}_{1},\hat{b}_{-1}^{\dagger}\right]$ $\displaystyle =$ $\displaystyle -\frac{1}{2}\left(\left[\hat{a}_{1},\hat{a}_{1}^{\dagger}\right]+i^{2}\left[\hat{a}_{2},\hat{a}_{2}^{\dagger}\right]\right)\ \ \ \ \ (21)$ $\displaystyle$ $\displaystyle =$ $\displaystyle 0\ \ \ \ \ (22)$ $\displaystyle \left[\hat{b}_{-1},\hat{b}_{-1}^{\dagger}\right]$ $\displaystyle =$ $\displaystyle \frac{1}{2}\left(\left[\hat{a}_{1},\hat{a}_{1}^{\dagger}\right]+\left(-i\right)i\left[\hat{a}_{2},\hat{a}_{2}^{\dagger}\right]\right)\ \ \ \ \ (23)$ $\displaystyle$ $\displaystyle =$ $\displaystyle 1\ \ \ \ \ (24)$ $\displaystyle \left[\hat{b}_{-1},\hat{b}_{1}^{\dagger}\right]$ $\displaystyle =$ $\displaystyle -\frac{1}{2}\left(\left[\hat{a}_{1},\hat{a}_{1}^{\dagger}\right]+i^{2}\left[\hat{a}_{2},\hat{a}_{2}^{\dagger}\right]\right)\ \ \ \ \ (25)$ $\displaystyle$ $\displaystyle =$ $\displaystyle 0 \ \ \ \ \ (26)$

so in general

$\displaystyle \left[\hat{b}_{i},\hat{b}_{j}^{\dagger}\right]=\delta_{ij} \ \ \ \ \ (27)$

To express the hamiltonian 2 in terms of the new operators, we need

 $\displaystyle \hat{b}_{1}^{\dagger}\hat{b}_{1}$ $\displaystyle =$ $\displaystyle \frac{1}{2}\left[\hat{a}_{1}^{\dagger}\hat{a}_{1}+\hat{a}_{2}^{\dagger}\hat{a}_{2}+i\left(\hat{a}_{2}^{\dagger}\hat{a}_{1}-\hat{a}_{1}^{\dagger}\hat{a}_{2}\right)\right]\ \ \ \ \ (28)$ $\displaystyle \hat{b}_{-1}^{\dagger}\hat{b}_{-1}$ $\displaystyle =$ $\displaystyle \frac{1}{2}\left[\hat{a}_{1}^{\dagger}\hat{a}_{1}+\hat{a}_{2}^{\dagger}\hat{a}_{2}-i\left(\hat{a}_{2}^{\dagger}\hat{a}_{1}-\hat{a}_{1}^{\dagger}\hat{a}_{2}\right)\right]\ \ \ \ \ (29)$ $\displaystyle \hat{b}_{1}^{\dagger}\hat{b}_{1}+\hat{b}_{-1}^{\dagger}\hat{b}_{-1}$ $\displaystyle =$ $\displaystyle \hat{a}_{1}^{\dagger}\hat{a}_{1}+\hat{a}_{2}^{\dagger}\hat{a}_{2}\ \ \ \ \ (30)$ $\displaystyle \hat{H}$ $\displaystyle =$ $\displaystyle \hbar\omega\sum_{i=-1}^{+1}\left(\frac{1}{2}+\hat{b}_{i}^{\dagger}\hat{b}_{i}\right) \ \ \ \ \ (31)$

Finally, we can write ${\hat{L}^{3}}$ in terms of the new operators.

 $\displaystyle -\hat{b}_{-1}^{\dagger}\hat{b}_{-1}+0\times\hat{b}_{0}^{\dagger}\hat{b}_{0}+\hat{b}_{1}^{\dagger}\hat{b}_{1}$ $\displaystyle =$ $\displaystyle i\left(\hat{a}_{2}^{\dagger}\hat{a}_{1}-\hat{a}_{1}^{\dagger}\hat{a}_{2}\right)\ \ \ \ \ (32)$ $\displaystyle \hat{L}^{3}$ $\displaystyle =$ $\displaystyle \hbar\sum_{m=-1}^{+1}m\hat{b}_{m}^{\dagger}\hat{b}_{m} \ \ \ \ \ (33)$

# Rubber band helium

References: Griffiths, David J. (2005), Introduction to Quantum Mechanics, 2nd Edition; Pearson Education – Problem 7.17.

We’ve looked at the helium atom using the variational principle. Although the helium atom using the correct Coulomb potential cannot be solved exactly, a variant known as ‘rubber band helium’ can be. Here we use simple harmonic oscillator potentials. The hamiltonian is then:

$\displaystyle H=-\frac{\hbar^{2}}{2m}\left(\nabla_{1}^{2}+\nabla_{2}^{2}\right)+\frac{1}{2}m\omega^{2}\left(r_{1}^{2}+r_{2}^{2}\right)-\frac{\lambda}{4}m\omega^{2}\left|\mathbf{r}_{1}-\mathbf{r}_{2}\right|^{2} \ \ \ \ \ (1)$

By introducing a change of variables, we can decouple the hamiltonian. Let

 $\displaystyle \mathbf{u}$ $\displaystyle \equiv$ $\displaystyle \frac{1}{\sqrt{2}}\left(\mathbf{r}_{1}+\mathbf{r}_{2}\right)\ \ \ \ \ (2)$ $\displaystyle \mathbf{v}$ $\displaystyle \equiv$ $\displaystyle \frac{1}{\sqrt{2}}\left(\mathbf{r}_{1}-\mathbf{r}_{2}\right) \ \ \ \ \ (3)$

The gradient operators then transform as

 $\displaystyle \nabla_{u}$ $\displaystyle =$ $\displaystyle \frac{1}{\sqrt{2}}\left(\nabla_{1}+\nabla_{2}\right)\ \ \ \ \ (4)$ $\displaystyle \nabla_{v}$ $\displaystyle =$ $\displaystyle \frac{1}{\sqrt{2}}\left(\nabla_{1}-\nabla_{2}\right)\ \ \ \ \ (5)$ $\displaystyle \nabla_{u}^{2}$ $\displaystyle =$ $\displaystyle \frac{1}{2}\left(\nabla_{1}^{2}+\nabla_{2}^{2}+2\nabla_{1}\cdot\nabla_{2}\right)\ \ \ \ \ (6)$ $\displaystyle \nabla_{v}^{2}$ $\displaystyle =$ $\displaystyle \frac{1}{2}\left(\nabla_{1}^{2}+\nabla_{2}^{2}-2\nabla_{1}\cdot\nabla_{2}\right)\ \ \ \ \ (7)$ $\displaystyle \nabla_{u}^{2}+\nabla_{v}^{2}$ $\displaystyle =$ $\displaystyle \nabla_{1}^{2}+\nabla_{2}^{2} \ \ \ \ \ (8)$

For the potential terms, we have

 $\displaystyle u^{2}+v^{2}$ $\displaystyle =$ $\displaystyle \frac{1}{2}\left[r_{1}^{2}+r_{2}^{2}+2\mathbf{r}_{1}\cdot\mathbf{r}_{2}+r_{1}^{2}+r_{2}^{2}-2\mathbf{r}_{1}\cdot\mathbf{r}_{2}\right]\ \ \ \ \ (9)$ $\displaystyle$ $\displaystyle =$ $\displaystyle r_{1}^{2}+r_{2}^{2}\ \ \ \ \ (10)$ $\displaystyle \left|\mathbf{r}_{1}-\mathbf{r}_{2}\right|^{2}$ $\displaystyle =$ $\displaystyle r_{1}^{2}+r_{2}^{2}-2\mathbf{r}_{1}\cdot\mathbf{r}_{2}\ \ \ \ \ (11)$ $\displaystyle$ $\displaystyle =$ $\displaystyle 2v^{2} \ \ \ \ \ (12)$

Thus the hamiltonian separates:

$\displaystyle H=\left[-\frac{\hbar^{2}}{2m}\nabla_{u}^{2}+\frac{1}{2}m\omega^{2}u^{2}\right]+\left[-\frac{\hbar^{2}}{2m}\nabla_{v}^{2}+\frac{1}{2}m\omega^{2}\left(1-\lambda\right)v^{2}\right] \ \ \ \ \ (13)$

which is the sum of two 3-d harmonic oscillators. The exact ground state energy of this system are then just the sum of the two separate oscillator energies:

$\displaystyle E_{0}=\frac{3}{2}\hbar\omega+\frac{3}{2}\hbar\omega\sqrt{1-\lambda} \ \ \ \ \ (14)$

To test the variational principle for this potential, we can start with the (known) ground state wave function for the 3-d harmonic oscillator as the test function.

$\displaystyle \psi=\left(\frac{m\omega}{\pi\hbar}\right)^{3/2}e^{-m\omega\left(r_{1}^{2}+r_{2}^{2}\right)/2\hbar} \ \ \ \ \ (15)$

This function is an eigenfunction of the first two terms in 1 with energy ${3\hbar\omega}$ so we have

$\displaystyle \left\langle H\right\rangle =3\hbar\omega+\left\langle V_{\lambda}\right\rangle \ \ \ \ \ (16)$

where

 $\displaystyle \left\langle V_{\lambda}\right\rangle$ $\displaystyle =$ $\displaystyle -\frac{\lambda}{4}m\omega^{2}\left(\frac{m\omega}{\pi\hbar}\right)^{3}\int e^{-m\omega\left(r_{1}^{2}+r_{2}^{2}\right)/\hbar}\left|\mathbf{r}_{1}-\mathbf{r}_{2}\right|^{2}d^{3}\mathbf{r}_{1}d^{3}\mathbf{r}_{2}\ \ \ \ \ (17)$ $\displaystyle$ $\displaystyle =$ $\displaystyle -\frac{\lambda}{4}m\omega^{2}\left(\frac{m\omega}{\pi\hbar}\right)^{3}\int e^{-m\omega\left(r_{1}^{2}+r_{2}^{2}\right)/\hbar}\left(r_{1}^{2}+r_{2}^{2}-2r_{1}r_{2}\cos\theta_{2}\right)d^{3}\mathbf{r}_{1}d^{3}\mathbf{r}_{2} \ \ \ \ \ (18)$

The term with ${\cos\theta_{2}}$ integrates to zero when we do the ${\theta_{2}}$ integral, so we’re left with two Gaussian integrals and we get

$\displaystyle \left\langle V_{\lambda}\right\rangle =-\frac{3}{4}\lambda\hbar\omega \ \ \ \ \ (19)$

Plugging this back into 16 we get

$\displaystyle \left\langle H\right\rangle =3\hbar\omega\left(1-\frac{\lambda}{4}\right) \ \ \ \ \ (20)$

This is actually the Taylor expansion with respect to ${\lambda}$ of 14 up to first order.

# Perturbing the 3-d harmonic oscillator

References: Griffiths, David J. (2005), Introduction to Quantum Mechanics, 2nd Edition; Pearson Education – Problem 6.30.

Here’s an example of perturbation theory applied to the 3-d harmonic oscillator. Using separation of variables, the stationary states in 3-d are just products of the 1-d states:

$\displaystyle \psi_{n}(x)=\left(\frac{m\omega}{\pi\hbar}\right)^{1/4}\frac{1}{\sqrt{2^{n}n!}}H_{n}\left(\sqrt{\frac{m\omega}{\hbar}}x\right)e^{-m\omega x^{2}/2\hbar} \ \ \ \ \ (1)$

where ${H_{n}}$ is the ${n}$th Hermite polynomial. The first two Hermite polynomials are

 $\displaystyle H_{0}\left(x\right)$ $\displaystyle =$ $\displaystyle 1\ \ \ \ \ (2)$ $\displaystyle H_{1}\left(x\right)$ $\displaystyle =$ $\displaystyle 2x \ \ \ \ \ (3)$

Thus the ground state in 3-d is (where the notation is ${\psi_{n_{x}n_{y}n_{z}}):}$

$\displaystyle \psi_{000}=\left(\frac{m\omega}{\pi\hbar}\right)^{3/4}e^{-m\omega x^{2}/2\hbar}e^{-m\omega y^{2}/2\hbar}e^{-m\omega z^{2}/2\hbar} \ \ \ \ \ (4)$

The first excited state is triply degenerate, since the excitation can be in any of the three coordinates. We get

 $\displaystyle \psi_{100}$ $\displaystyle =$ $\displaystyle \frac{2}{\pi^{3/4}}\left(\frac{m\omega}{\hbar}\right)^{5/4}xe^{-m\omega x^{2}/2\hbar}e^{-m\omega y^{2}/2\hbar}e^{-m\omega z^{2}/2\hbar}\ \ \ \ \ (5)$ $\displaystyle \psi_{010}$ $\displaystyle =$ $\displaystyle \frac{2}{\pi^{3/4}}\left(\frac{m\omega}{\hbar}\right)^{5/4}ye^{-m\omega x^{2}/2\hbar}e^{-m\omega y^{2}/2\hbar}e^{-m\omega z^{2}/2\hbar}\ \ \ \ \ (6)$ $\displaystyle \psi_{001}$ $\displaystyle =$ $\displaystyle \frac{2}{\pi^{3/4}}\left(\frac{m\omega}{\hbar}\right)^{5/4}ze^{-m\omega x^{2}/2\hbar}e^{-m\omega y^{2}/2\hbar}e^{-m\omega z^{2}/2\hbar} \ \ \ \ \ (7)$

Now suppose we introduce a perturbation given by

$\displaystyle H^{\prime}=\lambda x^{2}yz \ \ \ \ \ (8)$

for some constant ${\lambda}$. We can use non-degenerate perturbation theory to calculate the correction to the energy in the ground state. This gives

 $\displaystyle E_{000,1}$ $\displaystyle =$ $\displaystyle \left\langle \psi_{000}\right|H^{\prime}\left|\psi_{000}\right\rangle \ \ \ \ \ (9)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \lambda\int\int\int\left|\psi_{000}\left(x,y,z\right)\right|^{2}x^{2}yzdxdydz \ \ \ \ \ (10)$

The integrals over ${y}$ and ${z}$ are zero, since ${\psi_{000}}$ is an even function over all three coordinates and ${H^{\prime}}$ is odd in ${y}$ and ${z}$. Thus

$\displaystyle E_{000,1}=0 \ \ \ \ \ (11)$

and there is no correction to the ground state energy to first order.

For the first excited state, we need to use the degenerate theory which means we need to calculate the matrix elements ${W_{ab}=\left\langle n_{x}n_{y}n_{z}\right|H^{\prime}\left|m_{x}m_{y}m_{z}\right\rangle }$ for the 3 degenerate first excited states. These matrix elements are just the means of the coordinates or their squares, that is

$\displaystyle W_{ab}=\lambda\left\langle n_{x}\right|x^{2}\left|m_{x}\right\rangle \left\langle n_{y}\right|y\left|m_{y}\right\rangle \left\langle n_{z}\right|z\left|m_{z}\right\rangle \ \ \ \ \ (12)$

where the ${n}$s and ${m}$s must be chosen so that they make up one of the three degenerate states above. (That is, one each of the ${n}$s and ${m}$s must be 1 and the other two are 0.)

We’ve worked out some of these elements before, in particular in the 1-d case:

$\displaystyle \langle n|x|n'\rangle=\sqrt{\frac{\hbar}{2m\omega}}\left(\sqrt{n'+1}\delta_{n,n'+1}+\sqrt{n'}\delta_{n,n'-1}\right) \ \ \ \ \ (13)$

This shows that all the diagonal elements of ${W}$ are zero, so we need look only at the off-diagonal elements. The quantity ${\left\langle n_{x}\right|x^{2}\left|m_{x}\right\rangle =0}$ if ${n_{x}\ne m_{x}}$ since the integrand contains a factor of ${x^{3}}$ which is an odd function. Thus the only non-zero elements must have ${n_{x}=m_{x}}$, ${n_{y}\ne m_{y}}$ and ${n_{z}\ne m_{z}}$. The only elements that satisfy these conditions are ${W_{010,001}=\lambda\left\langle 0\right|x^{2}\left|0\right\rangle \left\langle 1\right|y\left|0\right\rangle \left\langle 0\right|z\left|1\right\rangle }$ and ${W_{001,010}=\lambda\left\langle 0\right|x^{2}\left|0\right\rangle \left\langle 0\right|y\left|1\right\rangle \left\langle 1\right|z\left|0\right\rangle }$. For the last two factors we can use the formula above, while for ${\left\langle 0\right|x^{2}\left|0\right\rangle }$ we can use an earlier result:

 $\displaystyle \left\langle 0\right|x^{2}\left|0\right\rangle$ $\displaystyle =$ $\displaystyle \frac{\hbar}{2m\omega}\ \ \ \ \ (14)$ $\displaystyle \left\langle 0\right|y\left|1\right\rangle$ $\displaystyle =$ $\displaystyle \sqrt{\frac{\hbar}{2m\omega}}\ \ \ \ \ (15)$ $\displaystyle \left\langle 1\right|z\left|0\right\rangle$ $\displaystyle =$ $\displaystyle \sqrt{\frac{\hbar}{2m\omega}} \ \ \ \ \ (16)$

Thus

$\displaystyle W=\lambda\left[\begin{array}{ccc} 0 & 0 & 0\\ 0 & 0 & \left(\frac{\hbar}{2m\omega}\right)^{2}\\ 0 & \left(\frac{\hbar}{2m\omega}\right)^{2} & 0 \end{array}\right] \ \ \ \ \ (17)$

The energy corrections are the eigenvalues of ${W}$, which are

$\displaystyle E_{1}=0,\pm\lambda\left(\frac{\hbar}{2m\omega}\right)^{2} \ \ \ \ \ (18)$

Thus the perturbation splits the degeneracy, leaving one energy the same and moving the other two on either side of the original energy.

# Statistical mechanics in quantum theory: 3-d harmonic oscillator

References: Griffiths, David J. (2005), Introduction to Quantum Mechanics, 2nd Edition; Pearson Education – Problem 5.37.

Here’s another example of working out the energy of a collection of particles. This time we’ll look at the 3-d harmonic oscillator and consider distinguishable particles only. In this case, the number of particles ${n_{j}}$ in energy level ${j}$ is

$\displaystyle n_{j}=d_{j}e^{-\alpha-\beta E_{j}} \ \ \ \ \ (1)$

where ${\alpha}$ and ${\beta}$ are the Lagrange multipliers and ${E_{j}}$ is the energy of that level. For the 3-d harmonic oscillator

$\displaystyle E_{j}=\left(j+\frac{3}{2}\right)\hbar\omega \ \ \ \ \ (2)$

and ${j}$ is the sum of the three quantum numbers ${j=j_{x}+j_{y}+j_{z}}$ in the three rectangular coordinates. The degeneracy of state ${j}$ was worked out to be

$\displaystyle d_{j}=\frac{1}{2}\left(j+1\right)\left(j+2\right) \ \ \ \ \ (3)$

The total number of particles is then

$\displaystyle N=\frac{e^{-\alpha-3\beta\hbar\omega/2}}{2}\sum_{j=0}^{\infty}\left(j+1\right)\left(j+2\right)e^{-\beta j\hbar\omega} \ \ \ \ \ (4)$

The sum can be evaluated directly in Maple, giving the result

 $\displaystyle N$ $\displaystyle =$ $\displaystyle e^{-\alpha-3\beta\hbar\omega/2}\frac{e^{3\beta\hbar\omega}}{\left(e^{\beta\hbar\omega}-1\right)^{3}}\ \ \ \ \ (5)$ $\displaystyle$ $\displaystyle =$ $\displaystyle e^{-\alpha}\frac{e^{3\beta\hbar\omega/2}}{\left(e^{\beta\hbar\omega}-1\right)^{3}}\ \ \ \ \ (6)$ $\displaystyle$ $\displaystyle =$ $\displaystyle e^{-\alpha}\frac{e^{-3\beta\hbar\omega/2}}{\left(1-e^{-\beta\hbar\omega}\right)^{3}} \ \ \ \ \ (7)$

However, to see how this is done using the method suggested by Griffiths, we start with the geometric series

$\displaystyle \sum_{j=0}^{\infty}x^{j}=\frac{1}{1-x} \ \ \ \ \ (8)$

Taking the first derivative:

 $\displaystyle \frac{1}{\left(1-x\right)^{2}}$ $\displaystyle =$ $\displaystyle \sum_{j=0}^{\infty}jx^{j-1}\ \ \ \ \ (9)$ $\displaystyle \frac{x}{\left(1-x\right)^{2}}$ $\displaystyle =$ $\displaystyle \sum_{j=0}^{\infty}jx^{j} \ \ \ \ \ (10)$

And the second derivative:

 $\displaystyle \frac{2}{\left(1-x\right)^{3}}$ $\displaystyle =$ $\displaystyle \sum_{j=0}^{\infty}j\left(j-1\right)x^{j-2}\ \ \ \ \ (11)$ $\displaystyle \frac{2x^{2}}{\left(1-x\right)^{3}}$ $\displaystyle =$ $\displaystyle \sum_{j=0}^{\infty}j\left(j-1\right)x^{j} \ \ \ \ \ (12)$

We can write the coefficient in 4 as

$\displaystyle \left(j+1\right)\left(j+2\right)=j\left(j-1\right)+4j+2 \ \ \ \ \ (13)$

and define

$\displaystyle x\equiv e^{-\beta\hbar\omega} \ \ \ \ \ (14)$

Then we have

 $\displaystyle \sum_{j=0}^{\infty}\left(j+1\right)\left(j+2\right)e^{-\beta j\hbar\omega}$ $\displaystyle =$ $\displaystyle \sum_{j=0}^{\infty}j\left(j-1\right)x^{j}+4\sum_{j=0}^{\infty}jx^{j}+2\sum_{j=0}^{\infty}jx^{j-1}\ \ \ \ \ (15)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{2x^{2}}{\left(1-x\right)^{3}}+4\frac{x}{\left(1-x\right)^{2}}+\frac{2}{1-x}\ \ \ \ \ (16)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{2}{\left(1-x\right)^{3}}\left[x^{2}+2x\left(1-x\right)+\left(1-x\right)^{2}\right]\ \ \ \ \ (17)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{2}{\left(1-x\right)^{3}} \ \ \ \ \ (18)$

Plugging this back into 4 gives the result above for ${N}$.

Using ${\beta=1/k_{B}T}$, we can solve for ${\alpha}$ to get

$\displaystyle e^{-\alpha}=N\frac{\left(1-e^{-\hbar\omega/k_{B}T}\right)^{3}}{e^{-3\hbar\omega/2k_{B}T}} \ \ \ \ \ (19)$

In terms of the chemical potential, we have ${\mu=-\alpha k_{B}T}$ so

$\displaystyle \mu=k_{B}T\left[\ln N+3\ln\left(1-e^{-\hbar\omega/k_{B}T}\right)+\frac{3\hbar\omega}{2k_{B}T}\right] \ \ \ \ \ (20)$

For the total energy, we have

 $\displaystyle E$ $\displaystyle =$ $\displaystyle \sum_{j=0}^{\infty}n_{j}E_{j}\ \ \ \ \ (21)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{\hbar\omega}{2}e^{-\alpha}\sum_{j=0}^{\infty}\left(j+1\right)\left(j+2\right)\left(j+\frac{3}{2}\right)e^{-\beta E_{j}} \ \ \ \ \ (22)$

This sum can be done in a similar way to the previous one (though we need another derivative of the geometric series). However, I can be lazy and get Maple to do the work, giving the result

 $\displaystyle E$ $\displaystyle =$ $\displaystyle \frac{3\hbar\omega}{2}e^{-\alpha-3\beta\hbar\omega/2}\frac{e^{3\beta\hbar\omega}\left(e^{\beta\hbar\omega}+1\right)}{\left(e^{\beta\hbar\omega}-1\right)^{4}}\ \ \ \ \ (23)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{3\hbar\omega}{2}e^{-\alpha-3\beta\hbar\omega/2}\frac{\left(1+e^{-\beta\hbar\omega}\right)}{\left(1-e^{-\beta\hbar\omega}\right)^{3}\left(1-e^{-\beta\hbar\omega}\right)} \ \ \ \ \ (24)$

Substituting for ${e^{-\alpha}}$ and ${\beta}$ from above, we get

$\displaystyle E=\frac{3\hbar\omega N}{2}\frac{\left(1+e^{-\hbar\omega/k_{B}T}\right)}{\left(1-e^{-\hbar\omega/k_{B}T}\right)} \ \ \ \ \ (25)$

In the limit of very low temperatures ${k_{B}T\ll\hbar\omega}$ and

$\displaystyle \mu\left(T\right)\rightarrow k_{B}T\ln N+\frac{3}{2}\hbar\omega\rightarrow\frac{3}{2}\hbar\omega \ \ \ \ \ (26)$

and

$\displaystyle E\rightarrow\frac{3}{2}\hbar\omega N \ \ \ \ \ (27)$

Thus all particles settle into the ground state, although even at absolute zero the energy is not zero.

At the other extreme, where ${k_{B}T\gg\hbar\omega}$

 $\displaystyle 1+e^{-\hbar\omega/k_{B}T}$ $\displaystyle \rightarrow$ $\displaystyle 2\ \ \ \ \ (28)$ $\displaystyle 1-e^{-\hbar\omega/k_{B}T}$ $\displaystyle \rightarrow$ $\displaystyle \frac{\hbar\omega}{k_{B}T}\ \ \ \ \ (29)$ $\displaystyle E$ $\displaystyle \rightarrow$ $\displaystyle 3Nk_{B}T \ \ \ \ \ (30)$

The equipartition theorem from classical statistical mechanics says that at thermal equilibrium, each degree of freedom in the system contributes ${\frac{1}{2}k_{B}T}$ to the total energy. In the 3-d harmonic oscillator, each particle has 3 translational degrees of freedom and 3 vibrational degrees of freedom making a total of 6, giving the energy above.