It is well known that if one takes the Schrödinger operator $$H = - \frac{1}{2} \frac{d^2}{dx^2} + \frac{1}{2} x^2$$ acting on some dense subspace of $L^2(\mathbb{R})$ then this has discrete spectrum $E_n = n + \frac{1}{2}$. This operator is important in quantum field theory because one can interpret $n+\frac{1}{2}$ as the energy of $n$ particles, each of energy $1$, plus a zero point energy of $\frac{1}{2}$.
This leads to the following question. Given a Schrödinger operator with potential $V$ (assumptions on $V$ left open) $$H = - \frac{1}{2} \frac{d^2}{dx^2} + V(x)$$ having only discrete spectrum $E_n$. Is it then the case that if $E_n = an +b$ then $V$ must be a quadratic function?
The Hamilton operator $H_0$ of a one-dimensional harmonic oscillator, $$H_0:= P^2/2+a^2 X^2/2+(b-a/2) \mathbf{1}, \quad a>0, \, b \in \mathbb{R}, \tag{1} \label{eq1}$$ with $(X \psi)(x)=x \psi(x)$, $(P\psi)(x)=-i\psi^\prime(x)$ acting on an appropriately chosen dense domain $D(H_0) \subset L^2(\mathbb{R})$ has the pure point spectrum $\sigma(H_0)= \{an+b, \;n=0,1,2,\ldots\}$.
The task is now to find all Hamilton operators of the form $$H^\prime = P^2/2 +V(X) \tag{2} \label{eq2} $$ having the same spectrum as $H_0$ (i.e. $\sigma(H^\prime) =\sigma(H_0)$).
The existence of such operators with $H^\prime \ne H_0$ is obvious, as any shift $X^2 \to (X-c)^2$ (with $c \in \mathbb{R}$) of the harmonic oscillator potential leaves the spectrum unaffected. Defining $$H_c:= P^2/2+a^2 (X-c)^2/2+(b-a/2) \mathbf{1}, \quad c \in \mathbb{R}, \tag{3} \label{eq3} $$ the shifted operators $H_c$ can be related to the original one by the unitary translation operator $$T_c = \exp(i c P) \tag{4} \label{eq4}$$ via the relation $H_c= T_c^\ast H_0 T_c$.
Conversely, let us assume that an operator $H^\prime$ of the form given in \eqref{eq2} is isospectral to $H_0$ given in \eqref{eq1}. The normalized eigenfunctions $\{\phi_0, \phi_1, \ldots \}$ of $H_0$ (forming a complete orthonormal system of $L^2(\mathbb{R})$) and the corresponding eigenfunctions $\{\chi_0, \chi_1,\ldots\}$ of $H^\prime$ are related by the unitary transformation $U$ defined by $\phi_n=U\chi_n$ for $n=0,1,2,\ldots$, such that $$ H^\prime = U^\ast H_0 U, \quad U^\ast U = U U^\ast = \mathbf{1}. \tag{5} \label{eq5}$$ As, by the assumed form of $H^\prime$ shown in \eqref{eq2}, the kinetic term $P^2/2$ should remain unaffected by this transformation, we infer $U=U(P)$, ruling out e.g. the Galilei transformed Hamiltonian $H^{\prime \prime}=(P-V)^2/2 +a^2(X-Vt)^2/2+(b-a/2) \mathbf{1} $, being isospectral to $H_0$ but not of the form \eqref{eq2} because of the presence of the additional term $-VP$ linear in $P$.
We thus end up with the relation $$V(X)= a^2 U^\ast(P) X^2 U(P)/2 +(b-a/2)\mathbf{1}, \quad U^\ast(P) U(P) = U(P) U^\ast(P) = \mathbf{1}. \tag{6} \label{eq6} $$ At this point, it turns out to be slightly more convenient to switch to the momentum representation, where a wave function $\psi(x)$ in $x$-space is related to the corresponding wave function $\tilde{\psi}(p)$ in momentum space by the Fourier transformation $$\tilde{\psi}(p)= \int\limits_{\mathbb{R}}dx \frac{e^{-ipx}}{\sqrt{2 \pi}} \psi(x) =: (\mathcal{F}\psi)(p). \tag{7} \label{eq7} $$ The momentum operator $\tilde{P} = \mathcal{F} P \mathcal{F}^ {-1}$ is now given by the multiplication operator $(\tilde{P} \tilde{\psi})(p)= p \tilde{\psi}(p)$, whereas the position operator $\tilde{X}= \mathcal{F} X \mathcal{F}^ {-1}$ becomes a differential operator, $(\tilde{X} \tilde{\psi})(p)= i \tilde{\psi}^\prime(p).$ Eq. \eqref{eq6} assumes the equivalent form $$V(\tilde{X})= a^2 U^\ast(\tilde{P}) \tilde{X}^2 U(\tilde{P})/2 +(b-a/2) \mathbf{1}, \tag{8} \label{eq8}$$ where $U^\ast(\tilde{P}) \tilde{X}^2 U(\tilde{P})$ acting on a momentum-space wave function $\tilde{\psi}$ takes now the simple form $$- \overline{U(p)} \frac{d^2}{dp^2} (U(p)\tilde{\psi}(p)),$$ where $U(p) = \exp (i f(p))$ for some real function $f(p)$. As $\exp{(i f(p)} )\tilde{\psi}(p)$ should be in the domain of $\tilde{X}^2$, we require differentiabilty of $f$ and $f^\prime$, such that eq. \eqref{eq8} assumes the final form $$ V(\tilde{X}) = a^2[\tilde{X}^2 -2 f^\prime(\tilde{P}) \tilde{X}-i f^{\prime \prime}(\tilde{P})+f^\prime(\tilde{P})^2 ]/2+ (b-a/2) \mathbf{1}. $$ As the left-hand side is independent of $\tilde{P}$, we obtain the condition $f^\prime(p) =c$, leaving \eqref{eq3} as the only possibility for \eqref{eq2} with $\sigma(H^\prime) = \sigma(H_0)$.