The hard way

We first follow the bruteforce method. We have a second order differential equation, and we bite the bullet and sit down to solve it. We can use the coordinate basis where \(\hat{x}\) and \(\hat{p}\) have the following representations: \[\begin{equation} \hat{x}= x, \text{ and }\, \hat{p}=-i\hbar \frac{d}{dx} \tag{7}. \end{equation}\]

Therefore, the coordianate part of the Schrödinger equation becomes: \[\begin{equation} \left(-\frac{\hbar^2}{2m}\frac{d^2}{dx^2} +\frac{1}{2} m \omega^2 x^2\right)\psi(x) = E\psi(x). \tag{8} \end{equation}\] We can define a couple of dimensionless quantities \(\chi=x\sqrt{\frac{m\omega}{\hbar}}\) and \(\epsilon=\frac{2E}{\hbar\omega}\) to get: \[\begin{equation} \left(-\frac{d^2}{d\chi^2} +\chi^2- \epsilon \right)\psi(x)=0. \tag{9} \end{equation}\] We should first try to understand the asymptotic solution where \(\chi^2\gg \epsilon\): \[\begin{equation} \left(-\frac{d^2}{d\chi^2} +\chi^2 \right)\psi(x)\simeq 0. \tag{10} \end{equation}\] This equation has a special solution called parabolic cylinder functions. However, since we are looking for the asymptotic solutions, we can make an educated guess of the form \(e^{-\alpha \chi^2}\) and plug it in to find that \(\alpha=1/2\) solves the differential equation at the first order. Or, we can try to split the \(-\frac{d^2}{d\chi^2} +\chi^2\) operator into two first order operators and drop a small term in the large \(\chi\) limit: \[\begin{equation} \left(-\frac{d^2}{d\chi^2} +\chi^2 \right)\psi(x)\simeq \left(-\frac{d}{d\chi} +\chi \right)\left(\frac{d}{d\chi} +\chi \right)\psi(x) \simeq 0. \tag{11} \end{equation}\]

We then combine the asymptotic solution with a yet-unknown function and propose a solution of the following form: \[\begin{equation} \psi(\chi)=e^{-\frac{\chi^2}{2}}h(\chi) \tag{12}, \end{equation}\] upto the normalization constant, which we will calculate later. Plugging this back into Eq. (9), we get the following second order differential equation. \[\begin{equation} \left(-\frac{d^2}{d\chi^2} +2\chi \frac{d}{d\chi} +1- \epsilon \right)h(\chi)=0. \tag{13} \end{equation}\] We can now try a power series of the form \[\begin{equation} h(\chi)=\sum_{j=0}^\infty c_j \chi^j \tag{14}. \end{equation}\] Inserting this back in, we get: \[\begin{equation} -\sum_{j=0}^\infty j(j-1) c_j \chi^{j-2}+\sum_{j=0}^\infty (2 j +1-\epsilon)c_j\chi^j=0 \tag{15}. \end{equation}\] Since the first two terms in the first summation vanish due to the \(j(j-1)\) coefficient, we can start the first summation index, \(j\), from \(2\), and redefine \(j\) as \(j+2\) and pull the starting point back to \(0\):

\[\begin{equation} -\sum_{j=0}^\infty (j+2)(j+1) c_{j+2} \chi^{j}+\sum_{j=0}^\infty (2 j +1-\epsilon)c_j\chi^j=\sum_{j=0}^\infty \left(-(j+2)(j+1) c_{j+2}+(2 j +1-\epsilon)c_j\right)\chi^j=0 \tag{16}. \end{equation}\] In order to set this to zero we need to have the recurrence equation: \[\begin{equation} c_{j+2}=\frac{2j +1-\epsilon}{(j+2)(j+1)} c_j \tag{17}. \end{equation}\] Note that this is problematic because the coefficients are not decaying fast enough. In fact, this relation implies that \(h(\chi)\propto e^{\chi^2}\), and even the prefactor \(e^{-\chi^2/2}\) will not be decaying fast enough to make the wave-function normalizable. The only way out of this is to truncate the series at some point. Remember that the only knob we have is \(\epsilon\), and we can set it an integer value such that when \(2j+1=\epsilon\), the series terminates. This is a profound finding because the physicality of the solution requires the quantization of the energy. Going back to the original parameters, \(E=\frac{\epsilon \hbar\omega}{2}\), we can write the energy as: \[\begin{equation} E=\hbar\omega(n+\frac{1}{2}) \tag{18}. \end{equation}\] There is another subtle problem: note that the recurrence formula relates \(c_0\) to \(c_2\), \(c_2\) to \(c_4\) so and so forth, and \(c_1\) to \(c_3\), \(c_3\) to \(c_5\). In other words, the only free coefficients are \(c_0\) and \(c_1\). As we discussed earlier, we can truncate the series by selecting \(\epsilon\) appropriately. However, we have only one degree of freedom in \(\epsilon\), and we can’t use that to terminate both of the odd and even series at the same time. That means only one of the coefficients \(c_0\) and \(c_1\) can be non-zero at the same time. This is also expected from the parity symmetry of the Hamiltonian: it stays invariant under \(x\to-x\), which implies that solutions should stay invariant up to the sign. Therefore, odd and even powers of \(\chi\) cannot mix in the energy eigenstates.

