Suppose we have a single particle in 1D, with wavefunction $\psi(x,t)$, obeying the Schrödinger equation in position space: $$i\hbar\frac{d\psi(x,t)}{dt}=\left(-\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2}+V(x)\right)\psi(x,t)$$ Suppose that the eigenenergies $E_n$ are nondegenerate, so $E_n\neq E_m$ for $n\neq m$.
Suppose that the probability density is constant in time everywhere: $$|\psi(x,t)|^2=\rho(x)$$ where $\rho(x)$ is not time-dependent. Is it necessarily true that $\psi(x,t)$ is an eigenfunction of the Hamiltonian, and how would you prove this? I have decomposed the density in some suggestive ways where it certainly looks like the density fluctuations would only cancel out for a stationary state, but I haven't actually been able to prove it.
I would certainly like to prove it in general, but I'm mostly interested in the case of the harmonic oscillator, so $V(x)=\frac{1}{2}m\omega^2x^2$. Can the result be proven in this specific case?
This property is definitely not true for a degenerate spectrum, but I don't care too much about that.
This is a follow-up to a previous question, where I was essentially asking the same question in the case of a many-body system. The property does not hold for many particles - one can come up with counterexamples, like if you have two uncorrelated particles that are exactly counter-oscillating, cancelling out the fluctuations. However, I'm almost certain that it's true in the case of a single particle in 1D with a nondegenerate energy spectrum.
Here's my thinking so far. We impose that $\frac{d}{dt}\left|\psi(x,t)\right|^2=0$. Expanding $\psi$ into the energy eigenbasis, we find that: $$ \sum_{nm}\frac{i}{\hbar}(E_n-E_m)\phi_n^*(x)\phi_m(x)\exp\left(\frac{it(E_n-E_m)}{\hbar}\right) = 0 $$ It seems to me like the only way this can be true is if we have an eigenstate, so every term is zero. But how can I prove that these terms don't all cancel out to give zero? Maybe using orthogonality of the basis or something?