2.5. Harmonic oscillator

The time-independent Schrödinger equation for the harmonic oscillator is

(2.76)\[\hat{H}\psi(x) = \left[ -\frac{\hbar^2}{2\mu}\frac{\mathrm{d}^2}{\mathrm{d}x^2} + \frac{1}{2}\mu\omega_\mathrm{e}^2 x^2\right]\psi(x) = E \cdot \psi(x)\]

where \(x\) is the displacement from the equilibrium position, \(\mu\) is the reduced mass, and \(\omega_\mathrm{e}\) is the oscillator frequency. The force constant \(k\) in the potential-energy term is expressed as \(k = \mu\omega_\mathrm{e}^2\). The goal here is to solve this equation for \(\psi(x)\) and \(E\).

First, let us simplify the equation. We divide both sides by \(\hbar\omega_\mathrm{e}\) and multiply by 2. This gives

(2.77)\[\left[ -\frac{\hbar}{\mu\omega_\mathrm{e}}\frac{\mathrm{d}^2}{\mathrm{d}x^2} + \frac{\mu\omega_\mathrm{e}}{\hbar} x^2\right]\psi(x) = 2\frac{E}{\hbar\omega_\mathrm{e}} \cdot \psi(x)\]

Next, we introduce two abbreviations:

(2.78)\[\epsilon = \frac{E}{\hbar\omega_\mathrm{e}} \qquad\qquad z = \frac{x}{\sqrt{\hbar/(\mu\omega_\mathrm{e})}}\]

\(\epsilon\) is a dimensionless scaled energy, since both the numerator \(E\) and the denominator \(\hbar\omega_\mathrm{e}\) have dimension of energy. Similarly, \(z\) is a dimensionless scaled displacement, since both \(x\) and \(\sqrt{\hbar/(\mu\omega_\mathrm{e})}\) has dimension of length. Substituting these into the differential equation and rearranging gives a much more compact equation

(2.79)\[\frac{\mathrm{d}^2\psi(z)}{\mathrm{d}z^2} + (2\epsilon - z^2)\psi(z) = 0\]

Next, we determine the solution of this equation in the limit of very large \(|z|\). In that limit, \(2\epsilon\) is negligible relative to \(z^2\), and we can replace it by any other small number. Let’s take \(1\), since that will result in a differential equation with known solution. We get the differential equation

(2.80)\[\frac{\mathrm{d}^2\psi}{\mathrm{d}z^2} + (1-z^2)\psi = 0\]

for which the general mathematical solution is known:

(2.81)\[\psi(z) = c_1\cdot\mathrm{e}^{-z^2/2} + c_2\cdot \mathrm{e}^{-z^2/2}\int_0^z \mathrm{e}^{t^2}\mathrm{d}t\]

with arbitrary constants \(c_1\) and \(c_2\). The first term is a Gaussian function and converges rapidly towards zero as \(z\rightarrow\pm\infty\). In contrast, the second term diverges exponentially due to the integral - it explodes as \(z\) grows in magnitude. This renders the overall function not normalizable. Therefore, for a physically valid wavefunction, the second term has to vanish, i.e. \(c_2 = 0\). The resulting function is

(2.82)\[\psi(z) = c_1\cdot\mathrm{e}^{-z^2/2}\]

Knowing the asymptotic behavior of \(\psi\) for large \(|z|\), we can factor it out and set up a more general function, valid also for small \(z\)

(2.83)\[\psi(z) = H(z)\cdot\mathrm{e}^{-z^2/2}\]

with a function \(H(z)\) yet to be determined. (Don’t confuse this function with the Hamiltonian!) For the differential equation, we need the second derivative

