How can one recover the classical frequency-modulation Bessel sidebands from a quantum emitter in a harmonic well? Consider an atom that emits at frequency $\omega_0$ that's located at a position $x$ that moves under the influence of a harmonic oscillator at frequency $\nu$. In both classical and quantum physics, the Doppler shift induced by the motion on the emitted radiation causes the appearance of sidebands on either side of the carrier-frequency emission at frequency $\omega_0$, and I would like to better understand the relationship between those two descriptions.

1. Classical description
In the classical side of things, you have a well-defined atomic position $x(t) = x_0\cos(\nu t)$, and you emit the signal $s(t)=s_0\cos(\omega_0t)$ as seen by the atom, but when observed from a fixed observation point $x=x_\mathrm{f}$, the path length between the emitter and the observation point will change in time and therefore so will the phase of the observed signal, which will go as
\begin{align}
s(t)
& = s_0 \cos(\omega_0t-k(x_\mathrm f-x(t)))
\\& = s_0 \cos(\omega_0t+kx_0\cos(\nu t)+\varphi),
\end{align}
where $\varphi=kx_\mathrm f$ is a fixed phase that one can set to zero, and $k=2\pi/\lambda$ is the wavevector of the emitted radiation.
Now, because $s(t)$ is no longer purely harmonic, its spectrum is no longer purely monochromatic, but luckily we can expand that pesky cosine as a harmonic series, 
\begin{align}
s(t)
& = s_0 \cos(\omega_0t+kx_0\cos(\nu t))
\\ & = 
\mathrm{Re}\left[
s_0 e^{-i\omega_0t}e^{-ikx_0\cos(\nu t)}
\right]
\\ & = 
\mathrm{Re}\left[
s_0 e^{-i\omega_0t}\sum_{n=-\infty}^\infty (-i)^nJ_n(kx_0)e^{-in\nu t}
\right]
\\ & = 
\mathrm{Re}\left[
\sum_{n=-\infty}^\infty s_0(-i)^nJ_n(kx_0) e^{-i(\omega_0+n\nu) t}
\right],
\end{align}
via the standard expansion of exponentiated cosines as Bessel functions. This gives, therefore, a spectrum with a strong carrier at frequency $\omega_0$, and then a bunch of sidebands separated by multiples of the atomic-position oscillation frequency $\nu$, with the strength of the sidebands given by Bessel functions $J_n(kx_0)$ evaluated at the amplitude of the oscillation.

2. Quantum description
Consider now the exact same system, but with both the emitter and its position quantized. In that spirit, then, we have a two-level system $\{|g\rangle,|e\rangle\}$ with atomic excitation frequency $\omega_0$ (and therefore atomic hamiltonian $\hat H_\mathrm{atom} = \frac12 \omega_0\hat{\sigma}_z$ as usual), interacting via some dipole coupling $\hat d = d_{eg}(|e\rangle\langle g| + |g\rangle\langle e|)=d_{eg}\hat{\sigma}_x$ with a classical laser mode at frequency $\omega_L$, so that the interaction hamiltonian reads
\begin{align}
\hat{H}_\mathrm{int}
& =
\vec E \cdot \hat{\vec{d}}
= E_0\hat{d} \cos(\omega_L t-k\hat x),
\end{align}
where now $\hat x=x_\mathrm g(\hat a+\hat a^\dagger)$, the position of the atom, is quantized and under the action of the harmonic hamiltonian $\hat H_\mathrm{motion} = \nu \, \hat{a}^\dagger \hat a$ with ground-state characteristic width $x_\mathrm{g}$. This gives us, then, a total hamiltonian of the form
\begin{align}
\hat H
& = \hat H_\mathrm{atom} + \hat{H}_\mathrm{int} + \hat H_\mathrm{motion} 
\\ & = \tfrac{1}{2}\omega_0\hat{\sigma}_z + \nu \, \hat{a}^\dagger \hat a + \Omega \hat{\sigma}_x\cos(\omega_L t-k\hat x)
\\ & = 
\tfrac{1}{2}\omega_0\hat{\sigma}_z 
+ \nu \, \hat{a}^\dagger \hat a 
+ \tfrac12 \Omega (\hat{\sigma}_++\hat{\sigma}_-)\left(
e^{ikx_\mathrm{g}(\hat a+\hat a^\dagger)}e^{-i\omega_L t}
+e^{-ikx_\mathrm{g}(\hat a+\hat a^\dagger)}e^{+i\omega_L t}
\right)
\end{align}
where $\Omega = E_0d_{eg}$ is the usual Rabi frequency, and one often defines $\eta=kx_\mathrm{g}$, the Lamb-Dicke parameter. In general, this hamiltonian is best understood in the rotating-wave approximation and in an interaction picture that takes into account the solvable dynamics of the atom and its motion when taken separately, which transforms it to the form
\begin{align}
\hat H_\mathrm{eff}
& = 
\tfrac12 \Omega\left(
e^{i\eta(\hat a+\hat a^\dagger)}\hat{\sigma}_+e^{-i\Delta t}
+e^{-i\eta(\hat a+\hat a^\dagger)}\hat{\sigma}_-e^{+i\Delta t}
\right),
\end{align}
where $\Delta=\omega_L-\omega_0$ is the detuning and the bosonic operators have been transformed with a time-dependent phase.
If we didn't have the quantized position, then the atomic excitation operators $\hat{\sigma}_\pm$ would mediate transitions between the atomic eigenstates, at the atomic frequency $\omega_0$, but now we have a quantized mode to exchange energy with, and the interaction can therefore operate internal transitions in the atomic motion via the coupling strengths
$$
\Omega_{mn} = \Omega \langle m|e^{i\eta(\hat a+\hat a^\dagger)}|n\rangle
$$
between states with well-defined translational-motion energy. These transitions form sidebands on either side of the carrier, which see everyday use as the primary objects in sideband cooling of ions in ion traps, and they occur, as with the classical FM sidebands, separated from the atomic frequency $\omega_0$ by multiples of the trap frequency $\nu$. 
(The quantum sidebands also have some differences with the classical sidebands, since in the Lamb-Dicke limit of $\eta\ll 1$ the red sidebands become weaker than blue sidebands, an effect that's absent from the classical description, but I'm not that interested in those features at the moment.)

