The key concept to look for is displaced number states. These are, quite simply, the number states $|n⟩$, moved by the displacement operator
$$D(\alpha)=\exp\left[\alpha a^\dagger-\alpha^*a\right]$$
to some point $\alpha=x+ip$ on the complex phase space. The ground state of a harmonic oscillator which has been displaced to a real $\alpha=x$ is, as you note, the coherent state
$$|\alpha⟩=D(\alpha)|0⟩.$$
The rest of the eigenstates are, similarly,
$$|\alpha,n⟩=D(\alpha)|n⟩.$$
The simplest way to prove this is to use the commutation/displacement relations between the displacement operator and the ladder operators,
$$aD(\alpha)=D(\alpha)(a+\alpha)\quad\text{and}\quad a^\dagger D(\alpha)=D(\alpha)(a^\dagger+\alpha^*).$$
You can then transform the eigenvalue equation $a^\dagger a|n⟩=n|n⟩$ into
$$
\begin{align}
(a^\dagger-\alpha^*)(a-\alpha)|\alpha,n⟩
&=\left(D(\alpha)a^\dagger D(\alpha)^\dagger\right)
\left(D(\alpha)aD(\alpha)^\dagger\right) D(\alpha)|n⟩
\\ &= D(\alpha)a^\dagger a|n⟩=nD(\alpha)|n⟩
\\&=n~|\alpha,n⟩,
\end{align}
$$
or alternatively
$$
\left[a^\dagger a-(\alpha a^\dagger +\alpha^*a)\right]~|\alpha,n⟩
=(n-|\alpha|^*)~|\alpha, n⟩.
$$
Now, as all states in Hilbert space, the displaced number states can be written in the basis of number states, as $|\alpha,n⟩=\sum_{m=0}^\infty c_m(\alpha)|m⟩$. The coefficients in this expansion are simply the matrix elements of the displacement operator in the number basis:
$$⟨m|\alpha,n⟩=⟨m|D(\alpha)|n⟩.$$
This can be calculated in a couple of ways, which either look messy or seem pulled out of thin air, but what matters is ultimately the result. They come out as
$$
⟨m|D(\alpha)|n⟩=\sqrt{\frac{n!}{m!}}\alpha^{m-n}e^{-\tfrac12|\alpha|^2}L_n^{(m-n)}(|\alpha|^2)\quad\text{when }m\geq n,
\tag 1
$$
where $L_n^{(m-n)}$ is a Laguerre polynomial.
The standard reference in the literature for this matrix element is Appendix B of
Ordered Expansions in Boson Amplitude Operators. K. E. Cahill and R. J. Glauber. Phys. Rev. 177 no. 5, 1857-1881 (1969).
It's worth noting that most of the literature simply quotes the result $(1)$ and attributes it to Cahill and Glauber, but the most papers cite both the paper above and a second back-to-back paper, also by Cahill and Glauber, which does not contain the result. So be careful and show that you've read what you cite!
My bachelor's thesis,
E. Pisanty. Estados coherentes generalizados y estructura analítica del operador de aniquilación [pdf] (Generalized coherent states and the analytical structure of the annihilation operator). Tesis de licenciatura, Universidad Nacional Autónoma de México, 2011.
was on this topic and it contains an alternative calculation, though it's all in Spanish (un?)fortunately.
The discovery and study of coherent states represents one aspect of one
of the biggest problems physicists have faced with the birth and the
subsequent development, supported by excellent experimental results, of the
quantum mechanics: the search for a correspondence between the new theory, conceived
for the analysis of microscopic systems, and classical physics, still fully valid
for the description of the macroscopic world.
The history of coherent states begins immediately after the advent of mechanics
quantum: their introduction on a conceptual level dates back to a
article published in 1926, in which Schrödinger reports the existence of a class of
states of the harmonic oscillator that show, in a certain sense, behavior
analogous to that of a classic oscillator: for these states it is verified that the energy
mean corresponds to the classical value and the position and momentum averages have
oscillatory forms in constant phase relation.
Returning to Schrödinger's article, the "almost classical" states from him
identified present, in addition to the characteristics already mentioned, an important aspect:
being represented by Gaussian wave packets that do not change shape in the
time, guarantee the minimization of the product among the uncertainties about
position and on the impulse, that is the condition closest to the possibility of measuring
simultaneously the aforesaid quantities with arbitrary precision, allowed
from classical physics.
So, starting from the following relations:
\begin{equation}
a^{\dagger}|n\rangle=\sqrt{n+1}|n+1\rangle \quad a|n\rangle=\sqrt{n}|n-1\rangle
\end{equation}
it is noted that, by virtue of the orthonormality of the states
stationary, the diagonal matrix elements of the position and momentum operators are
null in the representation of energy, which means that the expectation values of
position and momentum on any stationary state are zero instant by instant.
The stationary states just analyzed are characterized by distributions of
constant probabilities with respect to the position over time; the wait values of the position e
of the impulse are null at all times: this aspect is a fundamental one
difference with the states of the classic oscillator, for which, once the energy is defined
(as long as different from zero), the observables position and momentum evolve over time
according to sinusoidal functions and are always in phase quadrature with each other. Also, if yes
calculate the uncertainties on position and momentum for a steady state with n photons, yes
gets the uncertainty relation $\Delta x \Delta p=(n+1/2)\hbar$.
it is therefore possible to obtain the minimization of the product of the uncertainties on impulse and
position, which represents the maximum similarity with classical mechanics.
A state that is as similar as possible to the classical case must therefore
have the following characteristics:
1): The evolution over time of the position and momentum expectation values must
be of a simple periodic type, with a constant phase ratio between position and
impulse.
2): The wave functions must be as narrow as possible around the value
average of the position, so that the probability distribution with respect to the
position may tend, by varying appropriate parameters, to a delta function of
Dirac;
3): The product of the uncertainties on the position and on the impulse must be minimal.
So you can see this "classical" behavior, that admits these particular states, as an intrinsic property of QHO.
Furthermore I have say that: also the harmonic oscillator energy eigenstates actually behave like oscillators.
Maybe this answer is a bit too long but I hope it can help you.
Best Answer
+1 for non-dimensionalizing $\hbar=1$ and the rest superfluous parameters.
Define $K\equiv \frac{1 }{2}(p^2 + x^2)$. Then you see that the canonical transformation $x\mapsto x-F(t), \qquad p\mapsto p$ preserves the commutation relations; and hence shifts the spectrum of K, to that of H, $$ e^{-ipF} K e^{ipF}=H+ F^2/2, $$ so the spectrum of H is $n+(1 -F^2)/2$.
It is then evident, as you surmised (but for an exponent factor of 1/2), that the unnormalized "ground" state of H is a plain Gaussian, $$ \tfrac{1}{2} ( -\partial^2+x^2 -2Fx) e^{-(x-F)^2/2}= {(1-F^2)\over 2} e^{-(x-F)^2/2}, $$ etc, for the excited states, mutatis mutandis....
When you diagonalize to Dirac creation and annihilation operators, $\sqrt{2} a(t)\equiv x-F(t)+ip$, you hardly see substantial traces of the shift, beyond the ferocious dependence of all former "constants" on time, $H(t)=a^\dagger (t) a(t) + (1-F(t)^2)/2$.
You may then address your problem, the TDSE, $$ \Bigl (i\partial_t - n-(1-F(t)^2)/2\Bigr )~~ \psi_n(x,t)=0, $$ ... it is not trivial to solve.$^\natural$
$^\natural$For instance, it is evidently impossible to solve for facile Ansätze such as $e^{-(x-F(t))^2+ G(t)}$, etc. You might have to devise very special Fs, but perhaps I have missed your desideratum.