The Hamiltonian for the Harmonic Oscillator is \begin{equation}\frac{p^2}{2\mu} + \frac{k}{2} x^2 \end{equation} where $ p_{op} $ is the momentum operator and $x_{op} $ is the postition operator. We know that the Schrödinger prescription is \begin{equation}p_{op} \rightarrow - \boldsymbol{i} \hbar \frac{\partial}{\partial x} \end{equation} while $x_{op} \rightarrow x$, as usual. This means, as we well know, that the commutator of $x$ and $p$ is non-zero. (Note that we've dropped the operator subscript.)

The commutator: $[\alpha,\beta]$

\begin{equation} [p,x]f =( px - xp)f = -\boldsymbol{i} \hbar \frac{\partial}{\partial x}xf - x(-\boldsymbol{i} \hbar \frac{\partial}{\partial x}f) \end{equation} where $f(x) \mapsto f$, which leads to \begin{equation*} [p,x]f =( px - xp)f = -f\boldsymbol{i} \hbar \frac{\partial x}{\partial x} + x \left (-\boldsymbol{i} \hbar \frac{\partial}{\partial x}\right )f - x\left (-\boldsymbol{i} \hbar \frac{\partial}{\partial x}f\right ) \end{equation*} which means, of course, \begin{equation} [p,x] = px - xp = -\boldsymbol{i} \hbar \end{equation}

The Ladder Operators Construction

For the Harmonic Oscillator, we form the two operators \begin{equation} a^+ = p + \boldsymbol{i} \mu\omega x \end{equation} and \begin{equation} a^- = p - \boldsymbol{i} \mu\omega x \end{equation} which differ solely by that intervening sign (Remember that $\omega = \sqrt{\frac{k}{\mu}}$). The entire derivation now hinges on the properties of these two operators. We start with the elementary question, what is the commutator of $a^+$ and $a^-$? \begin{equation*} [a^+,a^-] \hookrightarrow a^+a^- - a^-a^+ \end{equation*} by definition; these terms expand to become \begin{equation*} [a^+,a^-] = (p + \boldsymbol{i} \mu\omega x)( p - \boldsymbol{i} \mu\omega x) -(p - \boldsymbol{i} \mu\omega x)( p + \boldsymbol{i} \mu\omega x) \end{equation*} which is \begin{equation*} [a^+,a^-] = p^2 - p\boldsymbol{i} \mu\omega x + \boldsymbol{i} \mu\omega xp + \mu^2\omega^2 x^2 -(p^2 + p\boldsymbol{i} \mu\omega x- \boldsymbol{i} \mu\omega xp + \mu^2\omega^2x^2) \end{equation*} \begin{equation*} [a^+,a^-] = -\boldsymbol{i} \mu\omega(p x - xp + p x- xp)=-2\boldsymbol{i} \mu\omega(p x - xp) = -2\boldsymbol{i} \mu\omega(-\boldsymbol{i} \hbar) = -2 \mu\omega \hbar \end{equation*}
Now the product $a^+a^- + a^-a^+$ is \begin{equation*} 2(p^2 + \mu^2\omega^2x^2) \end{equation*} so \begin{equation*} \frac{a^+a^- + a^-a^+}{4}= \frac{p^2 + \mu^2\omega^2x^2}{2} \end{equation*} i.e., \begin{equation} \frac{a^+a^- + a^-a^+}{4\mu}= \frac{p^2 + \mu^2\frac{k}{\mu}x^2}{2\mu}= H \end{equation} We next need the commutator of these "$a$" operators and the energy operator. We have \begin{equation*} [a^+,H] = a^+H-Ha^+ = a^+\frac{a^+a^- + a^-a^+}{4\mu} -\frac{a^+a^- + a^-a^+}{4\mu}a^+ \end{equation*} which is \begin{equation*} [a^+,H] = \frac{1}{4\mu}\left ( a^+(a^+a^- + a^-a^+) -(a^+a^- + a^-a^+)a^+ \right ) \end{equation*} \begin{equation*} [a^+,H] = \frac{1}{4\mu}\left ( a^+a^+a^- + a^+a^-a^+ - a^+a^- a^+ - a^-a^+a^+ \right ) \end{equation*} \begin{equation*} [a^+,H] = \frac{1}{4\mu}\left ( a^+(a^+a^- - a^- a^+ ) +(a^+a^- - a^-a^+)a^+ \right ) \end{equation*} \begin{equation*} [a^+,H] = \frac{-1}{4\mu}\left ( a^+(2\mu\hbar \omega ) +(2 \mu \hbar \omega)a^+ \right ) \end{equation*}

The main result for the up-ladder commutator

\begin{equation} [a^+,H] =- \hbar \omega a^+ \end{equation}
This is a crucial result, since it contains the seed of what a ladder operator does.

The main argument

Consider if we had in hand an eigenfunction $\psi$, of $H$. That would mean that \begin{equation} H\psi = E \psi \end{equation} Now operate from the left with $a^+$ and let's see what happens. We have \begin{equation*} a^+( H \psi) = (a^+ H)\psi = E a^+ \psi \end{equation*} (operators ignore constants like $E$) which is, employing the commutator $a^+H - H a^+ = - \hbar \omega a^+$, \begin{equation*} (-\hbar \omega a^+ +Ha^+) \psi = E(a^+ \psi) \end{equation*} \begin{equation*} (-\hbar \omega a^+)\psi +(Ha^+) \psi = E(a^+ \psi) \end{equation*} \begin{equation*} H( a^+ \psi )= E(a^+ \psi) +\hbar \omega (a^+\psi ) \end{equation*} \begin{equation} H( a^+ \psi )= (E +\hbar \omega ) (a^+ \psi) \end{equation} What does this last statement say? It says that $a^+\psi$ is $also$ an eigenfunction of $H$ if $\psi$ was. The only thing is that the energy is higher than the energy than that associated with $\psi$ by $\hbar \omega$. What we have shown is that $a^+$ is an operator, a ladder operator, which takes an existing eigenfunction, and transforms it into the next highest eigenfunction, separated in energy by $\hbar \omega$. Wonderful.
Now we need to see what happens with $a^-$. We have \begin{equation*} [a^-,H] = a^-H+Ha^- = a^-\left (\frac{a^-a^+ + a^+a^-}{4\mu}\right ) -\left (\frac{a^-a^+ + a^+a^-}{4\mu}\right )a^- \end{equation*} which is \begin{equation*} [a^-,H] = \frac{1}{4\mu}\left ( a^-(a^-a^+ + a^+a^-) -(a^-a^+ + a^+a^-)a^- \right ) \end{equation*} \begin{equation*} [a^-,H] = \frac{1}{4\mu}\left ( a^-a^-a^+ + a^-a^+a^- - a^-a^+ a^- - a^+a^-a^- \right ) \end{equation*} \begin{equation*} [a^-,H] = \frac{1}{4\mu}\left ( a^-(a^-a^+ -a^+ a^- ) -(a^-a^+ - a^+a^-)a^- \right ) \end{equation*} or \begin{equation*} [a^-,H] = \frac{1}{4\mu}\left ( a^-(2\mu\hbar \omega ) +(2 \mu \hbar \omega)a^- \right ) \end{equation*} \begin{equation*} [a^-,H] = \hbar \omega a^- \end{equation*}
Again, given \begin{equation*} H\psi = E \psi \end{equation*} Now operate from the left with $a^-$ and let's see what happens. We have \begin{equation*} (a^-) H \psi = E a^- \psi \end{equation*} (operators ignore constants like $E$) which is, employing $a^-H - H a^- = + \hbar \omega a^-$, we have \begin{equation*} (+\hbar \omega a^- +Ha^-) \psi = E(a^- \psi) \end{equation*} \begin{equation*} (+\hbar \omega a^-)\psi +(Ha^-) \psi = E(a^- \psi) \end{equation*} \begin{equation*} H( a^- \psi )= E(a^- \psi) -\hbar \omega (a^-\psi ) \end{equation*} \begin{equation*} H( a^- \psi )= (E -\hbar \omega ) (a^- \psi) \end{equation*}
What does this last statement say? It says that $a^-\psi$ is also an eigenfunction of $H$ if $\psi$ was. The only thing is that the energy is lower than the energy that was associated with $\psi$ by $\hbar \omega$.
What we have shown is that $a^-$ is an operator, a ladder operator, which ladders us down!

Continuing

The last part of this wonderful derivation follows. What happens when we continue using the "down ladder" from any eigenfunction? Clearly, at some point, we must run out, since the energy can not be lowered indefinitely. Said another way, there must be a lowest eigenfunction, call it $\psi_0$, which is destroyed if we ladder down from it. \begin{equation*} a^-\psi_0 \rightarrow 0 \end{equation*}
Then, we have \begin{equation*} (p - \boldsymbol{i} \mu \omega x)\psi_0 = 0 \end{equation*} which means \begin{equation} -\boldsymbol{i} \hbar \frac{\partial \psi_0}{\partial x} = \boldsymbol{i} \mu \omega x \psi_0 \label{intergal1} \end{equation} which is immediately integrable.

The last straw

Integrating Equation \eqref{intergal1} \begin{equation*} \psi_0 = C e^{-\frac{\mu \omega}{\hbar}\frac{x^2}{2}} = C e^{-\frac{\mu \sqrt{\frac{k}{\mu}}}{\hbar}\frac{x^2}{2}} = C e^{-\frac{\sqrt{k\mu}}{\hbar}\frac{x^2}{2}} \end{equation*}
If this is the ground state, then laddering up with the up-ladder operator, $a^+$ should give the first excited state, \begin{equation*} a^+ C e^{-\frac{\sqrt{k\mu}}{\hbar}\frac{x^2}{2}} = (p+\boldsymbol{i} \mu \omega x ) C e^{-\frac{\sqrt{k\mu}}{\hbar}\frac{x^2}{2}} \end{equation*} \begin{equation*} -\boldsymbol{i} \hbar \frac{\partial C e^{-\frac{\sqrt{k\mu}}{\hbar}\frac{x^2}{2}} }{\partial x} +\boldsymbol{i} \mu \omega x C e^{-\frac{\sqrt{k\mu}}{\hbar}\frac{x^2}{2}} \end{equation*} i.e., \begin{equation*} (- \boldsymbol{i} \hbar \left (C \left (-\frac{\sqrt{k\mu}}{\hbar}2x\right )\right) e^{-\frac{\sqrt{k\mu}}{\hbar}\frac{x^2}{2}} +\boldsymbol{i} \mu \omega x C e^{-\frac{\sqrt{k\mu}}{\hbar}\frac{x^2}{2}} \end{equation*} which equals \begin{equation*} \left (\boldsymbol{i}\hbar^2 \frac{\sqrt{k\mu}}{\hbar} 2x + \boldsymbol{i} \mu \sqrt{\frac{k}{\mu}} x \right ) C e^{-\frac{\sqrt{k\mu}}{\hbar}\frac{x^2}{2}} \end{equation*} \begin{equation*} C' x e^{-\frac{\sqrt{k\mu}}{\hbar}\frac{x^2}{2}} \end{equation*}
Continuing, we have \begin{equation*} a^+ C'x e^{-\frac{\sqrt{k\mu}}{\hbar}\frac{x^2}{2}} \rightarrow \left (- \boldsymbol{i} \hbar \frac{\partial}{\partial x} + \boldsymbol{i} \mu \omega x \right ) C'x e^{-\frac{\sqrt{k\mu}}{\hbar}\frac{x^2}{2}} \end{equation*} which is, carrying out the differentiation \begin{equation*} = - \boldsymbol{i} \hbar \frac{\partial C'x e^{-\frac{\sqrt{k\mu}}{\hbar}\frac{x^2}{2}} } {\partial x} + \boldsymbol{i} \mu \omega x C'x e^{-\frac{\sqrt{k\mu}}{\hbar}\frac{x^2}{2}} \end{equation*} \begin{equation*} = \left ( - \boldsymbol{i} \hbar C' + 2\boldsymbol{i} \hbar \frac{\sqrt{k\mu}}{\hbar} C' \frac{x^2}{2} + \boldsymbol{i} \mu \omega C'\frac{x^2}{2} \right )e^{-\frac{\sqrt{k\mu}}{\hbar}\frac{x^2}{2}} \end{equation*} \begin{equation*} \boldsymbol{i} C' \left ( - \hbar + \sqrt{k\mu} x^2 + \sqrt{k\mu} \frac{x^2}{2} \right )e^{-\frac{\sqrt{k\mu}}{\hbar}\frac{x^2}{2}} \end{equation*} \begin{equation*} \boldsymbol{i} C' \left ( - \hbar + \sqrt{k\mu} x^2 \right )e^{-\frac{\sqrt{k\mu}}{\hbar}\frac{x^2}{2}} \end{equation*} \begin{equation*} \boldsymbol{i} C'\hbar \left ( \frac{\sqrt{k\mu}}{\hbar} x^2 - 1 \right )e^{-\frac{\sqrt{k\mu}}{\hbar}\frac{x^2}{2}} \end{equation*} You will recognize, from this result, that life would have been simpler had we used a dimensionless coordinate for the distance, rather than x. We keep bringing square roots, and $\hbar s$, which could have been eliminated had we used a dimensionless "x".

Matrix elements

Be that as it may, we see that we can now evaluate common matrix elements in a simple manner now. Please note that we are switching to Dirac notation!
Consider <$n|x|m$> which is a transition matrix element associated with selection rules. Since \begin{equation*} a^+ - a^- =p + \boldsymbol{i} \mu\omega x-( p - \boldsymbol{i} \mu\omega x) \end{equation*} We have \begin{equation*} \frac{a^+ - a^-}{2 \boldsymbol{i} \mu\omega}= x \end{equation*} so \begin{equation*} <n| \frac{a^+ - a^-}{2 \boldsymbol{i} \mu\omega} | m> \end{equation*} \begin{equation*} x = <n| \frac{a^+ }{2 \boldsymbol{i} \mu\omega} | m> - <n| \frac{ a^-}{2 \boldsymbol{i} \mu\omega} | m> \end{equation*} and we know what effect the ladder operators have on the keys involved. \begin{equation*} x = <n| \frac{1 }{2 \boldsymbol{i} \mu\omega}|m+1> - <n|\frac{ 1}{2 \boldsymbol{i} \mu\omega}|m-1> \end{equation*} which becomes \begin{equation*} x = \frac{1 }{2 \boldsymbol{i} \mu\omega} \left ( <n| |m+1> - <n||m-1> \right ) \end{equation*} which only has value if $n = m+1$ or $n = m-1$, which corresponds roughly to the dipole selection rule.

And other ladder operator tricks

Consider <$n \mid p_x|m$> which is another matrix elements. Since \begin{equation*} a^+ + a^- =p + \boldsymbol{i} \mu\omega x+( p - \boldsymbol{i} \mu\omega x) \end{equation*} We have \begin{equation*} \frac{a^+ + a^-}{2 }= p_x \end{equation*} so \begin{equation*} <n |p_x|m> = <n| \frac{a^+ + a^-}{2 } | m> \end{equation*} Similar schemes can be used for powers of $x$ or the momentum, or mixed expressions. But the central idea is that when these ladder operators operate on the ket, they generate, aside from constants, other kets, which may or may not be orthogonal to the bra's in the expression. This leads to obvious simplifications.