Is the harmonic oscillator potential unique in having equally spaced discrete energy levels? I was wondering if the good old quadratic potential was the only potential with equally spaced eigenvalues. Obviously you can construct others, such as a potential that is infinite in some places and quadratic in others, but that's only trivially different. I am not referring to equally spaced as a limiting behavior either, I mean truly integer spaced.
Any ideas? If not, is there a proof for its uniqueness?
If there are other potentials with equally spaced eigenvalues, can one use them as starting points for a free-field QFT? It would be interesting to know if there is a deeper mathematical relation between all of these potentials and whether they could be used to study interacting systems.
 A: I believe my answer here to the question "How does one determine ladder operators systematically?" gives at least a partial answer to your question. It is a partial answer because I assume a little more than your bare question, but then, as we see by looking carefully at yuggib's answer observes you can obviously write down a hamiltonian with equally spaced eigenvalues and then characterise the whole family of such hamiltonians. It becomes clear that you need to talk about more than simply the Hamiltonian to answer your question: we need to define other observables and how they behave to define something akin to a "potential".
Let's look at yuggib's answer. You can obviously write down a hamiltonian with equally spaced eigenvalues. Then, as you counter, in your comment:

"Does this really qualify as a proof? Sure you can make operators with any set of discrete eigenvalues, but how do you know that they are not all just equivalent to each other and the harmonic potential? Also, I should say that I really mean a potential, because obviously you can write a Hamiltonian as a matrix and give it equally spaced eigenvalues, but how could you know what potential gave rise to it?"

For this question, I choose energy eigenstates bounded from below. Otherwise, you could get arbitrarily negative energy states and there would be no quantum ground state. This may or may not be more than what you want to assume, but I think it is physically reasonable. As I said, I'm giving a partial answer. So now our index set $I$ becomes in fact the set of semi positive integers $\mathbb{N}$. So, in the notation of yuggib's answer, choose an orthonormal basis $\left\{Y_j\right\}_{j=0}^\infty$ for our assumed separable (this is yet another assumption we must bring to bear) Hilbert space of quantum states with $P_j Y_k = \delta_{j\,k} Y_k$, where $\delta$ is clearly the Kronecker delta and $P_j$ are the projection operators onto the basis vectors. Then, the most general Hamiltonian with equispaced eigenvalues is as written in yuggib's answer with:
$$\lambda_i = E_0 + \sigma(i) \Delta$$
where $E_0$ is the ground state energy, $\Delta$ the energy spacing and $\sigma:\mathbb{N}\to\mathbb{N}$ a bijection between $\mathbb{N}$ and itself. So there are an infinite number of Hamiltonians with equispaced energy levels. All members of each family of such Hamiltonians defined by the family's ground state energy $E_0$ and spacing $\Delta$ are unitarily equivelent to one another: two members $\hat{H}_1,\,\hat{H}_2$ are equivalent by $\hat{H}_1 = U\,\hat{H}_2\,U^\dagger$ with $U$ some unitary operator.
So now, how to work this into something like a "potential"? My solution is then to abstractly define position and momentum observables $\hat{X},\,\hat{P}$ and we make our final three assumptions:


*

*They fulfill the canonical commutation relationship $\hat{X}\,\hat{P} - \hat{P}\,\hat{X}=i\,\hbar\,\mathrm{id}$;

*Measurements by these observables vary sinusoidally with time;

