Supersymmetric Quantum Mechanics (SUSYQM) is a method to solve most solvable potentials commonly seen in quanutum mechanics with relative ease. The idea behind this is similar to solving the harmonic oscilator with ladder operators - to factor the Hamiltonian in a day that makes it easy to compute the energy levels and stationary states.
Let \(\hbar=2m=1\), the overall goal of SUSYQM is to be able to solve for \(E_n\) and \(\psi_n\) in the Schrödinger equation: \[H\psi=\left(-\frac{d^2}{dx^2}+V(x)\right)\psi=E\psi\]
Using a similar idea as the harmonic oscilator solution, let’s define the operators \(A^\pm\) as \[A^\pm=\mp\frac d{dx}+W(x)\] where \(W\) is known as the superpotential. Define the Hamiltonians \(H^\pm\) as \[H^\pm=A^\mp A^\pm=-\frac{d^2}{dx^2}+W^2\pm W’=-\frac{d^2}{dx^2}+V^\pm\]
Let’s suppose we are given a Hamiltonian \(H^-\) with a ground state \(\psi_0^-\), i.e. \(H^-\psi_0^-=0\). It turns out with only the ground state information, we can recover \(V^-\) and \(W\) with \[V’=\frac{ {\psi_0^-}''}{\psi_0^-}\quad W=\frac{ {\psi_0^-}’}{\psi_0^-}\]
Let \(\psi_n^\pm\) be the normalizable eigenfunction to \(H^\pm\) with eigenvalue \(E_n^\pm\), with \(E_n^\pm<E_{n+1}^\pm\). Since by definition \(H^\pm\) is a product of conjugates, it is a semi-positive definite operator, so we also have \(E_n^\pm\geq0\).
We have the following way to relate both set of solutions: \[H^+\left(A^-\psi_n^-\right)=A^-A^+A^-\psi_n^-A^-H^-\psi_n^-=E_n^-\left(A^-\psi_n^-\right)\] \[H^-\left(A^+\psi_n^+\right)=A^+A^-A^+\psi_n^+A^+H^+\psi_n^+=E_n^+\left(A^+\psi_n^+\right)\] Hence as long as \(A^\pm\psi_n^\pm\) is renormalizable, it is a eigenfunction of \(H^\mp\) with eigenvalue \(E_n^\pm\).
The only case when \(A^\pm\psi_n^\pm\) is not renormalizable is when it vanishes. By construction we have \(A^-\psi_0^-=0\), this yields \(\psi_0^-\propto e^{-\int Wdx}\). If we have \(A^+\psi_0^+=0\), this implies that \(\psi_0+\propto e^{\int Wdx}\propto\frac1{\psi_0^-}\), which is not renormalizable, hence this is not possible and the only nonrenormalizable case is from \(\psi_0^-=0\). Hence we obtain \[E_{n+1}^-=E_n^+\]
Let’s consider the simple case where \(V^-\) is a infinite square well given by \[V^-(x)=\begin{cases}-1&0<x<\pi\\\infty&\text{otherwise}\end{cases}\] We have \(\psi_n^-\propto\sin\left((n+1)x\right)\) and \(E_n^-=(n+1)^2-1=n(n+2)\). Using the result above, we obtain \(W=-\cot x\) and \(V^+=2\csc^2-1\).
Using the isospectral relations above, the eigenvalues of \(H^+\) is given by \(E_n^+=n(n+2)\) and eigenfunctions are \[\psi_n^+=\left(\frac d{dx}-\cot(x)\right)\sin\left((n+1)x\right)=(n+1)\cos\left((n+1)x\right)-\cot(x)\sin\left((n+1)x\right)\]
Let’s introduce a parameter, \(a_n\) into our Hamiltonians, \(H^-\left(x,a_n\right)\) and \(H^+\left(x,a_n\right)\) where the ground state of \(H^-\) has energy \(0\) and let’s suppose that there exists some function \(g\) such that
\[H^+\left(x,a_k\right)+g\left(a_k\right)=H^-\left(x,a_{k+1}\right)+g\left(a_{k+1}\right)\]
If such a function exists, the potentials are shape invariant. We can also see this implies another relationship between the energies, now given by \(E_n^-\left(a_k\right)\) and \(E_n^+\left(a_k\right)\):
\[E_n^+\left(a_k\right)+g\left(a_k\right)=E_n^-\left(a_{k+1}\right)+g\left(a_{k+1}\right)\]
and this also implies that the wave functions \(\psi_n^-\left(x,a_{k+1}\right)\) and \(\psi_n^+\left(x,a_k\right)\) are equal.
Furthermore, this lets us compute \(E_n^-\left(a_k\right)\) with just the function \(g\). As an example, we can compute \(E_2^-\left(a_1\right)\) with
\[
\begin{align*}
E_2^-\left(a_1\right)&=E_1^+\left(a_1\right)\\\
&=E_1^-\left(a_2\right)+g\left(a_2\right)-g\left(a_1\right)\\\
&=E_0^+\left(a_2\right)+g\left(a_2\right)-g\left(a_1\right)\\\
&=E_0^-\left(a_3\right)+g\left(a_3\right)-g\left(a_2\right)+g\left(a_2\right)-g\left(a_1\right)\\\
&=g\left(a_3\right)-g\left(a_1\right)\\\
\end{align*}
\]
More generally, \(E_n^-\left(a_k\right)=g\left(a_{n+k}\right)-g\left(a_k\right)\), which is significantly easier than solving the eigenvalue problem.
The solution to the hydrogen atom can be simplified down to solving the Schrödinger equation for \(V=-\frac{ke^2}r+\frac{\ell(\ell+1)\hbar^2}{2mr}\), setting \(e=k=1\) as well, we have \(V=-\frac1r+\frac{\ell(\ell+1)}{r^2}\).
Consider the superpotential \(W(\ell,r)=\frac1{2(\ell+1)}-\frac{\ell+1}r\), with this, we obtain the following potentials.
\[V^+=\frac14\left(\frac1{\ell+1}\right)^2-\frac1r+\frac{(\ell+1)(\ell+2)}{r^2}\]
\[V^-=\frac14\left(\frac1{\ell+1}\right)^2\underbrace{-\frac1r+\frac{\ell(\ell+1)}{r^2}}_{\text{Coulomb potential}}\]
And we immediately notice the columb potential inside and also the shape invariance with \(a_n=\ell+n\) and \(g\left(a_n\right)=-\frac14\left(\frac1{a_n+1}\right)^2\). With this, we can easily find the energy levels of the hydrogen atom:
\[E_n(\ell)=g\left(a_{\ell+n}\right)-g\left(a_\ell\right)-\frac14\left(\frac1{\ell+1}\right)^2=-\frac14\left(\frac1{\ell+n+1}\right)^2+\frac14\left(\frac1{\ell+1}\right)^2-\frac14\left(\frac1{\ell+1}\right)^2=-\frac14\left(\frac1{\ell+n+1}\right)^2\]