Since the series in Eq. (14) will terminate at \(j=n/2\) due to the recurrence relation in Eq. (17), it is convenient to redefine the summation index \(s=\frac{n}{2}-j\), and rewrite the solution as a sum of finite number of terms:

\[\begin{equation} h_n(\chi)=\sum_{s=0}^{\lfloor n/2\rfloor}(-1)^s\frac{n!}{(n-2s)!s!} (2x)^{n-2s} \tag{19}, \end{equation}\] where \((-1)^s\) originates from flipping the sign of the numerator in Eq. (17), and powers of \(2\) originate from \(n/2\)’s in the denominator. These are Hermite polynomial, which can be written explicitly as \[\begin{eqnarray} h_0(\chi)&=&1\nonumber\\ h_1(\chi)&=&2x\nonumber\\ h_2(\chi)&=&4x^2-2\nonumber\\ h_3(\chi)&=&8x^3-12x\nonumber\\ h_4(\chi)&=&16x^4-48x^2+12\\ &\vdots&\tag{20} \end{eqnarray}\] Now we have to deal with the normalization of the wavefunction in Eq. (12). There is a very elegant way of doing this via the generating function. Let’s multiply Eq. (19) with \(\frac{t^n}{n!}\) and sum over \(n\):

\[\begin{eqnarray} g(\chi,t)&\equiv& \sum_{n=0}^\infty \frac{t^n}{n!} h_n(\chi)= \sum_{n=0}^\infty \sum_{s=0}^{\lfloor n/2\rfloor}(-1)^s\frac{t^n}{(n-2s)!s!} (2\chi)^{n-2s} =\sum_{n=0}^\infty \sum_{s=0}^{\infty}(-1)^s\frac{t^n}{(n-2s)!s!} (2\chi)^{n-2s}\tag{21}, \end{eqnarray}\] where we extended the summation upper limit since we will negative factorials give negative infinities killing all the terms for \(s>n/2\). Now define \(n-2s=m\) and do some shuffling: \[\begin{eqnarray} g(\chi,t)&=&\sum_{m=0}^\infty \sum_{s=0}^{\infty}(-1)^s\frac{t^{m+2s}}{m!s!} (2\chi)^{m} =\sum_{s=0}^{\infty}\frac{(-t^2)^s}{s!}\sum_{m=0}^\infty \frac{ (2t\chi)^{m}}{m!}=e^{-t^2+2t\chi}\tag{22}. \end{eqnarray}\] Now it becomes an easy task to compute the normalization factor. Consider the following: \[\begin{eqnarray} \int_{-\infty}^\infty d\chi e^{-\chi^2}g(\chi,t) g(\chi,q) &=&\int_{-\infty}^\infty d\chi e^{-\chi^2-t^2+2t\chi-q^2+2q\chi} =\int_{-\infty}^\infty d\chi e^{-(\chi -t-q)^2+2 q t}=\sqrt{\pi} e^{2 q t}=\sqrt{\pi}\sum_{n=0}^\infty \frac{(2qt)^n}{n!} \tag{23}, \end{eqnarray}\] and evaluate the integral in the series expansion: \[\begin{eqnarray} \int_{-\infty}^\infty d\chi e^{-\chi^2}g(\chi,t) g(\chi,q) &=&\int_{-\infty}^\infty d\chi e^{-\chi^2}\sum_{n=0}^\infty \frac{t^n}{n!} h_n(\chi) \sum_{m=0}^\infty \frac{q^m}{m!} h_m(\chi) =\sum_{n=0}^\infty \sum_{m=0}^\infty \frac{t^n}{n!}\frac{q^m}{m!} \int_{-\infty}^\infty d\chi e^{-\chi^2} h_n(\chi)h_m(\chi)\nonumber\\ &=&\sum_{n=0}^\infty \frac{(qt)^n}{(n!)^2} \int_{-\infty}^\infty d\chi e^{-\chi^2} h_n(\chi)h_n(\chi) +\sum_{n=0}^\infty \sum_{m=0, m\neq n}^\infty \frac{t^n}{n!}\frac{q^m}{m!} \int_{-\infty}^\infty d\chi e^{-\chi^2} h_n(\chi)h_m(\chi) \tag{24}. \end{eqnarray}\] By comparing the coefficients of \(q t\) terms in Eqs. (23) and (24), we first see that the cross terms should vanish. We also get the normalization constant:

