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.