How can one recover the classical frequency-modulation Bessel sidebands from a quantum emitter in a harmonic well?

Question

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.}

Answer 1

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.