\[\begin{eqnarray} \int_{-\infty}^\infty d\chi e^{-\chi^2} h_n(\chi)h_m(\chi)=2^n\sqrt{\pi} n!\delta_{n,m} \tag{25}. \end{eqnarray}\] The orhogonality of the eigenfunctions is not a coincidence since the differential operator we are dealing with can be transformed into a self-adjoint form, and the orthogonality is guaranteed due to the Sturm-Liouville theory[1]. Putting it all together we have the normalized energy eigenstates: \[\begin{equation} \psi_n(x)= \left(\frac{m\omega}{\pi\hbar}\right)^{\frac{1}{4}} \frac{1}{\sqrt{2^n n!}}h_n\left(\sqrt{\frac{m\omega}{\hbar}} x\right) e^{-\frac{m\omega x^2}{2\hbar}}, \quad n=0,1,2,...\quad. \tag{26} \end{equation}\]

It is operationally more practical to combine \(\hat x\) and \(\hat p\) operators into raising and lowering ladder operators. The harder method is based on the recurrence relations of the Hermite polynomials. Taking the derivative of the equality in Eq. (21) with respect to \(\chi\), we get: \[\begin{eqnarray} \frac{\partial}{\partial \chi}g(\chi,t)&=&2t e^{-t^2+2t\chi}=2 \sum_{n=0}^\infty \frac{t^{n+1}}{n!} h_n(\chi)=2 \sum_{m=1}^\infty \frac{t^{m}}{(m-1)!} h_{m-1}(\chi)=2 \sum_{m=1}^\infty \frac{m t^{m}}{m!} h_{m-1}(\chi)=2 \sum_{n=0}^\infty \frac{n t^{n}}{n!} h_{n-1}(\chi) \nonumber\\ &=&\sum_{n=0}^\infty \frac{t^n}{n!} h'_n(\chi)\tag{27}, \end{eqnarray}\] where we first defined \(m=n+1\), and then relabeled \(m\) as \(n\). We also added the vanishing \(n=0\) term in the summation to make the sum start from \(0\). Matching the coefficients of \(t^n\) terms, we get the first recurrence relation of the Hermite functions: \[\begin{eqnarray} 2n h_{n-1}(\chi)= h'_n(\chi)\tag{28}. \end{eqnarray}\]

Let’s try taking the derivative with respect to \(t\) to get: \[\begin{eqnarray} \frac{\partial}{\partial t}g(\chi,t)&=&(-2 t+ 2\chi) e^{-t^2+2t\chi}=-2 \sum_{n=0}^\infty \frac{t^{n+1}}{n!} h_n(\chi)+2 \sum_{n=0}^\infty \frac{t^{n}}{n!} \chi h_n(\chi)=\sum_{n=0}^\infty \frac{t^{n}}{n!}\left\{ 2\chi h_{n}(\chi) - 2n h_{n-1}(\chi) \right\} \nonumber\\ &=&\sum_{n=0}^\infty n\frac{t^{n-1}}{n!} h_n(\chi)=\sum_{n=0}^\infty \frac{t^{n}}{n!} h_{n+1}(\chi)\tag{29}. \end{eqnarray}\] Matching the coefficients of \(t^n\) terms, we get the second recurrence relation of the Hermite functions: \[\begin{eqnarray} h_{n+1}(\chi)= 2\chi h_{n}(\chi) - 2n h_{n-1}(\chi) \tag{30}. \end{eqnarray}\]

We can combine Eqs. (28) and (30) to get another flavor:

\[\begin{eqnarray} h_{n+1}(\chi)= \left(2\chi - \frac{d}{d\chi}\right) h_{n}(\chi) \tag{31}. \end{eqnarray}\]