3. My question
How can one recover classical-like expressions for the sidebands, and in particular the Bessel-function measures for their strengths as a function of the Lamb-Dicke parameter and the motional oscillation amplitude, as a suitable limit of the quantum Lamb-Dicke dynamics? 
I'm primarily interested in results that are as broad and robust as possible, and in flexible and insight-driven frameworks that I can adapt to analogous situations.

Oh, and also: RIP Danny Segal, the amazing physicist who taught me about Lamb-Dicke sidebands.
 A: I would formulate the quantum problem somewhat differently with
a quantized electromagnetic field, but in the
end with either method, the quantum transition rate between 2 energy eigenstates
of the unperturbed Hamiltonian will be given by Fermi's golden rule.
The energy delta function will enforce the photon frequency
$\omega$ to be  $\omega = k c = \omega_0+\nu(n-m)$,
where $n$ is the initial harmonic oscillator level and $m$ is the final.
The addition of the harmonic well, as you point out, gives an
additional factor of the oscillator
matrix element $\langle m|e^{i\eta(a+a^+)}|n\rangle$ in the Fermi's
golden rule matrix element. The operator is the momentum transfer operator,
so this is the amplitude for a harmonic oscillator state $n$ to have the
photon momentum added and to then overlap with the oscillator
state $|m\rangle$. In the position representation, this is the
Fourier transform of the product of the two oscillator wave functions.
For a more general well, calculating that integral
would likely be the easiest way to proceed.
For the harmonic oscillator case, this matrix element can be calculated
analytically. In position representation, the oscillator energy eigenstates are
gaussians time Hermite polynomials. The Fourier transform of the product
is therefore a gaussian times a polynomial in $k$. These integrals
are given in Gradsteyn and Ryzhik, Tables of integrals series and
products,
7.388.6 and 7.388.7. Alternatively,
you can work in the second quantized notation you used and write
\begin{equation}
e^{i\eta(a+a^+)} =e^{\frac{\eta^2}{2}} e^{i\eta a}e^{i\eta a^+}
\end{equation}
so that expanding the exponentials
\begin{equation}
\langle m|e^{i\eta(a+a^+)}|n\rangle =
e^{\frac{\eta^2}{2}} \frac{1}{\sqrt{m!n!}}
\sum_j\frac{(i\eta)^j}{j!}
\sum_{j'}\frac{(i\eta)^{j'}}{j'!}
\langle 0|
 (a)^{j+m}
 (a^+)^{j'+n}|0\rangle \,.
\end{equation}
This will be zero unless $j+m=j'+n$, and the remaining sum is 
a series in the form of a
confluent hypergeomtric function, Abramowitz and Stegun, Handbook of
mathematical functions, Eq. 13.1.2,
\begin{equation}
\langle m|e^{i\eta(a+a^+)}|n\rangle
= e^{\frac{\eta^2}{2}}\frac{1}{\sqrt{m!n!}}\frac{n!}{(n-m)!}
\left (i\eta\right)^{n-m}
M\left (n+1,n-m+1,-\eta^2\right)
\end{equation}
where for simplicity I have assumed that $n>m$, if $n<m$, interchange them.
This can be simplified by using the Kummer transformation, Abramowitz
and Stegun, Eq. 13.1.27
\begin{equation}
M(a,b,z)=e^zM(b-a,b,-z)
\end{equation}
and the special case Eq. 13.6.9 that
$M(-n,\alpha+1,x)=\frac{n!}{(\alpha+1)_n}L^{(\alpha)}_n(x)$, where
$L^{(\alpha)}_n(x)$ is a associated Laguerre polynomial. These Laguerre
polynomial equations are the result given in Gradshteyn and Ryzhik.
The classical limit is when both quantum numbers are much greater
than 1 so that the energy quantum is much smaller than the energy.
We can use the expansion
Abramowitz and Stegun, Eq. 13.3.7,
\begin{equation}
\frac{M(a,b,z)}{\Gamma(b)} = e^{\frac{z}{2}}
\left (\tfrac{1}{2}bz-az\right )^{\frac{1}{2}-\frac{b}{2}}
\sum_{n=0}^\infty A_n\left (\frac{z}{2}\right)^{\frac{n}{2}}
(b-2a)^{-\frac{n}{2}} J_{b-1+n}\left (\sqrt{2zb-4za}\right )
\end{equation}
where $A_0=1$, $A_1=0$, $A_2=\frac{b}{2}$,
$(n+1)A_{n+1}=(n+b-1)A_{n-1}+(2a-b)A_{n-2}$.
Applied to the matrix element, we see the higher order terms are down by
powers of $1/m$, so they can be ignored in the classical limit.
The matrix element becomes
\begin{eqnarray}
\langle m|e^{i\eta(a+a^+)}|n\rangle =
&=& \sqrt{\frac{n!}{m!}}
\left (i\sqrt{\frac{2}{n+m+1}}\right )^{n-m}
J_{n-m}\left (\sqrt{2}\eta\sqrt{n+m+1}\right)+...
\end{eqnarray}
Again taking $m$ and $n$ large, the factors in front of the Bessel function
go to 1, and this goes to
\begin{equation}
\langle m|e^{i\eta(a+a^+)}|n\rangle =
\rightarrow i^{n-m} J_{n-m}\left (\eta\sqrt{2}\sqrt{n+m+1}\right)
\end{equation}
To include the case where $m>n$, this becomes
\begin{equation}
\langle m|e^{i\eta(a+a^+)}|n\rangle =
\rightarrow i^{|n-m|} J_{|n-m|}\left (\eta\sqrt{2}\sqrt{n+m+1}\right)
\end{equation}
In the classical limit, we can identify
$\frac{1}{2} m\nu^2 x_0^2= \frac{\hbar \nu (n+m+1)}{2}$
as the energy of the classical oscillator, and with
$\eta = k\sqrt{\frac{\hbar}{2m\nu}}$, so that $\eta\sqrt{2}\sqrt{n+m+1}=kx_0$,
the result is
\begin{equation}
\langle m|e^{i\eta(a+a^+)}|n\rangle =
\rightarrow i^{|n-m|} J_{|n-m|}\left (kx_0)\right) \,,
\end{equation}
This agrees with your classical calculation with the 
$k$ value modified appropriately to conserve energy. The matrix element calculation is exact up until the truncation of the Bessel function expansion, so keeping more terms will bring back more quantum effects.
This shows that mainly you just need to calculate the matrix element
of $\langle m|e^{i\vec k \cdot \vec r}|n\rangle$ and take the large
quantum number limit in order to recover the classical result. 
For arbitrary potentials you might use a WKB method to get semiclassical
states and compute the Fourier transform with a stationary phase/steepest
descents method.
Alternatively, you could formulate the motion using a coherent state basis, which will
give a more direct connection between the classical and quantum results.