(2.84)\[\psi'' = \left[ H'\cdot\mathrm{e}^{-z^2/2} -zH\cdot\mathrm{e}^{-z^2/2} \right]' = \left[ H''-2zH'+(z^2-1)H \right] \cdot \mathrm{e}^{-z^2/2}\]

(For notational compactness, we use \('\) to indicate the derivative with respect to \(z\).) Substituting this into the differential equation gives

(2.85)\[\left[ H''-2zH'+(2\epsilon-1)H \right] \cdot \mathrm{e}^{-z^2/2} = 0\]

The second factor on the left-hand side, the Gaussian function, is positive for any \(z\), so we can divide it out and get

(2.86)\[H''-2zH'+(2\epsilon-1)H = 0\]

This is now a differential equation for \(H(z)\). To solve it, we use the power series method. First, expand \(H(z)\) into an infinite power series of the form

(2.87)\[H(z) = \sum_{q=0}^{\infty}a_q z^q \qquad q = 0, 1, 2, ...\]

where \(a_q\) are the coefficients to be determined. Inserting this into the equation, and after a bit of algebra work, we get the equation

(2.88)\[\sum_{q=0}^\infty \big[ (q+1)(q+2)a_{q+2}+(2\epsilon-(2q+1))a_q \big] z^q =0\]

Since the right-hand side is zero for any value of \(z\), the left-hand side must be zero as well, for any value of \(z\). This is only possible if each of the terms in square brackets is zero:

(2.89)\[(q+1)(q+2)a_{q+2}+(2\epsilon-(2q+1))a_q = 0\]

Rearranging this, we get an expression for \(a_{q+2}\) in terms of \(a_q\), \(q\), and \(\epsilon\):

(2.90)\[a_{q+2} = a_q \cdot\frac{(2q+1)-2\epsilon}{(q+1)(q+2)}\]

This is a recursion relation for the polynomial coefficients: given \(a_q\) for any \(q\), we can calculate \(a_{q+2}\). From \(a_0\), we can calculate \(a_2\), \(a_4\), etc. From \(a_1\), we get \(a_3\), \(a_5\), etc. Therefore, besides \(\epsilon\), the only two free parameters are \(a_0\) and \(a_1\), all other follow from the recursion relation.

Now, we know all the coefficients of the polynomial \(H(z)\), and therefore we have the general mathematical solution of the time-independent Schrödinger equation. To determine the conditions under which this solution is physical, we need to figure out whether it is normalizable. For large \(q\), we can neglect the small numbers and \(2\epsilon\), and we get

(2.91)\[a_{q+2} \approx a_q\cdot\frac{2}{q}\]

This polynomial is as fast as the exponential \(\mathrm{e}^{z^2} = \sum_{q=0}^\infty z^{2q}/q!\), since for the latter the coefficients grow as \(q!/(q+1)!\approx1/q\) for large \(q\). But in order for the overall \(\psi\) to be square integrable, \(H(z)\) has to grow slower than this. Only then will the overall \(\psi\) converge to zero as \(z\rightarrow\pm\infty\).

For this, the polynomial must terminate, i.e. there must be a maximum \(q\) (let’s designate it \(n\)) above which all coefficients are zero: \(a_{n+1}=a_{n+2}=\dots=0\). That is, the following requirements must be fulfilled:

(2.92)\[\begin{split}\begin{eqnarray} a_{n+2} &=& a_n\cdot\frac{(2n+1)-2\epsilon}{(n+1)(n+2)} = 0 \\ a_{n+1} &=& a_{n-1}\cdot\frac{(2n-1)-2\epsilon}{n(n+1)} = 0 \end{eqnarray}\end{split}\]

for some value of \(n\). Starting with the equation for \(a_{n+2}\), it can be zero only if the numerator of the fraction is zero. This yields an equation for \(\epsilon\):

(2.93)\[\epsilon = n+\frac{1}{2} \qquad n = 0, 1, 2 ...\]

Now recall that \(\epsilon\) is the dimensionless scaled energy \(\epsilon = E/(\hbar\omega_\mathrm{e})\). From this, we see that the solutions are physical (square integrable) only if the energy satisfies

(2.94)\[E_n = \hbar\omega_\mathrm{e}\left(n+\frac{1}{2}\right) \qquad \qquad n = 0,1,2,...\]

\(n\) is called the vibrational quantum number. It is sometimes denoted with the Greek letter uspilson, \(\upsilon\).

With the result for \(\epsilon\), the equation for \(a_{n+1}\) becomes

(2.95)\[a_{n+1} = a_{n-1}\cdot\frac{(2n-1)-(2n+1)}{n(n+1)} = a_{n-1}\cdot\frac{-2}{n(n+1)}=0\]

This equality can only hold if \(a_{n-1}\) is zero, and consequently \(a_{n-3}\) etc. all the way down to \(n=0\). We see therefore that, in order to give physically valid wavefunctions, \(\epsilon=n+1/2\) and either \(a_0\) or \(a_1\) needs to be zero.

The resulting polynomials \(H_n(z)\) (with added subscript \(n\) to indicate the maximum exponent) are called Hermite polynomials. Each of them contains a finite number of terms with powers of \(n\), \(n-2\), etc. down to 0 or 1. The few Hermite polynomials for lowest \(n\) are:

(2.96)\[H_0(z) = 1 \qquad H_1(z) = 2z \qquad H_2(z) = 4z^2-2 \qquad H_3(z) = 8z^3-12z\]

For higher values of \(n\), they can be calculated from these using the recursion relation \(H_{n+1}(z)=2zH_n(z)-2nH_{n-1}(z)\).

With the Hermite polynomials, the overall wavefunctions are

(2.97)\[\psi_n(x) = H_n(z) \cdot \mathrm{e}^{-z^2/2}\]

As the last step, we need to normalize these wavefunctions. With \(\int_{-\infty}^\infty H_n^2(x)\mathrm{e}^{-x^2}\mathrm{d}x = 2^n n! \sqrt{\pi}\), and inserting the definition of \(z\), we get

(2.98)\[\psi_n(x) = \frac{1}{\sqrt{2^n n!}} \left(\frac{\mu \omega_\mathrm{e}}{\pi\hbar}\right)^{1/4} H_n\left(\sqrt{\frac{\mu\omega_\mathrm{e}}{\hbar}}x\right)\cdot \mathrm{e}^{-\mu\omega_\mathrm{e}x^2/2\hbar}\]