Now consider the following operator acting on \(\psi_n(\chi)\) as it is defined Eq. (26): \[\begin{eqnarray} \frac{1}{\sqrt{2}}\left(\chi -\frac{d}{d\chi}\right)\psi_n(\chi)&=& \frac{1}{\sqrt{2}} \left(\frac{m\omega}{ \pi\hbar}\right)^{\frac{1}{4}} \frac{1}{\sqrt{2^n n!}} \left(\chi -\frac{d}{d\chi}\right)\left(h_n(\chi) e^{-\frac{\chi^2}{2}}\right) \nonumber\\ &=& \left(\frac{m\omega}{ \pi\hbar}\right)^{\frac{1}{4}} \frac{\sqrt{n+1}}{\sqrt{2^{n+1} (n+1)!}}e^{-\frac{\chi^2}{2}} \left(2\chi -\frac{d}{d\chi}\right)h_n(\chi)=\sqrt{n+1}\psi_{n+1}(\chi) \tag{32}, \end{eqnarray}\] where we used Eqs. (31).

Now consider another operator acting on \(\psi_n(\chi)\): \[\begin{eqnarray} \frac{1}{\sqrt{2}}\left(\chi +\frac{d}{d\chi}\right)\psi_n(\chi)&=& \frac{1}{\sqrt{2}} \left(\frac{m\omega}{ \pi\hbar}\right)^{\frac{1}{4}} \frac{1}{\sqrt{2^n n!}} \left(\chi +\frac{d}{d\chi}\right)\left(h_n(\chi) e^{-\frac{\chi^2}{2}}\right) = \left(\frac{m\omega}{ \pi\hbar}\right)^{\frac{1}{4}} \frac{1}{\sqrt{2^{n+1} n!}}e^{-\frac{\chi^2}{2}} h'_n(\chi) \nonumber\\ &=& \left(\frac{m\omega}{ \pi\hbar}\right)^{\frac{1}{4}} \frac{1}{2\sqrt{n} \sqrt{2^{n-1} (n-1)!}}e^{-\frac{\chi^2}{2}} 2 n h_{n-1}(\chi)=\sqrt{n}\psi_{n-1}(\chi) \tag{33}, \end{eqnarray}\] where we used Eqs. (28). The operators \(\frac{1}{\sqrt{2}}\left(\chi \pm \frac{d}{d\chi}\right)\) can be written interms of \(x\) and \(\hat p\) and they will be called \(a\) and \(a^\dagger\), and that would be how one solves the harmonic oscillator the hard way. Now let’s look into the method of operators.

The operational way

We can define the ladder operators as follows: \[\begin{equation} a=\frac{1}{\sqrt{2m\hbar \omega}} (m\omega x +ip), \quad a^\dagger =\frac{1}{\sqrt{2m\hbar \omega}} (m\omega x -ip) \iff \hat{x}=\sqrt{\frac{\hbar}{2m\omega}}(a+a^\dagger), \quad \hat{p}= -i\sqrt{\frac{m\hbar\omega}{2}}(a-a^\dagger) \tag{34}. \end{equation}\] The commutation relation \([x,p]=i\hbar\) turns to \[\begin{equation} [a,a^\dagger]=1. \tag{35} \end{equation}\] The Hamiltonian becomes \[\begin{equation} H \equiv \hbar\omega(a^\dagger a +\frac{1}{2}) \tag{36}. \end{equation}\] Comparing Eq. (36) with Eq. (18) we can associate \(a^\dagger a\) with number operator: \[\begin{equation} N= a^\dagger a \tag{37}, \end{equation}\] which returns the state number: \[\begin{equation} N|n \rangle = n|n \rangle. \tag{38} \end{equation}\] Let’s now figure out how \(a\) and \(a^\dagger\) act on eigenstate \(|n\rangle\). We can read the energy value by acting on the new state with \(H\): \[\begin{equation} Ha|n \rangle = \left( aH+[H,a]\right)|n \rangle=\left( aH-a\hbar\omega\right)|n \rangle =\hbar\omega\left(n-1+\frac{1}{2}\right)a|n \rangle \tag{39}, \end{equation}\] which shows that the state \(a|n \rangle\) can be indexed as \(n-1\), i.e., \(a|n \rangle =c |n-1\rangle\) where \(c\) is the normalization constant. The overall coefficient \(c\) can be calculated as \[\begin{equation} |a|n \rangle|^2 =\langle n | a^\dagger a|n \rangle= n \langle n|n \rangle=n =|c|^2 \implies c=\sqrt{n} \tag{40}. \end{equation}\] Therefore the lowering operator \(a\) does the following: \[\begin{equation} a|n \rangle =\sqrt{n}|n-1\rangle \tag{41}. \end{equation}\] As a consequence of Eq. (41), we see that the ground state, \(|0\rangle\), will be annihilated by the operator \(a\) \[\begin{equation} a|0\rangle =0 \tag{42}. \end{equation}\] Let’s repeat for \(a^\dagger\): \[\begin{equation} Ha^\dagger|n \rangle = \left( a^\dagger H+[H,a^\dagger]\right)|n \rangle=\left( a^\dagger H+a^\dagger\hbar\omega\right)|n \rangle =\hbar\omega\left(n+1+\frac{1}{2}\right)a|n \rangle \tag{43}, \end{equation}\] which shows that the state \(a^\dagger|n \rangle\) can be indexed as \(n+1\), i.e., \(a^\dagger|n \rangle =d |n+1\rangle\) where \(d\) is the normalization constant. The overall coefficient \(d\) can be calculated as \[\begin{equation} |a^\dagger|n \rangle|^2 =\langle n | a a^\dagger|n \rangle=\langle n | a^\dagger a +[a,a^\dagger] |n \rangle= (n+1) \langle n|n \rangle=n =|d|^2 \implies d=\sqrt{n+1} \tag{44}. \end{equation}\] Therefore the lowering operator \(a^\dagger\) does the following: \[\begin{equation} a^\dagger|n \rangle =\sqrt{n+1}|n+1\rangle \tag{45}. \end{equation}\]

By recursively applying \(a^\dagger\) on \(|0 \rangle\) we can get the \(n\)-th energy eigenstate, \(|n\rangle\): \[\begin{equation} |n \rangle = \frac{(a^\dagger)^n}{\sqrt{n!}} |0 \rangle \tag{46}. \end{equation}\]

Let’s take a look at certain expectation values. We can immediately see that the expected values of \(x\) and \(\hat p=i\hbar\frac{d}{dx}\) vanish since the integrands of \(\langle \psi_n|x|\psi_n\rangle\) and \(\langle \psi_n|\hat p|\psi_n\rangle\) are odd and the integration range is symmetric around the origin. Equivalently we can do the computation using the ladder operators: \[\begin{equation} \langle \hat{x}\rangle_n= \langle n|\hat{x}|n \rangle =\sqrt{\frac{\hbar}{2m\omega}} \langle n|(a+a^\dagger)|n \rangle =0 \tag{47}. \end{equation}\] Similarly for \(p\), we have: \[\begin{equation} \langle \hat{p} \rangle_n= \langle n|\hat{p}|n \rangle =-i\sqrt{\frac{m\hbar\omega}{2}} \langle n|a-a^\dagger|n \rangle =0 \tag{48}. \end{equation}\]

Now consider the quadratic operators: \[\begin{equation} \langle \hat{x}^2\rangle_n= \langle n|\hat{x}^2 |n \rangle =\frac{\hbar}{2m\omega} \langle n|(a+a^\dagger)^2|n \rangle=\frac{\hbar}{2m\omega} \langle n|(2a^\dagger a +[a,a^\dagger])|n \rangle =\frac{\hbar}{2m\omega} (2n+1) \tag{49}. \end{equation}\] Similarly for \(p\), we have: \[\begin{equation} \langle \hat{p}^2 \rangle_n= \langle n|\hat{p}^2 |n \rangle =-\frac{m\hbar\omega}{2} \langle n|(a-a^\dagger)^2|n \rangle =\frac{m\hbar\omega}{2} \langle n|2a^\dagger a +[a,a^\dagger]|n \rangle =\frac{\hbar}{2m\omega} (2n+1) \tag{50}. \end{equation}\] The uncertainity in \(x\) and \(p\) for state \(n\) are given as:

\[\begin{equation} \langle (\Delta x)^2 \rangle_n = \langle \hat x^2 \rangle_n - (\langle \hat x \rangle_n)^2 =\frac{\hbar}{2m\omega} (2n+1) \tag{51}, \end{equation}\] and \[\begin{equation} \langle (\Delta p)^2\rangle _n = \langle \hat p^2 \rangle_n - (\langle \hat p \rangle_n)^2 =\frac{\hbar m\omega}{2} (2n+1) \tag{52}. \end{equation}\]

The Heisenberg relation becomes

\[\begin{equation} (\Delta x)^2 (\Delta p)^2 = \frac{\hbar^2}{4}(2n+1)^2, \tag{53} \end{equation}\] which has the minimum value of \(\frac{\hbar^2}{4}\) at \(n=0\).