*Our observables are Hermitian operators.
So now we Writing a general quantum state as:
$$\psi = \sum\limits_{j\in\mathbb{N}} \psi_j e^{-i\,\left(j\,\omega_0+\frac{E_0}{\hbar}\right)\,t}$$
so that the mean of a general observable $\hat{A}$ is:
$$\left<\left.\psi\right|\right.\hat{A}\left.\left|\psi\right.\right> = \sum\limits_{j=0}^\infty a_{j,j}|\psi_j|^2  + 2\,{\rm Re}\left(\sum\limits_{j=0}^\infty\sum\limits_{k=j+1}^\infty a_{j,k} \psi_j \psi_k^* \exp(i\,\omega_0\,(k-j)\,t)\,\right)$$
we can readily see that observables with sinusoidally varying measurement means must have two symmetrically placed, complex conjugate off-leading-diagonal diagonal stripes. Moreover the displacement from the leading diagonal must be the same for both $\hat{X},\,\hat{P}$ if they are to fulfill the CCR. The simplest case is when the two stripes are immediately above and below the leading diagonal. In this case, the means of the observables will vary like $\cos(\omega_0\, t + \phi_0)+const.$: if the two stripes are displaced $N$ steps either side of the leading diagonal, then we have variation varies like $\cos(N\,\omega_0\, t + \phi_0)+const.$. The case with the stripes displaced $N$ from the leading diagonal yield analyses that are essentially the same as the below as discussed in my other answer: they essentially pertain to a quantum oscillator with $N$ times the energy spacing of the one we talk about here.
So now, without loss of generality, in our assumed basis we can write the Hermitian observables as:
$$\begin{array}{ll}\hat{X} = \sqrt{\frac{\hbar}{2}}\left(\tilde{X} + \tilde{X}^\dagger\right)  \\
\hat{P} = \sqrt{\frac{\hbar}{2}}\left(\tilde{P} + \tilde{P}^\dagger\right)\end{array}$$
where:
$$\tilde{X} = \left(\begin{array}{ccccc}0&x_1&0&0&\cdots\\0&0&x_2&0&\cdots\\0&0&0&x_3&\cdots\\\cdots&\cdots&\cdots&\cdots&\cdots\end{array}\right)\quad\tilde{P}=\left(\begin{array}{ccccc}0&p_1&0&0&\cdots\\0&0&p_2&0&\cdots\\0&0&0&p_3&\cdots\\\cdots&\cdots&\cdots&\cdots&\cdots\end{array}\right)$$
are both arbitrary lone-striped upper triangular matrices. You'll need to look at my other answer that I referenced above for full details, but by writing down the CCR we find that:
$$\begin{array}{ll}\hat{X} = \sqrt{\frac{\hbar}{2\,m\,\omega_0\,\cos\chi}}\left(a^\dagger\,e^{-i\,\xi} + a\,e^{i\,\xi}\right)  \\
\hat{P} =  i\,\sqrt{\frac{\hbar\,m\,\omega_0}{2\,\cos\chi}}\left(a^\dagger\,e^{-i\,(\xi+\chi)} - a\,e^{i\,(\xi+\chi)}\right)  \end{array}$$
where we have defined the arbitrary complex constant $\alpha = -i\,m\,\omega_0\,e^{i\,\chi}$ by writing it in terms of a second real positive magnitude $m$ with the dimensions of mass and an arbitrary phase factor $\chi$ and where we have also defined:
$$a = \left(\begin{array}{ccccc}0&\sqrt{1}&0&0&\cdots\\0&0&\sqrt{2}&0&\cdots\\0&0&0&\sqrt{3}&\cdots\\\cdots&\cdots&\cdots&\cdots&\cdots\end{array}\right)$$
and its Hermitian conjugate as the wonted ladder operators. 
One can quite straighforwardly show that our derived $\hat{X},\,\hat{P}$ must have continuous spectrums. So, if we change our co-ordinates so that we work in position co-ordinates, i.e. where $\hat{X}$ becomes the diagonal (multiplication) operator $\hat{X}f(x) = x f(x)$ for $f(x)\in \mathcal{H} = \mathbf{L}^2(\mathbb{R}^3)$, then we can argue as I do in my answer here that there is needfully such a co-ordinate system wherein $\hat{X}f(x) = x f(x)$ and $\hat{P} f(x) = -i\,\hbar\,\nabla f(x)$. So now the Hamiltonian in these co-ordinates is:
$$\begin{array}{lcl}\hat{H} &=& \hbar\,\omega_0 \left(a^\dagger\,a + \frac{1}{2}I\right) + \left(E_0 - \frac{\hbar\,\omega_0}{2}\right)I\\
 &=& \frac{1}{2\,m\,\cos\chi} \hat{P}^2 + \frac{1}{2\,\cos\chi}\,m\,\omega_0^2\,\hat{X}^2 - \frac{\omega_0\,\tan\chi}{2}(\hat{X}\hat{P} + \hat{P}\hat{X}) +\left(E_0 - \frac{\hbar\,\omega_0}{2}\right)I\end{array}$$
