Does a non-stationary state for a single particle with a nondegenerate spectrum necessarily have fluctuations in probability density? 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?
 A: Your expansion is not quite the most general and should more generally read
$$
\Psi(x,t)=\sum_m c_m e^{-i E_m t/\hbar} \psi_m(x)
$$
so that 
\begin{align}
\frac{d\rho}{dt}=0\quad\Rightarrow 0= \sum_{mn} c_m c_n^* (E_m-E_n) e^{-i(E_m-E_n)t/\hbar}
\psi_n(x)^*\psi_m(x)\tag{1}
\end{align}
(up to a factor of $i\hbar$).   For simplicity order the eigenvalues so that $E_0\le E_1\ldots $.
(1) must be valid for any $t$ which means that, for all those $m’,n’$ so that $\Delta E_{m’,n‘}=E_{m’}-E_{n‘}$ is the same you must have
$$
\Delta E_{m’,n’}\sum_{m’,n’} c_{m’}c_{n’}^* \psi_{m’}(x)\psi_{n’}(x)^*=0\, . 
\tag{2}
$$
Note that, if $c_{m’}c_{n’}^*\psi_{m’}(x)\psi_{n’}(x)^*$ appears in the sum (2), then its complex conjugate will not appear as then $\Delta E_{m’,n’}=-\Delta E_{n’,m’}$. In other words, each of the sums in (2) is of one of the
following types:


*

*The first is over $m’>n’$ with 
$\Delta E_{m’,n’}> 0$,

*the next is over $m’<n’$ with $\Delta E_{m’,n’}<0$, 

*the last is for $\Delta E_{m’,n’}=0$, which implies $m’=n’$ since the energies are assumed non-degenerate.


For the first case we have
$$
\Delta E_{m’n’}\sum_{m’>n’} c_{m’}c_{n’}^*\psi_{n’}(x)^*\psi_{m’}(x) = 0\, ,
\qquad 
\Delta E_{m’n’}\ne 0\, . \tag{3}
$$
 However, the set $\{\psi_{n’}(x)^*\psi_{m’}(x), m’>n’\}$ is a linearly independent basis set so there is no $c_{m’}c_{n’}^*$ to satisfy (2).  Hence there is no solution unless you admit $\Delta E_{m’n’}=0$, a contradiction.
The second case is done in the same way.  Only the third case leads to no contradiction, which shows the result.
A: This is likely not what you are asking for, but the answer to your question is yes for the oscillator, a highly symmetric and transparent quantum system: a time-invariant probability for a pure state entails a stationary state. 
This is manifest in the phase space formulation of QM, that treats the oscillator far more transparently and less abusively than the equivalent Hilbert space formulation. Unfortunately, however, it is hard for me now to provide an acceptable translation to Hilbert space as you evidently seek, and you may question whether it is worth learning a new language to ascertain a substantially narrow point  (some cognoscenti dismiss the oscillator as "basically semiclassical"). In any case, I provide this narrow truth for future readers;  it is all but impossible to say anything about general potentials beyond the oscillator. Formal details are accessible in our booklet Concise Treatise of QM in Phase Space where the two figures come from.
The probability density in coordinate space you are discussing is but the momentum integral of the Wigner function, a quasiprobability density in phase space,
$$
\rho (x;t)= \int dp  f(x,p;t) \\
{f(x,p;t)=\frac{1}{2\pi}\int\! dy~\psi^* \left (x-\frac{\hbar}{2} 
y ;t\right )~e^{-iyp} ~ \psi \left (x+\frac{\hbar}{2} y ;t\right )}.
$$
In the 30s, Wigner  showed that, for the oscillator, 
$$f(x,p;t)=f(x \cos t - p \sin t,  ~   p \cos t + x \sin t  ; ~0),$$
by solving the suitable dynamical equations for his function.
That is, the time evolution of an arbitrary pure state is rigid rotation in phase space, just as in classical mechanics. He thus all but trivialized the quantum description of the oscillator, substantially more recondite in Hilbert space. Rigid rotation, not radial pulsing, or any other motion. 

Fig 6 clockwise rigidly rotating oscillator pure state Wigner function for generic states (shoeboxes) and coherent states (yellow circles).
Now the momentum (marginal) integral of the Wigner functions schematically plotted on the x-axis is $\rho(x;t)$, oscillating back and forth on it for the (yellow) coherent states, for example. 
Evidently, then, for this $\rho(x)$ projection to be time-invariant, your state's Wigner function must be phase-space axially symmetric.
The only axially symmetric pure state Wigner functions are those resulting out of single-energy wavefunctions, that is stationary states, like the 3rd excited state Wigner function below:

Fig 3 Stationary 3rd excited state of the oscillator.
There are long and winding arguments for why the only axially symmetric Wigner functions are those corresponding to stationary states, the eigenstates of the oscillator hamiltonian, but it is a fact.
So, in this language, by virtue of the marvelously circumscribed universal motion of all oscillator states, you are done: time-invariance of the density projection excludes sneaky motions which wash out under integration over momenta.
If only one could generalize this to broader hamiltonians...