so that the Schrödinger equation in these co-ordinates is:
$$i\hbar\partial_t \psi = -\frac{1}{2\,m\,\cos\chi} \nabla^2 \psi + \frac{1}{2\,\cos\chi} m\,\omega_0^2 |\vec{x}|^2 \psi + i\,\hbar\,\omega_0\,\tan\chi\,\vec{x}\cdot\nabla \psi + \left(E_0 +i\,\frac{\tan\chi\,\hbar\,\omega_0}{2} - \frac{\hbar\,\omega_0}{2}\right)\psi$$
and more "traditional" Schrödinger equation is recovered when $\chi = 0$ and $E_0 = \hbar\,\omega_0/2$:
$$i\hbar\partial_t \psi = -\frac{1}{2\,m} \nabla^2 \psi + \frac{1}{2} m\,\omega_0^2 |\vec{x}|^2 \psi$$
and we see that we must have a quadratic potential. So, in summary, let's list the assumptions that lead to this conclusion:


*

*Separable quantum state space;

*Equispaced energy levels bounded from below;

*The existence of Hermitian position and momentum observables fulfilling the CCR and

*Measurements from the observables vary sinusoidally with time.


I have not yet analysed the case where we relax the assumption 4. From the above, we see that equispaced energy levels implies periodic variations of measurements with time, and it seems to me that this relaxation would likely yield a much more general potential in position co-ordinates.
A: Very interesting problem! I will propose it to my students.
However the answer is positive up to isomorphisms of Hilbert spaces if assuming that every eigenspace has dimension $1$.
I stress that the result below holds true even if the initial Hamiltonian $H$ is not of Schroedinger form.
Proposition. Let $H$ be an (unbounded) self-adjoint operator over the Hilbert space $\cal H$ such that for a real constant $\omega \neq 0$
$$\sigma(H)= \sigma_p(H) = \{\omega n \:|\: n=0,1,2,\ldots\}\tag{1}$$ and every eigenspace has dimension $1$.
Under these hypotheses there is a unitary map $U: {\cal H} \to L^2({\mathbb R}, dx)$ such that, for some uniquely fixed real numbers  $\alpha, \beta$
$$U H U^{-1} = \alpha H_0 + \beta I\tag{2}$$
where $H_0$ is the standard (self-adjoint) Hamiltonian of the harmonic oscillator.
Proof. Let $\psi_n$ be the unit eigenvector of the eigenvalue $\omega n$ defined up to a phase. The finite span $D$ of the $\psi_n$ is dense, since they define the spectral measure of $H$ from (1).
We can define the standard ladder operators over the common domain $D$
$$a \psi_n = \sqrt{n}\psi_{n-1}$$ (with $\psi_{-1}:=0$) and $$a^\dagger \psi_n = \sqrt{n+1}\psi_{n+1}\:.$$ Evidently,$D$ turns out to be invariant under $a$ and $a^\dagger$.
Next define, over the dense domain $D$, the symmetric (not yet selfadjoint) operators
$A = \frac{1}{\sqrt{2}}(a+a^\dagger)$ and $B=\frac{i}{\sqrt{2}}(a-a^\dagger)$.
We have $A^2+B^2 = (\omega H + \frac{1}{2}I)|_D$.
Here is the crucial observation, since the $\psi_n$ are a Hilbert basis of eigenvectors of $H$, $D$ is a core for $H$ and also for $A^2+B^2+ I$ since those vectors are analytic vectors for $H$ and thus also for trivially related symmetric operator $A^2+B^2+I^2$
By a known theorem by Nelson [1]  $A^2+B^2+I^2$  is essentially self-adjoint over $D$.
A famous theorem by Nelson (Nelson theorem about unitary representations [1]) assures that there exist a strongly-continuous unitary representation of the unique connected simply-connected Lie group whose Lie-algebra is generated by $A,B,I|_D$. It easy to see that this Lie algebra is nothing but that of Heisenberg group.
Finally, Nelson's theorem also establishes that the unique self-adjoint extensions of $A$, $B$ and $I|_D$ are the self-adjoint generators of the three fundamental one-parameter groups of Heisenberg group.
In other words, we have found that $exp(-iuA)$ and $exp(-ivB)$ have the same commutation relations of $exp(-iuX)$ and $exp(-ivP)$
Here Stone-von Neumann theorem enters the game (in the generalized version due to Mackay, without requiring the irreducibility of the representation).
There exist a unitary operator $U : {\cal H} \to \bigoplus_{k\in K}L^2(\mathbb{R}, dx)$ such that $$U(\omega H + \frac{1}{2}I)U^{-1} = \bigoplus_{k\in K} H_0$$
so that
$$H = \frac{1}{\omega} U^{-1}\left(\bigoplus_{k\in K} H_0- \frac{1}{2}I\right) U$$
It is clear that if $K$ contains more than one element, $\omega H + \frac{1}{2}I$ would have more than one eigenvector for a given eigenvalue contrarily to our hypotheses. So that $K$ includes only one element and we can write
$$H =  U^{-1} \left(\frac{1}{\omega}H_0- \frac{1}{2\omega}I\right) U\:.$$
QED
[1] Nelson, E.: Analytical vectors, Annals of Mathematics, 70, 572-615 (1969).
A: According to Inverse spectral theory as influenced by Barry Simon Fritz Gesztesy
(Submitted on 2 Feb 2010) page 4:

A particularly interesting open problem in inverse spectral theory
  concerns the characterization of the isospectral class of potentials $V$ with purely discrete spectra (e.g., the harmonic oscillator $V(x) = x^2$ ).

A: There are, in general, infinitely many operators with equally spaced eigenvalues. Suppose a self-adjoint operator $A$ has a purely discrete spectrum (i.e. it is either compact or with compact resolvent) and denote by $\{\lambda_i\}_{i\in\mathbb{I}}$ its real eigenvalues ($I\subseteq \mathbb{N}$): then by the spectral theorem it can be written (on its domain of definition) as
$$A=\sum_{i\in I} \lambda_i P_i$$ where $P_i$ is the orthogonal projector on the eigensubspace corresponding to $\lambda_i$.
Now this equality, read from right to left defines the operator $A$, choosing the eigenvalues and the (mutually disjoint) orthogonal projections. So playing with projections and eigenvalues you will define different operators, with equally spaced eigenvalues if you want. However it may be then necessary to prove the resulting operator is self-adjoint on a suitable domain (if the operator is unbounded).
A: It's been a while since I have done Quantum Mechanics, however your question reminded me of Supersymmetric Quantum Theory.  In simpler words it is a theorem that states that for every potential ($V(x)$) there exists a supersymmetic (SUSY) partner potential ($\tilde{V}(x)$) that has the same functional form for the energy levels and differ in in there energies by a shift.
Quantum Mechanics by Schwabl has a section devoted to this (19) and actually have the Harmonic Oscillator as an example (19.2.3).  However the answer was not satisfying because the SUSY potentials they got to be for it were...
$$V(x) = \frac{1}{2}\omega^2 x^2 -\frac{1}{2}\omega$$
$$\tilde{V}(x) = \frac{1}{2}\omega^2 x^2 +\frac{1}{2}\omega$$
Which are essentailly equivalent because adding constants to a Hamiltonian only shifts the energy.
So according to SUSY QM the Harmonic Oscillator can only be mapped to a Harmonic Oscillator with a constant energy shift.
Don't know if I completely answered your question, but I hope it helps!
A: The harmonic potential is not the only potential with evenly spaced energy levels.
Consider a potential $V(X)$ that is $\frac{1}{2} m \omega_1 x^2 - \frac{1}{2} \hbar \omega_1$ for $x > 0$ and $\frac{1}{2} m \omega_2 x^2 - \frac{1}{2}\hbar \omega_2$ for $x < 0$. If $\omega_2$ is a square-free multiple of $\omega_1$ then the energy levels will be spaced with levels $2 \hbar \omega_1$. The reason for this is that we can use the even $\omega_1$ energy eigenstates for the $x > 0$ part of the function and the corresponding $\omega_2$ energy eigenstates for the $x < 0$ part. At $x = 0$ we can set the functions to be equal, and the first derivative is $0$ there as well. We can not use the odd solutions because the first derivatives of the odd solutions will never be equal at $x = 0$ assuming that $\omega_2/ \omega_1$ is square-free. This is because the wave function is given by a Hermite polynomial $H_n(\sqrt{\frac{m \omega}{\hbar}} x)$ times the Gaussian curve, and for the odd solutions, the $x$ term in $H_n$ is always an integer.
In any case, let us now restrict our attention to even potentials $V(x)$, which. If $V(x)$ is even, then the harmonic potential is the only potential that has evenly spaced energy levels.
If the energy levels are evenly spaced with spacing $\Delta E$, then the time evolution operator
$$U(t) := e^{i \hat H t / \hbar}$$
will satisfy
$$U(t + \tfrac{2 \pi}{ \omega}) = e^{i \theta} U(t)$$
for all $t$, where $\omega = \Delta E/\hbar$ and $\theta$ is a physically irrelevant phase that can be set to $0$ by adding a constant to our Hamiltonian (such as setting the ground state energy to $0$).
This means that under time evolution, all state will return to their initial state every period $\frac{2 \pi}{\omega}$. Also note that by Erenfests theorem, the classical position expectation value will follow the classical equation of motion, so that too will return to its initial state after the same period.
I will now prove that the harmonic potential is the only symmetric 1D potential for which all solutions are periodic with the same period. Therefore, it is the only potential with evenly spaced energy levels.
First note that $V(x)$ must have only one minimum and no maximums. If there were a maximum, then there would exist "slow roll" solutions, rolling off of the hill, that can take arbitrarily long times to leave the hill, such that not all solutions would be periodic with the same period. Therefore, $V(x)$ must reach a minimum (WLOG at $x=0$) and increase as $|x|$ gets larger.
From the conservation of energy
$$\frac{1}{2} m v^2 = E - V(x)$$
$$v = \sqrt{\frac{2}{m}(E - V(x))}$$
$$v = \frac{dx}{dt}$$
$$\frac{dx}{\sqrt{\frac{2}{m}(E - V(x))}} = dt$$
$$\int \frac{dx}{\sqrt{\frac{2}{m}(E - V(x))}} = \int dt$$
We can consider the interval from the maximum value of $x$ that the particle reaches, $V^{-1}(E)$, to $0$.
$$\int_0^{V^{-1}(E)} \frac{dx}{\sqrt{\frac{2}{m}(E - V(x))}} = T$$
We need the right hand side to be independent of $E$ if all solutions are to have the same period. (From here on out we will suppress constants like $2$ and $m$, only wanting the integral to be independent of $E$.)
First we change variables to integrate with respect to $V$:
$$\int_0^{E} \frac{1}{V'}\frac{dV}{\sqrt{E - V}}$$
where $V' = \frac{dV}{dx}$.
Because $V$ only has one minimum, we can express $V'$ as some function $f(V)$.
$$\int_0^{E} \frac{1}{f(V)}\frac{dV}{\sqrt{E - V}}$$
In order for this integral to be independent of $E$, it must be dimensionless. The only way to accomplish this is for $f(V) \propto \sqrt{V}$. Indeed, the quantity
$$\int_0^{E} \frac{1}{\sqrt{V}}\frac{dV}{\sqrt{E - V}}$$
can also be explicitly evaluated to be independent of $E$. (It is $\pi$.)
This means that $V \propto x^2$, which is the harmonic potential. This is what we wanted to show.
