The Quantum Spectral Method: From Atomic Orbitals to Classical Self-Force

Majed Khalaf
Racah Institute of Physics, Hebrew University of Jerusalem, Jerusalem 91904, Israel

Ofri Telem
Racah Institute of Physics, Hebrew University of Jerusalem, Jerusalem 91904, Israel

We present the Quantum Spectral Method for the analytical calculation of observables in classical periodic and quasi-periodic systems. It is based on a novel application of the correspondence principle, in which classical Fourier coefficients are obtained as the $\hbar\rightarrow 0$ limit of quantum matrix elements. Our method is particularly suited for self-force calculations for inspiralling bodies, where it could, for the first time, provide fully analytical expressions. We demonstrate our method by calculating an adiabatic electromagnetic inspiral, along with its associated radiation, at all orders in the multipole expansion.

§ INTRODUCTION

Classical physical systems exhibiting periodic motion have been a centerpiece of physics throughout its centuries-old history - the quintessential examples being the harmonic oscillator and Keplerian motion. Early in the 20th century, these prototypes gave rise to two invaluable extensions: (a) their generalization to classical quasi-periodic motion using adiabatic invariants; and (b) their quantization in the form of the quantum harmonic oscillator, and the hydrogen-like atom, respectively.

In recent years, the subject of adiabatic and post-adiabatic (PA) corrections to periodic motion has gone through an experimental and computational renaissance, following the advent of gravitational-wave signals from inspiralling black holes <cit.>. The application of PA perturbation theory is particularly useful in the extreme mass-ratio regime in which a small compact object spirals into a supermassive black hole more than $10^6$ its mass. In this regime, the relative change in the fundamental frequency over a single period is very small, and the next-to-leading order in the PA expansion (1PA) is enough for the adequate description of the system. This seemingly innocuous statement requires some clarification; the application of next-to-leading order PA perturbation theory to gravitational inspirals is notoriously hard. It requires the calculation of the second-order self-force due to the metric back-reaction from the inspiralling body. This in turn necessitates the iterative solution of Einstein's equations with the inspiralling body (and the adiabatic metric perturbation) as a source.
Though adiabatic (0PA) codes are reasonably efficient, the computation of the second-order self-force tables, required for a full 1PA computation, is highly taxing. In particular, due to the non-linearity of general relativity, the effective source for the second-order gravitational field includes a term bilinear in the first-order gravitational field <cit.>. It is the multipole expansion of this bilinear term that leads to a significant computational bottleneck, due to its slow convergence in the multipole index $l$. Though some methods are known for mitigating this problem, modern codes for the generation of second-order self-force codes could take 2-3 days to complete for each data point, even on a 40-core machine <cit.>. Though we do not address this bottleneck directly in this work, our analytical first-order self-force computation could serve as a baseline for enhancing the convergence of the second-order computation, by allowing for a partial resummation of its bilinear contribution. This potential application to second-order self-force will be discussed in a follow-up work.

The aim of the current paper is to present a highly physical approach for the analytical computation of observables in classical (quasi)-periodic systems, of which the source integrals for inspirals are a prime example. Our method, which we call the "Quantum Spectral Method" (QSM), utilizes the correspondence principle between quantum and classical physics for the calculation of classical Fourier coefficients from quantum matrix elements. The marriage of adiabatic invariants and quantum mechanics is not new; in fact, it played a significant role in the development of old quantum theory <cit.>. Nevertheless, we believe that our particular application of the correspondence principle is new, and generates previously unknown analytical results.

The essence of our method is as follows. Consider a system periodic in an angle variable $\alpha$ with a single fundamental frequency $\Upsilon$. If the system is periodic in time, we can take $\alpha$ to be linear[Since $\alpha$ only enters as an angle, we can allow it to grow linearly, keeping in mind that $\alpha+2\pi\sim\alpha$.] $\alpha=\Upsilon (t-t_0)$; otherwise, if the system is quasi-periodic, $\alpha-\Upsilon (t-t_0)$ could be slowly changing. $\alpha$-dependent observables $\mathcal{O}(\alpha)$ are conveniently expressed as Fourier series

\begin{eqnarray}\label{eq:rtQSM3int}
\mathcal{O}(\alpha)&=&\,\sum_{\Delta n=0}^\infty\,\left[\mathcal{O}^c_{\Delta n}\,\cos\left(\Delta n\,\alpha\right)\,+\,\mathcal{O}^s_{\Delta n}\sin\left(\Delta n\,\alpha\right)\,\right]\,,
\end{eqnarray}

where the $\mathcal{O}_{\Delta n}$ are Fourier coefficients. These coefficients are the solutions to the Fourier-space classical equations of motion (EOM) of the system. In many next-to-trivial cases, the Fourier-space EOM do not have known analytical solutions, and one must resort to numerical integrals (see, e.g. <cit.>).

In this paper, we present a novel method for the analytic calculation of the Fourier coefficients $\mathcal{O}_{\Delta n}$. This method gives exact analytical solutions, even in cases where only numerical integration has been available before. We explicitly demonstrate our method by reproducing (a) time-dependent Keplerian motion; (b) the all-multipole expression for the electromagnetic (EM) radiation from a classical electron in Keplerian motion; and (c) an adiabatic EM inspiral of a classical electron, and its associated waveform.

Our method is based on the realization of $\mathcal{O}_{\Delta n}$ as the classical limit of sums of quantum matrix elements, taken between Hamiltonian eigenstates. This realization is a direct implementation of the correspondence principle, and is also verified explicitly. More specifically, we show that the classical Fourier coefficients are given by

\begin{eqnarray}\label{eq:masterin}
&&~~~~~~~~~~~~~~~~~\mathcal{O}^c_{\Delta n}={\rm Re}\,\mathcal{O}_{\Delta n}~~,~~\mathcal{O}^s_{\Delta n}={\rm Im}\,\mathcal{O}_{\Delta n}\nonumber\\[5pt]
&&\mathcal{O}_{\Delta n}\,\equiv\,\lim_{\hbar\rightarrow 0}\,\begin{cases}~~\sum_{\Delta \sigma}\left\langle n,\sigma-\Delta \sigma\right|\,\mathcal{O}\,\left|n,\sigma\right\rangle & \Delta n = 0\\2\sum_{\Delta\sigma}\left\langle n-\Delta n,\sigma-\Delta \sigma\right|\,\mathcal{O}\,\left|n,\sigma\right\rangle & \Delta n \neq 0\end{cases}\,,
\end{eqnarray}

where $\mathcal{O}$ inside the quantum matix element is regraded as an operator, $\left| n,\sigma\right\rangle$ are the Hamiltonian eigenstates,   $n$ is the principal quantum number, and $\sigma$ denotes the non-principal quantum numbers. Moreover, $(n,\sigma)=\hbar^{-1}(N,\Sigma)$, where $(N,\Sigma)$ are the classical action variables of the system.

The non-trivial applications of QSM presented in this paper all concern the motion of a bound classical electron in a Coulomb potential, with or without radiation-reaction (c.f. the study of the dissipative motion of a scattering electron in <cit.>). Nevertheless, the most important future application of the QSM is within the  PA approach <cit.> for gravitational inspiral problems  in the extreme-mass-ratio regime. These are crucial for the theoretical calculation of gravitational waveforms for the LISA future space-based gravitational wave detector <cit.>. In this regime, the computation of the gravitational self-force tables at first <cit.> and second-order <cit.> required for the accurate calculation of the inspiral <cit.> is a major computational bottleneck <cit.>. This is in contrast with the online step of waveform generation, which is largely a solved problem <cit.>. Using the QSM might allow to replace some or all of the numerical source integrals used in the calculation of self-force tables with analytical expressions, potentially allowing for a speedup of the computation at first and even second order.

Considering that there are exact and compact analytical expressions for the EM self-force, given primarily by the Abraham-Lorentz-Dirac (ALD) <cit.> and Landau-Lifshitz (LL) forces <cit.> (see <cit.> for a modern, rigorous derivation), it is important to clarify our choice of computing the self-force perturbatively, as part of the PA expansion. First, we ultimately envision the application of the QSM to the gravitational case, where the non-linearity of gravity necessitates the use of PA perturbation theory. The EM example presented here then serves as an important first step towards this goal. Secondly, the non-linearity of the ALD force makes it impractical for the computation of the inspiralling the trajectory over the course of a large number of cycles, as required for EMRI. Finally, the QSM is applicable beyond EM (e.g., in cases with forces mediated by scalar fields). Consequently, in order to not limit our scope, we opt not to utilize the LL equation, as it is typically used when the external forces are electromagnetic.

The place of the QSM within the landscape of quantum-to-classical methods is illustrated in Fig. <ref>.
The grand majority of the current quantum-to-classical methods are based on taking the classical limit of perturbative Quantum Field Theory (QFT) scattering amplitudes, for example in the KMOC method <cit.>. In relativistic Effective Field Theory (EFT)
<cit.>, the classical scattering angle is then translated to bound state data via the Boundary-to-Bound map <cit.>. In contrast, the non-relativistic EFT (NREFT) <cit.> and On-Shell <cit.> approaches utilize the classical limit of QFT scattering amplitudes to perform EFT matching to an effective two-body Hamiltonian, which is subsequently used either to calculate inspirals, or as input to an effective one-body Hamiltonian <cit.>. The Classical Bethe-Salpeter <cit.> approach, on the other hand, aims to go directly from QFT amplitudes to bound state data, via an all-order resummation in the soft limit. Lastly, there has been growing interest in the self-force EFT approach - the effective field theory for the calculation of two-body gravitational interactions in the extreme mass ratio limit. Within this framework, an effective action for self-forced motion is constructed using a curved-background EFT (e.g., <cit.>) or by the resummation of post-Minkowskian Feynman diagrams using a flat background EFT<cit.>.

The Landscape of Quantum-to-Classical Methods. The methods indicated in the plot are KMOC <cit.>, relativistic EFT
<cit.>, NREFT <cit.>, On-Shell methods <cit.>, Boundary-to-Bound <cit.>, self-force EFT <cit.>, Classical Bethe-Salpeter <cit.>, original Bethe-Salpeter <cit.>, and the QSM (this work).

Compared to the existing literature, our approach traces an original path; in our work, the classical limit is taken at the level of the quantum bound states, and so is always non-perturbative in the coupling. This makes the QSM particularly suited for the quantum-to-classical map in the context of self-force/PA perturbation theory. The all-order, probe limit quantum-to-classical map was previously explored by one of the current authors in <cit.>, in the context of scattering, rather than bound-state motion. The current work extends and formalizes the ideas of <cit.>, and applies them to bound states and radiation reaction.

The paper is structured as follows. In Section <ref> we present the master equation of the QSM and its proof for motion in a $1/r$ potential. We then present three applications of the QSM in increasing levels of sophistication. First, we reproduce the known series solution for time-dependent Keplerian motion in Section <ref>. Next, in Section <ref>, we present the first application of the QSM to a problem whose analytical solution has not been previously known: the all-order emission (in the multipole expansion) of EM waves from a classical electron in a Keplerian orbit. In Section <ref>, we use the QSM to compute an adiabatic (0PA) EM inspiral trajectory and its accompanying EM wave signal. Remarkably, we find that in the adiabatic limit (0PA), the loss rate for energy and angular momentum is the classical limit of the rate for quantum spontaneous emission.
Finally, Section <ref> presents our conclusions and blueprint for a future application of the QSM to systems with multiple fundamental frequencies, non-periodic systems, and to gravitational self-force problems. To enhance readability, we relegate some of the more technical aspects of our computations to our appendices. In particular, Appendix <ref> provides a more detailed proof of the QSM master formula; Appendix <ref> is a compendium of quantum matrix elements for hydrogen-like atoms; Appendix <ref> presents all of the classical limits of the matrix element from Appendix <ref>; Appendix <ref> is the derivation of the 0PA classical loss rate as the classical limit of spontaneous emission; and Appendix <ref> contains auxilliary computations.

§ THE QUANTUM SPECTRAL METHOD WITH A SINGLE FUNDAMENTAL FREQUENCY

Consider a classical system undergoing (quasi)-periodic motion. For simplicity, we consider the case in which the system has a single fundamental frequency  - though the QSM works just as well for systems with multiple fundamental frequencies. In this case, every classical observable $\mathcal{O}(\alpha)$ of the system can be expressed as a simple Fourier series (<ref>), where $\mathcal{O}_{\Delta n}$ are Fourier coefficients that depend on the system's dynamics. We denote the integer "harmony" by $\Delta n$ for reasons that will become clear momentarily.
The ”master equation" of the QSM is the statement that the Fourier coefficients $\mathcal{O}_{\Delta n}$ are the direct classical limits of quantum matrix elements, taken between eigenstates of the quantum Hamiltonian of the system. Specifically, consider a set of eigenstates of the Hamiltonian,

\begin{eqnarray}
\end{eqnarray}

Here, $n$ is the principal quantum number, and $\sigma$ are other quantum numbers. The QSM ”master equation" is then given by (<ref>), and we repeat it here for emphasis:

\begin{eqnarray}\label{eq:master}
\mathcal{O}_{\Delta n}\,=\,\lim_{\hbar\rightarrow 0}\,\begin{cases}~~\sum_{\Delta \sigma}\left\langle n,\sigma-\Delta \sigma\right|\,\mathcal{O}\,\left|n,\sigma\right\rangle & \Delta n = 0\\2\sum_{\Delta\sigma}\left\langle n-\Delta n,\sigma-\Delta \sigma\right|\,\mathcal{O}\,\left|n,\sigma\right\rangle & \Delta n \neq 0\end{cases}\,,
\end{eqnarray}

where $(n,\sigma)=\hbar^{-1}(N,\Sigma)$, and $N,\Sigma$ are the classical action variables of the system. In particular, $N$ is related to the classical energy $E$.
Furthermore, the classical fundamental frequency corresponding to $\alpha$, which we denote by $\Upsilon$, is related to the energy eigenstates as

\begin{eqnarray}\label{eq:masterfreq}
\Upsilon=\lim_{\hbar\rightarrow 0}\frac{E_n-E_{n-\Delta n}}{\hbar\Delta n}\,.
\end{eqnarray}

In the rest of this section we outline a direct proof of the QSM master equation for the particular case of motion in a $1/r$ potential. The classical version of the system is Keplerian motion, while the quantum version is the hydrogen atom. Note that while this is an energy conserving, non-radiating system, in the last section the QSM is directly applied to a radiating system as part of PA perturbation theory.

§.§ Proof of the QSM for a $1/r$ Potential

Our proof of the QSM for a $1/r$ potential relies on the use of coherent states, and taking their classical limit. The use of coherent states for the quantum-to-classical map is not new, and has been used before <cit.> in the context of the Eikonal exponentiation for spinning particles <cit.>. Here we use it directly for the hydrogen-like atom.

In the classical limit, the specifics of our choice of hydrogen-atom  coherent states become insignificant. Nevertheless, for concreteness we choose a particular set of coherent states, detailed below. Given a particular choice of coherent states, any classical observable is expressible by definition as

\begin{eqnarray}
\mathcal{O}(\alpha)\,=\,\lim_{\hbar\rightarrow 0}\,\left\langle {\rm coh}\right|\,\mathcal{O}\,\left|{\rm coh}\right\rangle\,.
\end{eqnarray}

The master equation of the QSM is then derived by choosing a particular set of coherent states and explicitly taking the classical limit. The Hamiltonian of the system is the famous Hamiltonian of the hydrogen-like atom

\begin{eqnarray}
\end{eqnarray}

where $K=\frac{Zq^2}{4\pi\varepsilon_0}$, $\mu$ is the mass, and we take $\varepsilon_0=1$.
We choose the coherent states defined by Crawford <cit.> based on the famous construction by Klauder <cit.>. Specializing to coherent states corresponding to classical elliptical motion in the XY plane, with the major axis along X, these are[Note our change of sign in $m_1$ and $m_2$ with respect to Crawford's definition. This is reflective of our definition (<ref>) of $\eta$, while Crawford chooses the positive sign of the square root.]

\begin{eqnarray}\label{eq:cohcraw}
\left|{\rm coh}\right\rangle=\left|N,\eta,\alpha\right\rangle=&&\mathcal{A}(N/\hbar)\sum_{n,l,m}\,g_n(N,\alpha)\,f_{l,m}\left(\tfrac{n-1}{2},\eta\right)\left|n,l,m\right\rangle\nonumber\\[5pt]
g_n(N,\alpha)\equiv &&\,\frac{n\,(N/\hbar)^{n-1}\exp\left[{-i\frac{\alpha E_{n}}{\Upsilon\hbar}}\right]}{\sqrt{(n-1)!}}\nonumber\\[5pt]
f_{l,m}(j,\eta)\equiv &&\,\sum_{m_1,m_2}\frac{(2j)!}{\sqrt{(j-m_1)!(j+m_1)!(j-m_2)!(j+m_2)!}}\,\frac{(-1)^{j-m_2}\eta^{2j-m_1-m_2}}{(1+\eta^2)^{2j}}\,C^{lm}_{jm_1;j,m_2}\,.
\end{eqnarray}

Here $\left|n,l,m\right\rangle$ are the usual quantum states of the hydrogen atom, $E_n=-\frac{\mu K^2}{2\hbar^{2}n^2}$ are its energy levels, and $C^{lm}_{jm_1;j,m_2}$ are Clebsch-Gordan coefficients. The normalization is $\mathcal{A}(N/\hbar)=e^{-N/2\hbar}(1+3\hbar^{-1}N+\hbar^{-2}N^2)^{-1/2}$. Here $N$, $\Upsilon$ and $\eta$ are related to the classical $E,e$ via

\begin{eqnarray}\label{eq:Eeta}
E=-\frac{\mu K^2}{2N^2}~~~,~~~\Upsilon=\frac{2E}{N}~~~,~~~
\eta=\frac{1-\sqrt{1-e^2}}{e}\,,
\end{eqnarray}

where $e$ is the eccentricity, defined in the next section. The state is constructed so that $t\rightarrow t+\Delta t$ upon Hamiltonian time evolution. Taking the expectation value with respect to these coherent states, we have

\begin{eqnarray}\label{eq:rtQSM2}
\mathcal{O}(\alpha)=&&\lim_{\hbar\rightarrow 0}\left\langle N,\eta,\alpha\right|\mathcal{O}\left|N,\eta,\alpha\right\rangle=\sum_{\Delta n,\Delta l,\Delta m}\,\lim_{\hbar\rightarrow 0}\,\mathcal{A}^2(N)\sum_{n,l,m}\,g^*_{n'}(N,\alpha)\,g_n(N,\alpha)\,\nonumber\\[5pt]
&&\times f^*_{l',m'}\left(\tfrac{n-\Delta n-1}{2},\eta\right)f_{l,m}\left(\tfrac{n-1}{2},\eta\right)\left\langle n',l',m'\right| \mathcal{O}\left|n,l,m\right\rangle\,,
\end{eqnarray}

where $(n',l',m')=(n,l,m)-(\Delta n,\Delta l,\Delta m)$. As we take the $\hbar\rightarrow 0$ limit, the $g_n,\,g_{n'}$ and $f_{l,m},\,f_{l',m'}$ factors single out only the leading contribution from the $n,\,l,\,m$ sums, namely the leading order Laplace approximation becomes exact. This fixes

\begin{eqnarray}\label{eq:fixclas}
&&n\rightarrow N\hbar^{-1}~~,~~l,m\rightarrow L\hbar^{-1}\,.
\end{eqnarray}

In this way we get

\begin{eqnarray}\label{eq:rtQSM3}
\mathcal{O}(\alpha)&=&\,\sum_{\Delta n}\,\exp\left[{-i\Delta n\,\alpha}\right]\,\mathcal{O}_{\Delta n}\nonumber\\[5pt]
\mathcal{O}_{\Delta n}&\equiv&\lim_{\hbar\rightarrow 0}\,\sum_{l',m'}\,\left\langle n'l'm'\right| \mathcal{O}\left|nlm\right\rangle\,,
\end{eqnarray}

where the matrix element is taken between $(n,l,m)=\hbar^{-1}(N,L,L)$ and $(n',l',m')=(n,l,m)-(\Delta n,\Delta l,\Delta m)$. The time dependence in this expression comes from $g^*_{n'}\,g_n$. For a more thorough derivation of (<ref>), see Appendix <ref>. As a classical Fourier coefficient, ${\lim_{\hbar\rightarrow 0}\rm Re}\left\langle\mathcal{O}_{\Delta n}\right\rangle$ is symmetric under $\Delta n\rightarrow-\Delta n$, while $\lim_{\hbar\rightarrow 0}{\rm Im}\left\langle\mathcal{O}_{\Delta n}\right\rangle$ is antisymmetric, and so finally, we get the master equation (<ref>).

§ FIRST APPLICATION OF THE QSM: TIME-DEPENDENT KEPLERIAN MOTION

In this section we demonstrate the QSM by applying it to one of the oldest problems in physics: the time dependent motion of a classical object in a $1/r$ potential - i.e. Keplerian motion. For concreteness we take this object to be a classical electron moving in the Coulomb potential of a classical nucleus.

§.§ The Time-Dependent Solution to the Kepler Problem

The EOM governing the motion in a central Coulomb potential is

\begin{eqnarray}\label{eq:EOM}
\mu\ddot{\vecc{r}}=-K\frac{\vecc{r}}{r^3}\,,
\end{eqnarray}

where $\mu$ is the reduced mass, $K=\frac{Zq^2}{4\pi\varepsilon_0}$, and we take $\varepsilon_0=1$. Owing to spherical symmetry and angular momentum conservation, the motion is planar. The EOM is separable into a radial and an azimuthal part:

\begin{eqnarray}\label{eq:EOMrph}
\mu\ddot{r}+\frac{K}{r^2}-\frac{L^2}{\mu r^3}=0~~,~~\dot{\varphi}=\frac{L}{\mu r^2}\,,
\end{eqnarray}

where $L$ is a separation constant - the conserved angular momentum of the system. Once the radial part is solved, one can solve for the azimuthal part as well. To do this, we define the Kepler orbital function

\begin{eqnarray}\label{eq:orbitf}
\mathfrak{r}(p,e,\psi)\equiv\frac{p}{1+e\cos(\psi)}\,.
\end{eqnarray}

For fixed $p$ and $e$, $\mathfrak{r}(\psi)$ is an ellipse with semi-latus-rectum $p$ (with units of length) and eccentricity $e$. Substituting the ansatz $r(t)=\mathfrak{r}(p,e,\psi(t))$ in (<ref>), we get the EOM for $\psi(t)$,

\begin{eqnarray}\label{eq:eqpsir}
\dot{\psi}=\frac{L}{\mu r^2}~~~,~~~p=\frac{L^2}{K\mu}~~~,~~~e=\sqrt{1+\frac{2EL^2}{K^2\mu}}\,.
\end{eqnarray}

The pair $(p,e)$ is collectively called the "orbital elements" of the Keplerian motion[In the presence of a non-planar perturbing force, it is also crucial to keep note of the orientation of the Keplerian ellipse with respect to a reference plane. In this paper we only deal with a planar perturbing force, and so we can focus only on $(p,e)$].
Note that for pure Keplerian motion, the EOM (<ref>) and (<ref>) for $\varphi$ and $\psi$ are exactly of the same form and so they can be identified. This would not generically be the case for perturbed Keplerian motion[In some cases, however, this identification still holds even for perturbed Keplerian motion - see footnote 7 below and the discussion in Section <ref> for more details.].
The time dependence of $\psi(t)$ in equation (<ref>) is solved explicitly by changing variables to the Mean Anomaly[Sometimes denoted by $M$ in the literature] $\alpha(t)$, which is the action-angle suitable for the Kepler problem. This change of variables is given by

\begin{eqnarray}\label{eq:AAKep}
\alpha\equiv\beta-e\sin\beta~~~,~~~\beta\equiv2\arctan\left[\sqrt{\frac{1-e}{1+e}}\tan(\psi/2)\right]\,.
\end{eqnarray}

The angle $\beta$ is known as the Eccentric Anomaly[Sometimes denoted by $E$ in the literature.], and the first of these equations is known as the Kepler Equation. $\Upsilon$ is the fundamental frequency of the Kepler problem, and is given by

\begin{eqnarray}\label{eq:Upsilon}
\Upsilon\equiv\sqrt{\frac{K(1-e^2)^3}{\mu p^3}}\,.
\end{eqnarray}

One can check that this definition is consistent with (<ref>). Since $r$ and $\varphi$ (up to integer multiples of $2\pi$) are periodic in $\alpha$, we can represent them as Fourier series. Following <cit.>, we have

\begin{eqnarray}\label{eq:reKwpclas}
r(\alpha)&=&\frac{p}{1-e^2}\left[1+\frac{e^2}{2}-2e\sum_{\Delta n=1}^\infty\,\frac{1}{\left(\Delta n
\right)^2}\frac{d J_{\Delta n}(\Delta n \,e)}{de}\cos\left(\Delta n \alpha\right)\right]\nonumber\\[5pt]
\psi(\alpha)&=&2\arctan\left\{\sqrt{\frac{1+e}{1-e}}\tan\left[\beta(\alpha)/2\right]\right\}\nonumber\\[5pt]
\beta(\alpha)&=&\alpha+2\sum_{\Delta n=1}^\infty\,\frac{1}{\Delta n
}J_{\Delta n}(\Delta n \,e)\,\sin\left(\Delta n \alpha\right)\,.
\end{eqnarray}

Here $J$ is a Bessel function of the first kind and we choose the symbol $\Delta n$ to label the integer Fourier harmonics, for reasons that would become clear in the next section. Note again that for pure Keplerian motion $\varphi=\psi$.
All that is left is to find the time dependence of $\alpha$. Going through the change of variables in (<ref>), it is straightforward to check that

\begin{eqnarray}\label{eq:alphadotKep}
\dot{\alpha}=\Upsilon\,,
\end{eqnarray}

and so $\alpha=\Upsilon(t-t_0)$. Equation (<ref>) is also true for perturbed Keplerian motion in the PA expansion <cit.>, as long as $\alpha$ is properly chosen[In other words, $\alpha$ is the angle variable after a near-identity transformation (NIT) with no frequency corrections - see section 2.6.4 of <cit.>.]. Note that in this case $\Upsilon$ has a slow time dependence due to the adiabatic change in $p(t)$ and $e(t)$.

§.§ Time-Dependent Keplerian Motion with the QSM

Using the master equation for the radius $r(\alpha)$, we have

\begin{eqnarray}\label{eq:masterKep}
r(\alpha)&=&\,\sum_{\Delta n=0}^\infty\,r^c_{\Delta n}(N,L)\,\cos\left(\Delta n\,\alpha\right)\nonumber\\[5pt]
r^c_{\Delta n}(N,L)\,&=&\,\lim_{\hbar\rightarrow 0}\,\begin{cases}~~\sum_{\Delta l}\left\langle n',l',m'\right|\,r\,\left|n,l,m\right\rangle & \Delta n =0\\
2\sum_{\Delta l}\left\langle n',l',m'\right|\,r\,\left|n,l,m\right\rangle & \Delta n \neq 0\end{cases}\,.
\end{eqnarray}

Here $(n,l,m)=\hbar^{-1}(N,L,L)$ and $(n',l',m')=(n,l,m)-(\Delta n,\Delta l,\Delta l)$. The matrix element is real and so we only have cosines in our Fourier series. The $r$ matrix element and its classical limit are calculated in Appendix <ref>, with the result

\begin{eqnarray} \label{ClassLimRadialtext}
&&r_{\Delta n}(N,L)=\frac{p}{1-e^2}\begin{cases}1+\frac{\,e^2}{2}& \Delta n=0\\-\frac{2e}{\Delta n^2}\frac{d}{de}J_{\Delta n}(e\Delta n)&\Delta n\neq 0\end{cases}\,.
\end{eqnarray}

Here $p,\,e$ are again the semi latus rectum and eccentricity from  (<ref>). Substituting this in (<ref>), we reproduce the classical Fourier series for $r(\alpha)$ (<ref>). The series for $\beta(\alpha)$ can be reproduced in a similar manner.

§ ALL-ORDER EM RADIATION WITH THE QSM

In our treatment of the Kepler problem, we were fortunate enough to have the known Fourier-series solution (<ref>) to compare with. In this section we present a first application of the QSM to a classical problem whose all-order analytical solution, we believe, is not previously known. In fact, the conventional treatment involves a reduction to a series of radial integrals that require numerical evaluation.

The problem we address in this section is the calculation of the retarded EM field $A^\mu_{\rm ret}$ generated by a classical electron moving along a (quasi-) Keplerian orbit. By quasi-Keplerian, we mean that the trajectory is given by (<ref>), but allowing for a slow change in $\alpha-\Upsilon (t-t_0)$ over time. When $\dot{\alpha}$ is constant, the motion is purely Keplerian.

A classical non-relativistic[We leave the consideration of relativistic effects to future work. This involves considering fine-structure corrections to the hydrogen orbitals on the QSM side. This is not hard to do, and we omit it here for simplicity. Note that our results will be all-order in the multipole expansion, which also involves $v/c$.] electron on a quasi-Keplerian orbit generates an electric current density given by

\begin{eqnarray}
\end{eqnarray}

where $v^\mu_{\rm{Kep}}(t')=\left(1,\partial_{t'}\vec{r}_{\rm{Kep}}\right)$ and $q$ is the electric charge of the electron. Here we use the label $\rm{"Kep"}$ to indicate the position and velocity of the source, which is in a (quasi-) Keplerian orbit. In this section we assume that the Keplerian/quasi-Keplerian orbit of the electron is known, namely that $\vec{r}_{\rm Kep}(t')$ are given as an input. The Lorenz gauge EM field generated by this current is the solution to the sourced wave equation,

\begin{eqnarray}
\square\,A^\mu_{\rm ret}=J^\mu\,.
\end{eqnarray}

The solution is an integral of the retarded EM Green's function with the source:

\begin{eqnarray}\label{eq:greensol}
A^\mu_{\rm ret}(t,\vec{x})=\int\,d^4x'\,G^{{\rm ret};\mu}_{~~~~~\nu}(t,\vec{x};t',\vec{x}')\,J^\nu(t',\vec{x}')\,.
\end{eqnarray}

The retarded Green's function for the EM field is famously <cit.>

\begin{eqnarray}~\label{eq:GreenretPS}
G^{\rm ret}_{\mu\nu}(t,\vec{x};t',\vec{x}')=\frac{\Theta(t-t')}{4\pi R}\,\delta(t-t'-R)\,,
\end{eqnarray}

where $R=|\vec{x}-\vec{x}'|$. This is the Green's function used to derive the Liénard-Wiechert retarded potential in electromagnetism. For our purposes it's convenient to use the Fourier representation of the delta function

\begin{eqnarray}~\label{eq:GreenretFour}
G^{\rm ret}_{\mu\nu}(t,\vec{x};t',\vec{x}')=\frac{\Theta(t-t')}{2\pi}\,\int_{-\infty}^{\infty}\,d\omega \,e^{-i\omega(t-t')}\,\frac{e^{i\omega R}}{4\pi R}\,.
\end{eqnarray}

We now expand $e^{i\omega R}/(4\pi R)$ in multipoles and get

\begin{eqnarray}~\label{eq:GreenretMulti}
G^{\rm ret}_{\mu\nu}(t,\vec{x};t',\vec{x}')=\frac{\Theta(t-t')}{2\pi}\,\int_{-\infty}^{\infty}\,d\omega \,e^{-i\omega(t-t')}\,\left\{i\omega\,\sum_{l=0}^\infty\,j_{l}(\omega r_<)\,h^{(1)}_l(\omega r_>)\,\sum_{m=-l}^l\,Y_l^{m*}(\theta',\varphi')Y_l^{m}(\theta,\varphi)\right\}\,,
\end{eqnarray}

where $\{r_<,r_>\}=\{{\rm min}(r,r'),{\rm max}(r,r')\}$, $j_{l_\gamma}$ is a spherical Bessel function, and $h^{(1)}_{l_\gamma}$ is a spherical Hankel function of the first kind. Using the QSM, we can express $A^\mu$ as the classical limit of an expectation value over our coherent states,

\begin{eqnarray}\label{eq:greencoh}
A^\mu_{\rm ret}(t,\vec{x})&=&\lim_{\hbar\rightarrow 0}\,\int\,d^4x'\,\left\langle N,\eta,\alpha'(t')\right|G^\mu_{~\nu}(t,\vec{x};t',\vec{x}')\,J^\nu\left|N,\eta,\alpha'(t')\right\rangle\,\nonumber\\[5pt]
&=&\frac{q}{\mu}\,\int\,dt'\,\lim_{\hbar\rightarrow 0}\,\left\langle N,\eta,\alpha'(t') \right|\,G^\mu_{~\nu}(t,\vec{x};t',\vec{r}_{\rm{Kep}})\,p^\nu_{\rm{Kep}}\left|N,\eta,\alpha'(t')\right\rangle\,.
\end{eqnarray}

Note that the expectation value is only over $\vec{r}_{\rm{Kep}}$ and $p^\mu_{\rm{Kep}}=\mu\,(1,\vec{v}_{\rm{Kep}})$ regarded as position/momentum operators acting on the coherent states, while $t, \vec{x}$ and $t'$ are c-numbers. The time dependence of the coherent states is in terms of $\alpha'(t')$, considered as input for the calculation.
As in the previous section, in the classical limit the coherent states merely act to link the quantum numbers $n,l,m$ to their classical counterparts, while endowing the result with the correct fundamental frequencies. In other words, we have

\begin{eqnarray}\label{eq:greencoh2}
A^\mu_{\rm ret}(t,\vec{x})&=&\frac{iq}{ \mu}\,\sum_{l_{\gamma}=0}^\infty\sum_{m_{\gamma}=-l_{\gamma}}^{l_\gamma}\,\sum_{\Delta n,\Delta l,\Delta m}\,Y_{l_\gamma}^{m_\gamma*}(\theta,\varphi)\,\times\nonumber\\[5pt]
&&\int_{-\infty}^\infty\,\frac{\Theta(t-t')}{2\pi}\,dt'\,\int_{-\infty}^{\infty}\,d\omega\,\omega\,h^{(1)}_{l_{\gamma}}(\omega r)\,\exp\left[{-i\Delta n\,\alpha'(t')-i\omega(t-t') }\right]\,\,\mathcal{M}^{\mu}_{\Delta,l_\gamma,m_\gamma}(\omega,N,L)\,\nonumber\\[5pt]
\mathcal{M}^{\mu}_{\Delta,l_\gamma,m_\gamma}(\omega,N,L)&\equiv&\lim_{\hbar\rightarrow 0}\,\left\langle n'l'm'\right|\,j_{l_{\gamma}}(\omega\, r_{\rm{Kep}})Y_{l_\gamma}^{m_\gamma}(\theta_{\rm{Kep}},\varphi_{\rm{Kep}})\,p^\mu_{\rm{Kep}}\,\left|nlm\right\rangle\,,
\end{eqnarray}

where the matrix element is taken between $(n,l,m)=\hbar^{-1}(N,L,L)$ and $(n',l',m')=(n,l,m)-(\Delta n,\Delta l,\Delta m)$.

All that is left now is to calculate the classical matrix element. Its $0_{th}$ component is given by

\begin{eqnarray}\label{eq:Amumatel0}
&&\mathcal{M}^{0}_{\Delta,l_\gamma,m_\gamma}(\omega,N,L)=\mu\,\underbrace{\left\{\lim_{\hbar\rightarrow 0}\left\langle n' l'\right|\,j_{l_{\gamma}}(\omega r_{\rm{Kep}})\left| n l\right\rangle\right\}}_{{\rm Eq.}\,\text{\eqref{eq:rad1}}}~\underbrace{\left\{\lim_{\hbar\rightarrow 0}\left\langle l' m'\right| Y_{l_\gamma}^{m_\gamma}(\theta_{\rm{Kep}},\varphi_{\rm{Kep}})\left| l l \right\rangle\right\}}_{{\rm Eq.}\,\text{\eqref{eq:ang1}}}\,.
\end{eqnarray}

On the other hand, we have

\begin{eqnarray}\label{eq:Amumatelvec}
\overrightarrow{\mathcal{M}}_{\Delta,l_\gamma,m_\gamma}(\omega,N,L)&=&-i\underbrace{\left\{\lim_{\hbar\rightarrow 0}\left\langle n' l'\right|\,r^{-1}j_{l_{\gamma}}(\omega r_{\rm{Kep}})\left| n l\right\rangle\right\}}_{{\rm Eq.}\,\text{\eqref{eq:rad2}}}~\underbrace{\left\{\lim_{\hbar\rightarrow 0}\,\hbar \sum_{q=-1}^1\,\left\langle l' m'\right| Y_{l_\gamma}^{m_\gamma}(\theta_{\rm{Kep}},\varphi_{\rm{Kep}})(\overrightarrow{\nabla}_{\Omega})_q\left| l l \right\rangle\,\vec{\varepsilon}_q\,\right\}}_{{\rm Eq.}\,\text{\eqref{eq:ang2}}}\nonumber\\[5pt]
&&-i\underbrace{\left\{\lim_{\hbar\rightarrow 0}\hbar\left\langle n' l'\right|\,j_{l_{\gamma}}(\omega r_{\rm{Kep}})\,\partial_r\left| n l\right\rangle\right\}}_{{\rm Eq.}\,\text{\eqref{eq:rad3}}}~\underbrace{\left\{\lim_{\hbar\rightarrow 0}\sum_{q=-1}^1\,\left\langle l' m'\right| Y_{l_\gamma}^{m_\gamma}(\theta_{\rm{Kep}},\varphi_{\rm{Kep}})(\hat{r})_q\left| l l \right\rangle\,\vec{\varepsilon}_q\,\right\}}_{{\rm Eq.}\,\text{\eqref{eq:ang3}}}\,,
\end{eqnarray}

where $\vec{\varepsilon}_0=\hat{z},\,\vec{\varepsilon}_{\pm}=\tfrac{1}{\sqrt{2}}\left(\mp\hat{x}+i\hat{y}\right)$ and $(\vec{v})_q=\vec{v}\cdot\vec{\varepsilon}_{q}^{\,*}$. In the above equations, we added tags to indicate the relevant equations in the appendix where the corresponding classical matrix elements are calculated. Equation (<ref>) with the matrix elements (<ref>) and  (<ref>) is generic and exact, and we shall use it in the next section to calculate the waveform corresponding to an adiabatic inspiral, in a saddle point approximation.

For the rest of this section, we will focus on the special case of exact Keplerian motion. In this case $\alpha'=\Upsilon t'$ with $\Upsilon$ constant in time, and we can explicitly carry the $t'$ and $\omega$ integrals, obtaining

\begin{eqnarray}\label{eq:greencoh3}
&&A^\mu_{\rm ret}=\frac{iq}{\mu}\,\sum_{l_{\gamma}=0}^\infty\sum_{m_{\gamma}=-l_{\gamma}}^{l_\gamma}\,\sum_{\Delta n}\,\omega_{\Delta n}\,h^{(1)}_{l_{\gamma}}(\omega_{\Delta n} r)\,\,\exp\left[-i\Delta n\, \alpha\right]\,Y_{l_\gamma}^{m_\gamma*}(\theta,\varphi)\,\sum_{\Delta l,\Delta m}\, \mathcal{M}^{\mu}_{\Delta,l_\gamma,m_\gamma}(\omega_{\Delta n},N,L)\,,
\end{eqnarray}

where $\alpha=\Upsilon t$ and $ \omega_{\Delta n}\equiv \Upsilon \Delta n$. In a similar way, one can derive the advanced field for a classical electron in a Keplerian orbit:

\begin{eqnarray}\label{eq:greencoh3adv}
&&A^\mu_{\rm adv}=\frac{iq}{\mu}\,\sum_{l_{\gamma}=0}^\infty\sum_{m_{\gamma}=-l_{\gamma}}^{l_\gamma}\,\sum_{\Delta n}\,\omega_{\Delta n}\,\,h^{(2)}_{l_{\gamma}}(\omega_{\Delta n} r)\,\,\exp\left[-i\Delta n\, \alpha\right]\,Y_{l_\gamma}^{m_\gamma*}(\theta,\varphi)\,\sum_{\Delta l,\Delta m}\, \mathcal{M}^{\mu}_{\Delta,l_\gamma,m_\gamma}(\omega_{\Delta n},N,L)\,.
\end{eqnarray}

The only difference from (<ref>) is the substitution $h^{(1)}_{l_\gamma}\rightarrow h^{(2)}_{l_\gamma}$. For future reference, we can also write down the regular field

\begin{eqnarray}\label{eq:Areg}
A^\mu_{\rm reg}&=&\frac{1}{2}\left(A^\mu_{\rm ret}-A^\mu_{\rm adv}\right)=\nonumber\\[5pt]
&=&\frac{iq}{\mu}\,\sum_{l_{\gamma}=0}^\infty\sum_{m_{\gamma}=-l_{\gamma}}^{l_\gamma}\,\sum_{\Delta n}\,\omega_{\Delta n}\,\,j_{l_{\gamma}}(\omega_{\Delta n} r)\,\,\exp\left[-i\Delta n\, \alpha\right]\,Y_{l_\gamma}^{m_\gamma*}(\theta,\varphi)\,\sum_{\Delta l,\Delta m}\, \mathcal{M}^{\mu}_{\Delta,l_\gamma,m_\gamma}(\omega_{\Delta n},N,L)\,.
\end{eqnarray}

This is the field that is responsible for the EM self-force (see <cit.> for a review), as we shall see in the next section.

In Fig. <ref>, we present waveforms –$A_t$–  corresponding to various Keplerian orbits. Specifically, we consider the four combinations of circular ($e=0$) versus elliptic ($e=0.5$), and "fast" ($p\mu/K=20$) versus "slow" ($p\mu/K=10^4$) orbits. Note that the asymmetry between the maximum and minimum of the waveforms is due to the observation point being on the x-axis. The qualitative features of the waveforms can be understood as follows. For circular orbits, it can be shown that the Fourier components of (<ref>) vanish unless $\Delta n= m_{\gamma}$ can be satisfied, whereas the Fourier components of elliptic curves are generally non-vanishing. Moreover, unlike "fast" orbits, "slow" orbits are well-described using the dipole-approximation ($l_\gamma=1$), as can be seen in Fig. <ref>. Hence, "slow" circular orbits have sinusoidal waveforms, in contrast to more involved waveforms in the general case. Finally, the broadening (narrowing) of the minimum (maximum) of the waveforms is due to the Doppler effect.

Keplerian Waveforms. This figure shows the  waveform –$A_t$– radiated over one period by an electron undergoing Keplerian motion in four cases of the orbital parameters $(e,\Bar{p})$: (i) $(0,20)$; (ii) $(0,10^4)$; (iii) $(0.5,20)$; and (iv) $(0.5,10^4)$, where $\Bar{p}=p\mu/K$. The observation point is on the x-axis, far away from the electron's orbit. For each case, the horizontal and vertical axes are normalized by the orbital period $T=2\pi/\Upsilon$ and the maximum of the waveform, respectively.

The multipole contributions to the Fourier coefficients of (<ref>), for "fast" and "slow" elliptic orbits, are presented in Fig. <ref>. As mentioned above, in the "fast" orbit case (Fig. <ref>), the Fourier coefficients receive substantial contributions from higher multipoles, whereas, in the "slow" orbit case (Fig. <ref>), the dipole is by far the dominant contribution.

As can be seen from the relative error plots in Fig. <ref>, the Fourier coefficients of  (<ref>) are indeed equal to the standard coefficients obtained from classical electromagnetism, expressed using (numerical) Bessel integrals.

Multipole contributions and relative error for $(e,\Bar{p})=(0.5,20)$. We show the first 10 leading multipole contributions,  ignoring the monopole $(l_\gamma=0)$ contribution.

Multipole contributions and relative error for $(e,\Bar{p})=(0.5,10^4)$. We show the first 3 leading multipole contributions,  ignoring the monopole $(l_\gamma=0)$ contribution.

Multipole Contributions and Equality with the Classical Integrals. This figure shows the contributions of the dominant multipoles to the Fourier coefficients of $A_t$ in (<ref>), where we exclude $h^{(1)}_{l_{\gamma}} $ for aesthetic purposes. We take two sets of orbital elements $(e,\Bar{p})$, where $\Bar{p}=p\mu/K$. The vertical axis is normalized by the maximum in each panel. The relative error of the Fourier coefficients with respect to the (numerical) classical  integrals is shown in the bottom inserts.

§ EM SELF-FORCE WITH THE QSM

§.§ Adiabatic EM Inspirals

The natural application of the QSM is in the calculation of the self-force on an inspiralling body. Here and below, we calculate the self-force for the case of an inspiralling classical electron. Before delving into the details of the self-force itself, we first briefly describe how it is used to calculate inspiralling trajectories in PA perturbation theory. Consider the EOM for a classical electron in a central potential, with an additional small self-force $\vec{F}$:

\begin{eqnarray}\label{eq:EOMF}
\mu\ddot{\vecc{r}}+K\frac{\vecc{r}}{r^3}=\varepsilon \vec{F}\,.
\end{eqnarray}

The $\varepsilon$ in front of the force $\vec{F}$ is a power-counting parameter. Expanding in it gives us PA perturbation theory. The self-force itself may also depend on $\varepsilon$, both explicitly and through its dependence on the trajectory of the inspiralling electron.
For the systematic solution of (<ref>), order-by-order in $\varepsilon$, we first have to choose a convenient parameterization for the inspiralling trajectory. The standard choice is using osculating orbits <cit.> - a parameterization of the true trajectory of an inspiralling object as transitioning between different Keplerian orbits that are momentarily tangent (”osculating") to the true trajectory. At time $t$, the corresponding osculating orbit is given in terms of the orbital elements $(e(t),\,p(t))$ and the mean anomaly $\alpha(t)$ via (<ref>), with $\psi(t)$ related to $\alpha(t)$ using  (<ref>). Unlike in pure Keplerian motion, the orbital elements $(e(t),p(t))$ (equivalently $E(t)$ and $L(t)$) are all taken to be functions of time. The time dependence of $e(t),\,p(t)$, and $\alpha(t)$ is derived from the osculation (tangentiality) conditions between the osculating and true orbits, namely:

\begin{eqnarray}\label{eq:osc}
\vec{r}_{\rm Kep}\left[p(t),e(t),\alpha(t)\right]&=&\vec{r}\nonumber\\[5pt]
\dot{\vec{r}}_{\rm Kep}\left[p(t),e(t),\alpha(t)\right]&=&\dot{\vec{r}}\,,
\end{eqnarray}

where $\vec{r}$ satisfies the full EOM (<ref>), which in turn depends on the self-force $\vec{F}$. These relations greatly simplify in our case of interest, in which $\vec{F}$ is the force due to electromagnetic radiation-reaction. This force is purely dissipative, and so it is enough to know $dE/dt$ and $dL/dt$ to determine the complete time evolution of the system. In particular, the time evolution of $(p(t),e(t))$ in this case is fixed to be

\begin{eqnarray} \label{0PAUp}
\end{eqnarray}

where $dE/dt$ and $dL/dt$ are explicitly calculated below, at the adiabatic (0PA) order. We still need to solve for the time evolution of $\alpha(t)$ from the osculation conditions (<ref>), but this would become trivial below once we specialize to the adiabatic order.

Importantly, the method of osculating orbits is not an approximation but rather a parameterization of the full trajectory. This method becomes particularly useful as long as the energy and angular momentum loss rates are slow, so that $\frac{T\dot{E}}{E}\sim \frac{T\dot{L}}{L}\ll 1 $, where $T$ is the rotation period - this is the case for us when $Z\gg 1$. Within this framework, we can apply PA perturbation theory <cit.> in the power counting parameter $\varepsilon$. In this work, we limit ourselves to the leading order in $\varepsilon$, also called the adiabatic (0PA) order. At this order, the dissipation rates of $\dot{E}$ and $\dot{L}$ are calculated for the osculating Keplerian orbit, without back-reaction.

The self-force acting on a classical electron is nothing but the Lorentz force from its own EM (self)-field. However, care must be taken to use a correctly regulated self-field, which does not diverge at the position of the electron itself. By Dirac's famous prescription <cit.> (see <cit.> for a curved space generalization), this regulated field is given by $A^\mu_{\rm reg}=\frac{1}{2}\left(A^\mu_{\rm ret}-A^\nu_{\rm adv}\right)$, rather than the retarded field $A^\mu_{\rm ret}$ used in the previous section. The regular field $A^\mu_{\rm reg}$ has to be evaluated at the momentary position of the particle, which at 0PA is taken to be the osculating orbit itself, without back-reaction. In other words, for every value of the orbital elements $(e,p)$ and the mean anomaly $\alpha$, we can find the corresponding point $(t,\vec{x})$ in space-time where we evaluate $A^\mu_{\rm reg}$. This essentially defines the function $A^\mu_{\rm reg}(e,p,\alpha)$, the regular field at a point on the osculating orbit.

Finally, $A^\mu_{\rm reg}(e,p,\alpha)$ is used to compute the energy and angular momentum loss rates via

\begin{eqnarray}\label{eq:dEdL}
\frac{dE}{dt}&=&q\,\vec{v}\cdot \vec{E}^{\rm reg}\nonumber\\[5pt]
\frac{d\vec{L}}{dt}&=&q\,\vec{r}\,\times\,\left[\vec{E}^{\rm reg}+\vec{v}\times \vec{B}^{\rm reg}\right]\,,
\end{eqnarray}

where $E^{\rm reg}_i=\partial_{[0}A^{\rm reg}_{i]}$ and $B^{\rm reg}_i=-\frac{1}{2}\epsilon_{ijk}\partial^{[j}A^{\rm reg;k]}$, evaluated at $(e(t),p(t),\alpha(t))$.
Equations (<ref>) are supplemented by an equation for $d\alpha/dt$, which is derived from the osculation condition (<ref>). Together, these three equations constitute a full set of differential equations that determine the time evolution of the inspiral. However, the right-hand side of (<ref>) is a highly nonlinear function of $\alpha$, which in practice makes the system challenging for numerical integration. Reference <cit.> addressed this problem by performing a clever change of variables called a near-identity transformation (NIT). At 0PA and under a particular choice of the NIT[This is the choice made in section 2.6.4 of <cit.> and section 6.2.5 of <cit.>], the evolution equations become

\begin{eqnarray} \label{0PA}
&&{\left(\frac{dE}{dt}\right)}_{0\text{PA,NIT}}= \left\langle\frac{d E}{dt}\right\rangle_{\alpha}~~~,~~~{\left(\frac{dL}{dt}\right)}_{0\text{PA,NIT}}=\left\langle\frac{d L}{dt}\right\rangle_{\alpha}~~~,~~~{\left(\frac{d\alpha}{dt}\right)}_{0\text{PA,NIT}}=\Upsilon(t)\equiv\sqrt{\frac{K\left(1-e^2(t)\right)^3}{\mu p^3(t)}}\,.\nonumber\\
\end{eqnarray}

where $\left\langle\mathcal{O}(\alpha)\right\rangle_{\alpha}\equiv\frac{1}{2\pi}\int^{2\pi}_0\,\mathcal{O}\,d\alpha$, and so the $\alpha$ dependence is eliminated from the right-hand side.
To calculate the right-hand side of (<ref>), all we have to do is calculate $A^\mu_{\rm reg}$ generated by an electron in an osculating Keplerian orbit, and evaluated on the very same orbit (averaged over $\alpha$). This is done via the QSM in a manner similar to the previous section. Remarkably, it turns out that these expressions are no other than the quantum rates for spontaneous emission.

§.§ 0PA Radiation as the Classical Limit of Spontaneous Emission

We are now ready to compute $dE/dt$ and $dL/dt$ due to the EM self-force on the classical electron. Here we immediately see the power of the method of osculating orbits. In this framework, we only need to compute $dE/dt$ and $dL/dt$ for a source moving not on the full trajectory, but only on the osculating orbit[At 0PA, we only need to compute $dE/dt$ and $dL/dt$ for the osculating orbit. At higher PA orders, we also need the partial derivatives of $dE/dt$ and $dL/dt$ with respect to $p,e,\alpha$, evaluated again on the osculating orbit.]. Consequently, we can directly substitute (<ref>) in (<ref>) and get

\begin{eqnarray}\label{eq:dE2}
&&\left\langle\frac{dE}{dt}\right\rangle_{\alpha}=\frac{iq^2}{2\pi \mu^2}\,\,\lim_{\hbar\rightarrow 0}\,\sum_{\Delta n,\Delta l,\Delta m}\,\omega_{\Delta n}\,\int_{0}^{2\pi}\,d\alpha\,\exp\left[-i\Delta n \,\alpha\right]\,\sum_{l_{\gamma}=0}^\infty\sum_{m_{\gamma}=-l_{\gamma}}^{l_\gamma}\,\times\nonumber\\[5pt]
&&~~~~~~~\left\{i\omega_{\Delta n}\,\,j_{l_{\gamma}}(\omega_{\Delta n} r_{\rm{Kep}})\,Y_{l_\gamma}^{m_\gamma*}(\theta_{\rm{Kep}},\varphi_{\rm{Kep}})\,p^i_{\rm{Kep}}\,\mathcal{M}_{i}-\,\partial_i\left[j_{l_{\gamma}}(\omega_{\Delta n} r_{\rm{Kep}})\,Y_{l_\gamma}^{m_\gamma*}(\theta_{\rm{Kep}},\varphi_{\rm{Kep}})\right]\,p^i_{\rm{Kep}}\,\mathcal{M}_{0}\right\}\,,\nonumber\\
\end{eqnarray}

where $\left(r_{\rm{Kep}},\theta_{\rm{Kep}},\varphi_{\rm{Kep}}\right)$ are evaluated at $\alpha$. Using the QSM for the $t$-dependent terms (this time with a coherent state defined with respect to the time $t$ rather than $t'$), we get

\begin{eqnarray}\label{eq:dE3}
&&\left\langle\frac{dE}{dt}\right\rangle_{\alpha}=\frac{iq^2}{2\pi \mu^2}\,\lim_{\hbar\rightarrow 0}\,\sum_{\Delta n,\Delta l,\Delta m}\,\omega_{\Delta n}\,\int_{0}^{2\pi}\,d\alpha\,\exp\left[-i\Delta n \,\alpha\right]\,\sum_{l_{\gamma}=0}^\infty\sum_{m_{\gamma}=-l_{\gamma}}^{l_\gamma}\,\times\nonumber\\[5pt]
&&~~~~~~~\left\{i\omega_{\Delta n}\,\left\langle{\rm coh}\right|j_{l_{\gamma}}\,Y_{l_\gamma}^{m_\gamma }\,p^i\left|{\rm coh}\right\rangle ^{*}\,\mathcal{M}_{i}-\,\left\langle{\rm coh}\right|\partial_i\left[j_{l_{\gamma}}\,Y_{l_\gamma}^{m_\gamma}\right]\,p^i\left|{\rm coh}\right\rangle ^{*} \,\mathcal{M}_{0}\right\}\,,
\end{eqnarray}

where we suppressed the arguments of $j_{l_\gamma}$ and $Y_{l_\gamma}^{m_\gamma}$ for brevity. In the classical limit we have

\begin{eqnarray}\label{eq:dE3t}
&&\left\langle\frac{dE}{dt}\right\rangle_{\alpha}=\frac{iq^2}{2\pi \mu^2}\,\lim_{\hbar\rightarrow 0}\sum_{\Delta n,\Delta l,\Delta m}\sum_{\Delta \tilde{n}}\,\omega_{\Delta n}\,\int_{0}^{2\pi}\,d\alpha\,\exp\left[-i(\Delta n-\Delta \tilde{n}) \alpha\right]\,\sum_{l_{\gamma}=0}^\infty\sum_{m_{\gamma}=-l_{\gamma}}^{l_\gamma}\,\times\nonumber\\[5pt]
&&~~~~~~~\left\{i\omega_{\Delta n}\,\left\langle \tilde{n}'l'm'\right|j_{l_{\gamma}}\,Y_{l_\gamma}^{m_\gamma}\,p^i\left|nlm\right\rangle ^{*}\,\mathcal{M}_{i}-\,\left\langle\tilde{n}'l'm'\right|\partial_i\left[j_{l_{\gamma}}\,Y_{l_\gamma}^{m_\gamma }\right]\,p^i\left|nlm\right\rangle ^{*}\,\mathcal{M}_{0}\right\}\,,
\end{eqnarray}

where $\tilde{n}'=n-\Delta\tilde{n}$. Note that in the classical limit, the matrix element is nonzero only for $l'=m'=m-m_\gamma + \left\{0,\pm 1\right\}$, so $\Delta l,\,\Delta m$ are fixed and there is no need to introduce $\Delta \tilde{l},\,\Delta\tilde{m}$.
The $\alpha$ integration gives $2\pi\delta_{\Delta n,\Delta \tilde{n}}$, and so

\begin{eqnarray}\label{eq:dE4}
\left\langle\frac{dE}{dt}\right\rangle_{\alpha}&=&\frac{iq^2}{\mu^2}\,\lim_{\hbar\rightarrow 0}\sum_{\Delta n,\Delta l,\Delta m}\omega_{\Delta n}\,\sum_{l_{\gamma}=0}^\infty\sum_{m_{\gamma}=-l_{\gamma}}^{l_\gamma}\,\left\{ \,i\omega_{\Delta n}\,\left|\left\langle n' l' m'\right| j_{l_{\gamma}}(\omega_{\Delta n} r)\,Y_{l_\gamma}^{m_\gamma}\,p^i\left|nlm\right\rangle\right|^2-\right.\nonumber\\[5pt]
&&~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\left.\mu\,\left\langle  n' l' m'\right|\partial_i\left[j_{l_{\gamma}}(\omega_{\Delta n} r)\,Y_{l_\gamma}^{m_\gamma}\right]\,p^i\left|nlm\right\rangle^*\,\left\langle n'l'm'\right|j_{l_{\gamma}}(\omega_{\Delta n} r)\,Y_{l_\gamma}^{m_\gamma}\left|nlm\right\rangle\right\}\,.\nonumber\\
\end{eqnarray}

Finally, using (<ref>) in Appendix <ref>, we can cast this expression in the suggestive form

\begin{eqnarray} \label{SE1main}
\left\langle\frac{dE}{dt}\right\rangle_{\alpha}=-\lim_{\hbar\rightarrow 0}\sum_{\Delta n>0,\Delta l,\Delta m}\left(E_{n}-E_{n'}\right)\,\Gamma_{s.e.}
\end{eqnarray}

where $\Gamma_{s.e.}$ is the transition rate given by

\begin{eqnarray} \label{SE3main}
\Gamma_{s.e.}=\frac{q^{2}\left(E_{n}-E_{n'}\right)}{8\pi^2\hbar^2\mu^{2}}\int d^{2}\Omega_{k}\,\left.\left\{ \sum_{\sigma}\left|\left\langle n',l',m'\right|e^{-i\vec{k}\cdot\vec{x}}\left(\vec{\varepsilon}^{~\sigma}_{\vec{k}}\cdot\vec{p}\right)\left|n,l,m\right\rangle \right|^{2}\right\} \right|_{k=\left(E_{n}-E_{n'}\right)/\hbar}\,.
\end{eqnarray}

This is simply the expression for spontaneous emission, as calculated via Fermi's golden rule. In a similar manner, one can check that

\begin{eqnarray} \label{SE2main}
\left\langle\frac{dL}{dt}\right\rangle_{\alpha}=-\lim_{\hbar\rightarrow 0}\,\sum_{\Delta n>0,\Delta l,\Delta m}\,\hbar\,\left(l-l'\right)\,\Gamma_{s.e.}\,.
\end{eqnarray}

In practice, when calculating the all-multipole EM self-force, we make use of the explicit expression

\begin{eqnarray} \label{SE4main}
\Gamma_{s.e.}= -\frac{2q^{2}}{\hbar \mu^2}\,\Upsilon \Delta n\,\sum_{l_{\gamma}=0}^\infty\sum_{m_{\gamma}=-l_{\gamma}}^{l_\gamma}\,\mathcal{M}^*_\mu\mathcal{M}^\mu\,+\,\mathcal{O}(\hbar^0)\,,
\end{eqnarray}

which can be derived either from (<ref>) or directly from (<ref>).

Full Versus Dipole Inspiral. This figure shows the adiabatic (0PA) trajectory (green) of the electron obtained from (<ref>) - (<ref>). The dipole trajectory (blue), obtained by retaining only the dipole-order contributions in (<ref>), is also shown. The initial values of the orbital parameters are $(e,\Bar{p})=(0.5,300)$. Moreover, in this figure alone, we took an exaggerated value of $Z$ to make the evolution more discernible. We note that both trajectories start together from the same position, and run for the same duration. As seen in this figure, the full trajectory spirals deeper than the dipole one.

§.§ The 0PA Trajectories and Waveforms

To obtain the full 0PA trajectory, we can now integrate (<ref>) numerically using the expressions (<ref>) and (<ref>) for the period-averaged energy and angular momentum loss rates, and get $(e(t),p(t),\alpha(t))$. The distance from the center and azimuthal angle are then given by

\begin{eqnarray} \label{0PAdist}
r(t)&=&\frac{p(t)}{1-e^2(t)}\left\{1+\frac{e^2(t)}{2}-2e(t)\sum_{\Delta n=1}^\infty\,\frac{1}{\left(\Delta n
\right)^2}\frac{d J_{\Delta n}\left(\Delta n \,e(t)\right)}{de}\cos\left(\Delta n\, \alpha(t)\right)\right\}\nonumber\\[5pt]
\varphi(t)&=&2\arctan\left\{\sqrt{\frac{1+e(t)}{1-e(t)}}\tan\left(\beta(t)/2\right)\right\}\nonumber\\[5pt]
\beta(t)&=&\alpha(t)+2\sum_{\Delta n=1}^\infty\,\frac{1}{\Delta n
}J_{\Delta n}\left(\Delta n \,e(t)\right)\,\sin\left(\Delta n\, \alpha(t)\right)\,.
\end{eqnarray}

Once the inspiralling trajectory is known, the associated waveform can be extracted using the expression (<ref>). Note that this time we cannot use the simplified expression (<ref>), which is only valid for a Keplerian source. In practice, (<ref>) is cumbersome to evaluate, and so, as is the case for gravitational inspirals, one usually resorts to approximate methods for waveform generation. To approximate (<ref>), we first focus on the far-field regime, where $h^{(1)}_{l_\gamma}(\omega r)$ attains its asymptotic value $(-i)^{l_\gamma+1} e^{i\omega r}/(\omega r)$. Furthermore, in the far field regime, all time and distance scales are hierarchically larger than the inspiral period, and so we can use the stationary phase (saddle point) approximation for the exponential in (<ref>). The stationary phase condition over $
\omega$ and $t'$ gives

\begin{eqnarray}
t'=t-r~~~~\omega=\Delta n \dot{\alpha}'=\Delta n \Upsilon(t')\,.
\end{eqnarray}

Finally, we have

\begin{eqnarray}\label{eq:greencoh3'}
&&A^\mu_{\rm ret,stationary}=\frac{q}{\mu r}\,\sum_{l_{\gamma}=0}^\infty\sum_{m_{\gamma}=-l_{\gamma}}^{l_\gamma}\,\sum_{\Delta n}\,\exp\left[-i\Delta n\, \alpha(t_{ret})\right]\,(-i)^{l_\gamma}\,Y_{l_\gamma}^{m_\gamma*}(\theta,\varphi)\,\sum_{\Delta l,\Delta m}\, \mathcal{M}^{\mu}_{\Delta,l_\gamma,m_\gamma}(\Delta n \Upsilon(t_{ret}),N,L)\,,\nonumber\\
\end{eqnarray}

where $t_{ret}=t-r$. This far-field approximation to the emitted EM field is very easy to evaluate numerically.
A sample inspiralling trajectory is presented in Fig. <ref> (with an exaggeratedly small value of $Z$ to make the evolution more discernible), and a typical waveform is shown in Fig. <ref>. The corresponding dipole-order inspiral and waveforms, obtained by retaining only the dipole ($l_\gamma=1$) contributions in (<ref>), are also shown. At the beginning of the inspiral, the dipole-approximation is valid, and as such the trajectories (Fig.  <ref>) and waveform (Fig.  <ref>) match. As the inspiral continues, the dipole-approximation breaks down, and a phase mismatch develops (Fig.  <ref>). At later times, a frequency and amplitude mismatch also develops (Fig.  <ref>). This outlines the importance of higher multipoles for the accurate description of inspiral dynamics.

The (normalized) waveform produced during the first 5 periods in the inspiral. The full (green) waveform does not appear because it coincides with the dipole one in this duration.

The (normalized) waveform produced from $t=1100\,
T_{\rm{init}}$ to $t=1105\, T_{\rm{init}}$.

The (normalized) waveform produced from $t=2100\,
T_{\rm{init}}$ to $t=2105\, T_{\rm{init}}$.

Full Versus Dipole Waveform. This figure shows (in green) the adiabatic (0PA) waveform – $A_{t}$ –  of the electron obtained from (<ref>) - (<ref>) for $Z=4\pi$. The dipole waveform (blue), obtained by retaining only the dipole-order contributions in (<ref>), is also shown. The observation point is on the x-axis, far away from the electron's inspiral.
The initial orbital parameters of the electron's trajectory are given by $(e_0,\Bar{p}_0)=(0.5,300)$.
The horizontal axis denotes the retarded time, $t-r$, normalized by the period of the initial Keplerian orbit, $T_{\rm{init}}$. The vertical axis is normalized by the maximum of the initial Keplerian waveform.

§ CONCLUSIONS AND OUTLOOK

In this paper we presented the Quantum Spectral Method (QSM) for the calculation of classical Fourier coefficients from their corresponding quantum matrix elements. The QSM allows for the exact analytical calculation of classical quantities that so far have only been calculated numerically.

We demonstrated the QSM for the case of a classical electron moving in a $1/r$ potential. By taking the classical limit of hydrogen-like atom matrix elements, we were able to analytically compute to all multipole orders: (a) time dependent Keplerian motion, (b) the electromagnetic (EM) field radiated by a classical electron in a (quasi)-Keplerian orbit , and; (c) an adiabatic EM inspiral and its associated waveform. Anecdotally, in the cases studied, our analytical expression allowed for a faster evaluation in Mathematica than their traditional numerical counterparts. We leave the systematic study of the potential runtime improvement of self-force problems for future work.

More generally, in the QSM, the classical Fourier coefficients emerge as the classical limit of quantum matrix elements, obtained using the full quantum description of bound states. In this way, the QSM is inherently non-perturbative in the coupling constant, and so it is naturally suitable for self-force/post-adiabatic (PA) calculations. This is in contrast to other quantum-to-classical approaches which start from scattering amplitudes, and are thus inherently perturbative in the coupling constant. The perturbativity of these quantum-to-classical approaches makes them well suited for post-Minkowskian (PM) <cit.> and post-Newtonian (PN) <cit.> calculations, but not for the PA approach <cit.>, which is perturbative in the mass ratio, but not in the coupling. It would be very interesting to compare between our approach and the perturbative quantum-to-classical approach, in the spirit of <cit.>. For example, the results of this work could be compared to the perturbative EM calculations in <cit.>, and the future gravitational application of the QSM should be compared to a resummation of self-force Effective Field Theory <cit.>.

Though the current demonstration of the QSM is within the well-studied case of a classical electron in a $1/r$ potential, the validity and usefulness of the method could go far beyond this scope. In particular, it should be straightforward to generalize our method to systems with multiple fundamental frequencies, and to higher orders in the PA expansion. Moreover, the QSM should be equally applicable for non-periodic systems, with the Fourier series replaced by a Fourier transform. Finally, in our future work, we intend to apply the QSM to  EM and gravitational inspirals in Schwarzschild and Kerr backgrounds. In this case, the eigenstates of the single-particle Hamiltonian become solutions of the Regge-Wheeler and Teukolsky equations <cit.>, given by Heun functions (see, e.g. <cit.>). Importantly, the explicit computation of Kerr quasinormal modes is not required for the application of the QSM, since they go over to integer multiples of the classical fundamental frequencies for bound Schwarzschild/Kerr geodesics. There is also no need to explicitly construct coherent states for BH backgrounds, as one could directly use the QSM master formula (<ref>). At the adiabatic level (0PA), the relevant first order self-force can be computed using the Green's function method, analogous to Section <ref>. At the first post-adiabatic order (1PA), the second-order self-force involves an effective gravitational source <cit.>, nevertheless, we expect the QSM to be applicable in conjunction with metric reconstruction schemes <cit.>. The recoil of the heavy black hole, discussed in <cit.>, would be accounted for in a similar way to the standard PA treatment <cit.>. Furthermore, the QSM is well suited for the incorporation of classical spin into the PA expansion <cit.>, by taking the classical limit of quantum spin-orbit coupling. This would serve as a cross-check to the existing results for spinning bodies in the PN and PM expansions <cit.>. Overall, with these advantages in mind, we find it worthwhile to consider incorporating the QSM as a module in self-force/gravitational inspiral algorithms.

§ CLASSICAL LIMIT OF COHERENT STATES

§.§ Hydrogen Atom Coherent States

We take a slightly modified version of Crawford's coherent states <cit.>. In this section we express them more systematically as:

\begin{eqnarray}\label{eq:cohapp1}
\left|{\rm coh}\right\rangle=\left|N,\eta,\alpha\right\rangle=&&\mathcal{A}(N/\hbar)\sum_{n}\,g_n(N,\alpha)\,\delta_{2j+1,n}\,\left|j,\eta\right\rangle\nonumber\\[5pt]
g_n(N,\alpha)\equiv &&\,\frac{n\,(N/\hbar)^{n-1}\exp\left[{-i\frac{\alpha E_{n}}{\Upsilon\hbar}}\right]}{\sqrt{(n-1)!}}\,,
\end{eqnarray}

with the normalization $\mathcal{A}(N/\hbar)=e^{-N/2\hbar}(1+3\hbar^{-1}N+\hbar^{-2}N^2)^{-1/2}$.
Here $\left|j,\eta\right\rangle$ is a (XY planar) coherent state realizing the $SO(4)$ symmetry of the hydrogen atom <cit.>. It is defined as

\begin{eqnarray}\label{eq:cohapp2}
\left|j,\eta\right\rangle&=&\sum_{m_1,m_2=-j}^j\,a^j_{m_1,m_2}\left(\eta\right)\left|j,m_1;j,m_2\right\rangle\,\nonumber\\[5pt]
\end{eqnarray}

The $\left|j,m_1;j,m_2\right\rangle$ are eigenstates of $J^2_{\pm},\,J^z_{\pm}$, where $\vec{J}_{\pm}=\frac{1}{2}\left(\vec{L}\pm \vec{A}\right)$. Here $\vec{L}$ is the angular momentum operator and $\vec{A}$ is the Runge-Lenz operator. They are related to the usual hydrogen atom eigenstates via Clebsch-Gordan coefficients,

\begin{eqnarray}
\left|j,m_1;j,m_2\right\rangle=\sum_{l,m}\,C^{lm}_{jm_1;j,m_2}\,\left|n,l,m\right\rangle\,.
\end{eqnarray}

Here again $2j+1\equiv n$.

§.§ Classical Limit

To take the classical limit of an expectation value, we take a step-by-step approach, first using (<ref>) and then unpacking the $SO(4)$ coherent state. We have

\begin{eqnarray}
\lim_{\hbar\rightarrow 0}\,\left\langle {\rm coh}\right|\mathcal{O}\left|{\rm coh}\right\rangle=\sum_{\Delta n}\,\lim_{\hbar\rightarrow 0}\,\sum_{n}\,g^*_{n'}\left(N,\alpha\right)\,g_{n}\left(N,\alpha\right)\,\left\langle \frac{n'-1}{2},\eta\right|\mathcal{O}\left|\frac{n-1}{2},\eta\right\rangle\,,
\end{eqnarray}

where $n'=n-\Delta n$. We can use Poisson summation for the $n$ sum,

\begin{eqnarray}
\lim_{\hbar\rightarrow 0}\,\left\langle {\rm coh}\right|\mathcal{O}\left|{\rm coh}\right\rangle&=&\sum_{\Delta n}\,\lim_{\hbar\rightarrow 0}\,\sum_{\tilde{n}}\,\int\,dn\,\mathcal{A}^2(N/\hbar)\,{\left(\frac{N}{\hbar}\right)}^{2n-\Delta n-1}\,\frac{n(n-\Delta n)\exp\left[{-i\frac{\alpha (E_{n}-E_{n-\Delta n})}{\Upsilon\hbar}}\right]}{\sqrt{(n-1)!(n-\Delta n-1)!}}\,\exp\left(-2\pi i n\tilde{n}\right)\times\nonumber\\[5pt]
&&\left\langle \frac{n-\Delta n-1}{2},\eta\right|\mathcal{O}\left|\frac{n-1}{2},\eta\right\rangle\,.
\end{eqnarray}

In the $\hbar\rightarrow 0$ limit, the Laplace (saddle point) method for the $n$ integral becomes exact. The integral only receives contributions from $n=\hbar^{-1}N,\,\tilde{n}=0$, and the prefactor neatly becomes $1$, up to a time-dependent phase. We then have

\begin{eqnarray}
\lim_{\hbar\rightarrow 0}\,\left\langle {\rm coh}\right|\mathcal{O}\left|{\rm coh}\right\rangle&=&\sum_{\Delta n}\,\exp\left[-i\Delta n \alpha\right]\,\lim_{\hbar\rightarrow 0}\,\left\langle j',\eta\right|\mathcal{O}\left|j,\eta\right\rangle\,,
\end{eqnarray}

where we fixed $j=\frac{\hbar^{-1} N-1}{2}$ and $j'=j-\Delta n/2$. Next, we take the classical limit of $\left\langle j',\eta\right|\mathcal{O}\left|j,\eta\right\rangle$, using (<ref>). Again doing Poisson summation on the $m_1,\,m_2$ sums, and taking the saddle point, we have

\begin{eqnarray}
\lim_{\hbar\rightarrow 0}\,\left\langle j',\eta\right|\mathcal{O}\left|j,\eta\right\rangle&=&\sum_{\Delta m_1}\,\sum_{\Delta m_2}\,\lim_{\hbar\rightarrow 0}\,\left\langle j',m'_1;j',m'_2\right|\mathcal{O}\left|j,m_1;j,m_2\right\rangle\,,
\end{eqnarray}

where $m_1=m_2=\frac{N}{2\hbar}\frac{1-\eta^2}{1+\eta^2}=\frac{L}{2\hbar}$ and $m'_{1,2}=m_{1,2}-\Delta m_{1,2}$. Finally, we take the classical limit of the Clebsch-Gordan and get

\begin{eqnarray}
\lim_{\hbar\rightarrow 0}\,\left\langle j',m';j',m'\right|\mathcal{O}\left|j,m;j,m\right\rangle&=& \lim_{\hbar\rightarrow 0}\,\sum_{\Delta l}\,\left\langle n',l',m'\right|\mathcal{O}\left|n,l,m\right\rangle\,,
\end{eqnarray}

where $l=m=\hbar^{-1}L$, $l'=l-\Delta l$ and $m'=m-\Delta m_1-\Delta m_2$. Putting everything together, we have

\begin{eqnarray}\label{eq:derclas}
\lim_{\hbar\rightarrow 0}\,\left\langle {\rm coh}\right|\mathcal{O}\left|{\rm coh}\right\rangle&=&\sum_{\Delta n}\,\exp\left[-i\Delta n \alpha\right]\,\sum_{\Delta l,\Delta m}\,\lim_{\hbar\rightarrow 0}\,\left\langle n,l',m'\right|\mathcal{O}\left|n,l,m\right\rangle\,,
\end{eqnarray}

where again $(n,l,m)=\hbar^{-1}(N,L,L)$ and $(n',l',m')=(n,l,m)-(\Delta n,\Delta l,\Delta m)$. Here we combined the sums over $\Delta m_{1,2}$ to a sum over $\Delta m$. A simple check of our overall normalization is done by setting $\mathcal{O}=1$.

§ THE QUANTUM MATRIX ELEMENTS OF THE HYDROGEN ATOM

§.§ Angular Matrix Elements

§.§.§ Spherical Harmonic

The angular matrix element $\left\langle l',m'\right|Y_{l_\gamma}^{m_\gamma}\left(\hat{r}\right)\left|l,l\right\rangle$ can be evaluated using the Wigner
$3j$ symbols, yielding

\begin{equation} \label{SphericalScalarElem}
\left\langle l',m'\right|Y_{l_\gamma}^{m_\gamma}\left(\hat{r}\right)\left|l,l\right\rangle=\left(-1\right)^{m'}\delta_{m',l+m_\gamma}\sqrt{\frac{\left(2l+1\right)\left(2l_\gamma+1\right)\left(2l'+1\right)}{4\pi}}
\left(\begin{array}{ccc}
l' & l & l_\gamma\\
0 & 0 & 0
\end{array}\right)\left(\begin{array}{ccc}
l' & l & l_\gamma\\
-\left(l+m_\gamma\right) & l & m_\gamma
\end{array}\right).
\end{equation}

§.§.§ Spherical Harmonic times $\hat{r}$

Vectorial matrix elements are easily computed using a tensor operator basis, namely, an arbitrary vector $\vec{v}$ is represented as

\begin{eqnarray}
\vec{v}=\sum_{q=-1}^1\,(\vec{v})_q\,\vec{\varepsilon}^{\,q}\,\,,
\end{eqnarray}

where $\vec{\varepsilon}_{0}=\hat{z},\,\vec{\varepsilon}_{\pm}=\tfrac{1}{\sqrt{2}}\left(\mp\hat{x}+i\hat{y}\right)$ and $(\vec{v})_q=\vec{v}\cdot\vec{\varepsilon}^{\,q*}$. We start with the matrix element $\left\langle l',m'\right|Y_{l_\gamma}^{m_\gamma}\left(\hat{r}\right)\,\left(\hat{r}\right)_{q}\left|l,l\right\rangle $.
We can calculate this by adding first the angular momentum of $Y_{l_\gamma}^{m_\gamma}$
with that of $\left| l,l\right\rangle$ using the contraction formula,
thereby transforming the element to

\begin{equation}
\label{SphericaltimesrElem}
\begin{gathered}
\left\langle l',m'\right|Y_{l_{\gamma}}^{m_\gamma}\left(\hat{r}\right)_{q}\left|l,l\right\rangle =\left(-1\right)^{l+m_{\gamma}}\sqrt{\frac{\left(2l_{\gamma}+1\right)\left(2l+1\right)}{4\pi}}\\
\times\sum_{c=m_{\gamma}}^{l_{\gamma}}\sqrt{2\left(l+c\right)+1}\left(\begin{array}{ccc}
l+c & l & l_{\gamma}\\
0 & 0 & 0
\end{array}\right)\left(\begin{array}{ccc}
l+c & l & l_{\gamma}\\
-\left(l+m_{\gamma}\right) & l & m_{\gamma}
\end{array}\right)\left\langle l',m'\right|\left(\hat{r}\right)_{q}\left|l+c,l+m_{\gamma}\right\rangle .
\end{gathered}
\end{equation}

We then use the relation <cit.>

\begin{equation}
\begin{gathered}
\left\langle l',m'\right|\left(\hat{r}\right)_{q}\left|l_{c},m_{c}\right\rangle =\delta_{m',m_{c}+q}\left\{ \left(\sqrt{\frac{\left(l_{c}-m_{c}+1\right)\left(l_{c}+m_{c}+1\right)}{\left(2l_{c}+1\right)\left(2l_{c}+3\right)}}\delta_{q,0}+\sqrt{\frac{\left(l_{c}\pm m_{c}+1\right)\left(l_{c}\pm m_{c}+2\right)}{2\left(2l_{c}+1\right)\left(2l_{c}+3\right)}}\delta_{q,\pm1}\right)\delta_{l',l_{c}+1}\right.\\
+\left.\left(\sqrt{\frac{\left(l_{c}-m_{c}\right)\left(l_{c}+m_{c}\right)}{\left(2l_{c}+1\right)\left(2l_{c}-1\right)}}\delta_{q,0}-\sqrt{\frac{\left(l_{c}\mp m_{c}-1\right)\left(l_{c}\mp m_{c}\right)}{2\left(2l_{c}+1\right)\left(2l_{c}-1\right)}}\delta_{q,\pm1}\right)\delta_{l',l_{c}-1}\right\} .
\end{gathered}
\end{equation}

§.§.§ Spherical Harmonic times Gradient

Finally, the matrix element $\left\langle l',m'\right|Y_{l_\gamma}^{m_\gamma}\left(\hat{r}\right)\left(\vec{\nabla}_{\Omega}\right)_{q}\left|l,l\right\rangle $
can be calculated using the expression <cit.>

\begin{equation}
\label{SphericaltimesGradientElem}
\begin{gathered}
\left\langle l',m'\right|Y_{l_\gamma}^{m_\gamma}\left(\hat{r}\right)\left(\vec{\nabla}_{\Omega}\right)_{q}\left|l,l\right\rangle =\left(-1\right)^{m_\gamma+q}\sqrt{\frac{\left(2l'+1\right)\left(2l_\gamma+1\right)}{4\pi\left(2l+1\right)}} \\
\times\sum_{s=\pm1}\left(l+\frac{1+3s}{2}\right)\sqrt{\frac{l+\frac{1+s}{2}}{2l+1+2s}}C_{l'\,0\,l_\gamma\,0}^{l+s\,0}C_{l'\,m'\,l_\gamma\,-m_\gamma}^{l+s\,m'-m_\gamma}C_{l+s\,m'-m_\gamma\,1\,-q}^{l\,l},
\end{gathered}
\end{equation}

where the $C$ coefficients are the Clebsch-Gordan coefficients.

§.§ Radial Matrix Elements

For hydrogen-like atoms, we can express the radial wave-function using the Kummer confluent hypergeometric function $\,_{1}F_{1}$ as follows

\begin{equation} \label{HydrogenKummer}
R_{n,l}\left(r\right)=\frac{1}{\left(2l+1\right)!}\sqrt{\frac{\left(n+l\right)!}{\left(n-l-1\right)!2n}}\left(\frac{2Z}{na_{0}}\right)^{l+3/2}\exp\left(-\frac{Zr}{na_{0}}\right)r^{l}\, \,_{1}F_{1}\left(-n+l+1;2l+2;\frac{2Zr}{na_{0}}\right).
\end{equation}

where $a_{0}=\frac{4\pi \hbar^{2}}{\mu  q^{2}}$ is the Bohr radius, and $Z$ is the atomic number. By using Gordon's integral <cit.>,

\begin{equation} \label{IntIdenKummer}
\int_{0}^{\infty}e^{-sx}r^{\rho-1}\,_{1}F_{1}\left(a;b;pr\right)\,_{1}F_{1}\left(c;d;qr\right)\,dr =s^{-\rho}\,\Gamma\left(\rho\right)F_{2}\left(\rho,a,c,b,d;\frac{p}{s},\frac{q}{s}\right),
\end{equation}

one can explicitly calculate matrix elements of the form $\left\langle n',l'\right|r^{j}\left|n,l\right\rangle$ and $\left\langle n',l'\right|r^{j}\partial_r\left|n,l\right\rangle$.

§.§.§ ”Keplerian" Matrix Element

\begin{eqnarray}~\label{eq:quantkep}
&&\left\langle n',l\right| r \left| n ,l\right\rangle=\frac{\hbar^{2}}{\mu  K}\frac{\left(-1\right)^{n'-l}2^{2l+2}\left(nn'\right)^{l+2}\left(n-l-1\right)}{\left(2l+1\right)!\left(n+n'\right)^{2l+4}}\left(\frac{n-n'}{n+n'}\right)^{n-n'-2}\sqrt{\frac{\left(n+l\right)!\left(n'+l\right)!}{\left(n-l-1\right)!\left(n'-l-1\right)!}}\times\nonumber\\[5pt]
\end{eqnarray}

and specifically

\begin{eqnarray}\label{eq:nteqn}
&&\left\langle n,l\right| r \left| n ,l\right\rangle=\frac{\hbar^{2}}{\mu  K}\frac{3n^2-l(l+1)}{2}\,.
\end{eqnarray}

§.§.§ General Transition

\begin{equation} \label{RadialMat1}
\begin{gathered}
\left\langle n',l'\right|r^{j}\left|n,l\right\rangle =\left(\frac{\hbar^{2}}{\mu  K}\right)^{j}\sqrt{\frac{\left(n+l\right)!\left(n'+l'\right)!}{\left(n-l-1\right)!\left(n'-l'-1\right)!}}
\times F_{2}\left(l+l'+j+3,-n+l+1,-n'+l'+1,2l+2,2l'+2;\frac{2n'}{n+n'},\frac{2n}{n+n'}\right),
\end{gathered}
\end{equation}

where $F_{2}$ is the Appell hypergeometric function of two variables, defined as:

\begin{equation}
\end{equation}

Next, we wish to calculate $\left\langle n',l'\right|r^{j}\partial_{r}\left|n,l\right\rangle$. Using Eq. <ref>, it is straightforward to show that

\begin{equation}
\begin{gathered}
\partial_{r}R_{n,l}\left(r\right)=-\frac{Z}{na_{0}}R_{n,l}\left(r\right)+\frac{l}{r}R_{n,l}\left(r\right)
\times\exp\left(-\frac{Zr}{a_{0}n}\right)r^{l}\,\,_{1}F_{1}\left(-n+l+2;2l+3;\frac{2Zr}{a_{0}n}\right).
\end{gathered}
\end{equation}

Hence, we have that

\begin{equation}\label{eq:partialrME}
\begin{gathered}
\hbar\left\langle n',l'\right|r^{j}\partial_{r}\left|n,l\right\rangle =-\frac{\mu K}{ N}\left\langle n',l'\right|r^{j}\left|n,l\right\rangle +L\left\langle n',l'\right|r^{j-1}\left|n,l\right\rangle
-\frac{2\mu K}{N}\frac{n-l-1}{2l+2}\left(\frac{a_{0}}{Z}\right)^{j}\\
\times\sqrt{\frac{\left(n+l\right)!\left(n'+l'\right)!}{\left(n-l-1\right)!\left(n'-l'-1\right)!}} 2^{l+l'+2}\frac{\left(l+l'+j+2\right)!}{\left(2l+1\right)!\left(2l'+1\right)!}\frac{n^{l'+j+1}\left(n'\right)^{l+j+1}}{\left(n+n'\right)^{l+l'+j+3}}\\
\times F_{2}\left(l+l'+j+3,-n+l+2,-n'+l'+1,2l+3,2l'+2;\frac{2n'}{n+n'},\frac{2n}{n+n'}\right).
\end{gathered}
\end{equation}

§.§.§ Bessel Matrix Elements

Using the recursion relations for spherical Bessel functions:

\begin{equation}\label{eq:rad2Q}
\left\langle n',l'\right|r^{-1}j_{l_\gamma}\left(\omega_{\Delta n}\,r\right)\left|n,l\right\rangle=\frac{\omega_{\Delta n}}{2l_\gamma+1}~\left\{\,\left\langle n',l'\right|j_{l_\gamma+1}\left(\omega_{\Delta n} \,r\right)\left|n,l\right\rangle+\left\langle n',l'\right|j_{l_\gamma-1}\left(\omega_{\Delta n} \,r\right)\left|n,l\right\rangle\,\right\}\,.
\end{equation}

On the other hand, Taylor expanding the spherical Bessel gives

\begin{eqnarray} \label{eq:rad1Q}
\left\langle n',l'\right|j_{l_\gamma}\left(\omega_{\Delta n} \,r\right)\left|n,l\right\rangle&=&2^{l_\gamma}\sum_{j=0}^{\infty}\frac{\left(-1\right)^{j}\left(j+l_\gamma\right)!}{j!\left(2j+2l_\gamma+1\right)!}\,\omega_{\Delta n}^{2j+l_\gamma} \left\langle n',l'\right|r^{2j}\left|n,l\right\rangle\nonumber\\[5pt]
\left\langle n',l'\right|j_{l_\gamma}\left(\omega_{\Delta n} \,r\right)\partial_r\left|n,l\right\rangle&=&2^{l_\gamma}\sum_{j=0}^{\infty}\frac{\left(-1\right)^{j}\left(j+l_\gamma\right)!}{j!\left(2j+2l_\gamma+1\right)!}\,\omega_{\Delta n}^{2j+l_\gamma} \left\langle n',l'\right|r^{2j}\partial_r\left|n,l\right\rangle\,.
\end{eqnarray}

§ CLASSICAL LIMITS OF QUANTUM MATRIX ELEMENTS

§.§ Angular Matrix Elements

§.§.§ Spherical Harmonic

By explicit calculation, the classical limit of (<ref>) is given by

\begin{equation}\label{eq:ang1}
\begin{gathered}
\lim_{\hbar\rightarrow 0}\left\langle l',m'\right|Y_{l_\gamma,m}^{m_\gamma}\left(\hat{r}\right)\left|l,l\right\rangle =\delta_{l',m'}\delta_{-\Delta l,m_\gamma}\,f_{l_\gamma,m_\gamma}\,,
\end{gathered}
\end{equation}

\begin{equation}
f_{l_\gamma,m_\gamma}\equiv\frac{\cos\left[\tfrac{\pi(l_\gamma-m_\gamma)}{2}\right]}{2\pi}\sqrt{\frac{\left(2 l_\gamma+1\right) \Gamma\left(\frac{l_\gamma+m_\gamma+1}{2}\right) \Gamma\left(\frac{l_\gamma-m_\gamma+1}{2}\right)}{\Gamma\left(\frac{l_\gamma+m_\gamma}{2}+1\right) \Gamma\left(\frac{l_\gamma-m_\gamma}{2}+1\right)}}
\end{equation}

§.§.§ Spherical Harmonic times $\hat{r}$

Again by explicit calculation, the classical limit of (<ref>) is

\begin{equation}\label{eq:ang2}
\begin{gathered}
\lim_{\hbar\rightarrow 0}\left\langle l',m'\right|Y_{l_\gamma}^{m_\gamma}\left(\hat{r}\right)\left(\hat{r}\right)_{q}\left|l,l\right\rangle =\delta_{l',m'}\delta_{-\Delta l,m_\gamma+q}\,\frac{1}{\sqrt{2}}\left(\delta_{q,1}-\delta_{q,-1}\right)\,f_{l_\gamma,m_\gamma}\,.
\end{gathered}
\end{equation}

§.§.§ Spherical Harmonic times Gradient

Finally, the classical limit of (<ref>) is

\begin{equation}\label{eq:ang3}
\begin{gathered}
\lim_{\hbar\rightarrow 0}\,\hbar\left\langle l',m'\right|Y_{l_\gamma}^{m_\gamma}\left(\hat{r}\right)\left(\vec{\nabla}_{\Omega}\right)_{q}\left|l,l\right\rangle =-\frac{L}{\sqrt{2}}\delta_{l',m'}\delta_{-\Delta l,m_\gamma+q}\,\left(\delta_{q,1}+\delta_{q,-1}\right)\,f_{l_\gamma,m_\gamma}\,.
\end{gathered}
\end{equation}

§.§ Radial Matrix Elements

§.§.§ ”Keplerian" Matrix Element

Starting from the matrix element (<ref>), we can use (<ref>) to simplify the part in square brackets. We get

\begin{eqnarray}~\label{eq:quantkepcl}
(-1)^{l+n'}\frac{\Gamma(n-l-1) \,\Gamma(2 l+2)}{\Gamma(n'+l+1)\, \Gamma(\Delta n)}\left[\,f_{\Delta n-1}\,-\,\frac{(n+l+1)(n-l)}{(n+n')^{2}}\frac{\Delta n}{ (\Delta n+1)}\,f_{\Delta n+1}\right]\,.
\end{eqnarray}

In the $\hbar\rightarrow 0$ limit, and using the explicit form of $f_{\Delta k}$, we have

\begin{eqnarray}~\label{eq:quantkepcl2}
\left[\ldots\right]&=&-
(-1)^{l+n'}\frac{\Gamma(n-l-1) \,\Gamma(2 l+2)}{\Gamma(n'+l+1)}{\left(\frac{2}{e\Delta n}\right)}^{\Delta n-1}\left[J_{\Delta n-1}(e\Delta n)\,-\,J_{\Delta n+1}(e\Delta n)\right]\nonumber\\[5pt]
(-1)^{l+n'}\frac{\Gamma(n-l-1) \,\Gamma(2 l+2)}{\Gamma(n'+l+1)}{\left(\frac{2}{e\Delta n}\right)}^{\Delta n}e\frac{d}{de}J_{\Delta n}(e\Delta n)\,.
\end{eqnarray}

Note that the Bessel function is evaluated at $\sqrt{1-\frac{L^2}{N^2}}\,\Delta n=e\Delta n$. Substituting this expression in (<ref>) and taking $\hbar\rightarrow 0$, we have

\begin{eqnarray}~\label{eq:resquantclskep}
\lim_{\hbar\rightarrow 0}\left\langle n',l\right| r \left| n ,l\right\rangle&=&-\frac{N^2}{\mu K}\frac{e}{\Delta n^2}\frac{d}{de}J_{\Delta n}(e\Delta n)\nonumber\\[5pt]
&=&-\frac{p}{1-e^2}\frac{e}{\Delta n^2}\frac{d}{de}J_{\Delta n}(e\Delta n)\,.
\end{eqnarray}

For $n'= n$ we can take the classical limit of (<ref>) directly and get

\begin{eqnarray}~\label{eq:resquantclskepeq}
\lim_{\hbar\rightarrow 0}\left\langle n,l\right| r \left| n ,l\right\rangle=\frac{p}{1-e^2}\left(1+\frac{e^2}{2}\right)\,.
\end{eqnarray}

§.§.§ General Transition

To calculate $\lim_{\hbar\to0}\left\langle n',l'\right|r^{j}\left|n,l\right\rangle$, from (<ref>), we make use of the Taylor expanded Appell $F_2$, (<ref>), with the parameters

\begin{eqnarray}
&&a=\frac{2L}{\hbar}-\Delta l+j+3~~~,~~~b_1=\frac{N-L}{\hbar}-1
~~~,~~~b_2=\frac{N-L}{\hbar}+\Delta l-\Delta n-1\nonumber\\[5pt]
&&c_1=\frac{2L}{\hbar}+2~~~,~~~c_2=\frac{2L}{\hbar}-2\Delta l+2~~~,~~~\Delta=\hbar \frac{\Delta n}{2N}\,.
\end{eqnarray}

The $\hbar\rightarrow 0$ limit is then straightforward, using the limit of the Pochhammer symbol (<ref>).
In the physical cases that we encounter in this work, at least one of $j\pm \Delta l+1\ge 0$ is positive, so we specialize to such cases. Assuming $j\pm \Delta l+1\ge 0$, the limit is then

\begin{eqnarray}
\left\langle n',l'\right|r^{j}\left|n,l\right\rangle =&&\left(\frac{p}{1-e^2}\right)^{j}\left(-\eta\right)^{-\Delta n-\Delta l}\left(\frac{\eta e}{2}\right)^{j+1}\sum_{m=0}^{\infty}\,\eta^{-2m}\times\nonumber\\[5pt]
j-\Delta l+1+s\\
\end{array}\right)\left(\frac{\eta e\Delta n}{2}\right)^s\right]\,\left[\sum_{r=0}^\infty\,\tfrac{1}{r!}\left(\begin{array}{c}
j+\Delta l+1+r\\
\Delta l+\Delta n + m
\end{array}\right)\left(-\frac{\eta e\Delta n}{2}\right)^r\right]\,\nonumber\\
\end{eqnarray}

where $\Delta n=n-n'$ and $\Delta l=l-l'$. The $r$ and $s$ sums can be carried analytically, giving rise to our final expression

\begin{eqnarray} \label{ClassLimRadialLaguerre2}
\left\langle n',l'\right|r^{j}\left|n,l\right\rangle =&&\left(\frac{p}{1-e^2}\right)^{j}\left(-\eta\right)^{-\Delta n-\Delta l}\left(\frac{\eta e}{2}\right)^{j+1}\times\nonumber\\[5pt]
&&\sum_{m=0}^{\infty}L_{m+\Delta l+\Delta n }^{j+1-m-\Delta n}\left(\frac{\eta e\Delta n}{2}\right)L_{m}^{j+1-m-\Delta l}\left(-\frac{\eta e \Delta n}{2}\right)\eta^{-2m}\,.
\end{eqnarray}

Otherwise, if $(j+\Delta l+1)(j-\Delta l+1)< 0$, then the classical limit is obtained similarly using (<ref>), albeit a less compact expression is obtained.

Next, we calculate $\lim_{\hbar\to0}\left\langle n',l'\right|r^{j}\partial_r\left|n,l\right\rangle$, from (<ref>), using the same asymptotic form of the Appell $F_2$ function from (<ref>). Assuming $j \ge \rm{Max}$ $\left\{\Delta l,-\Delta l-1\right\} $, we get

\begin{equation} \label{RadialDerivativeClassLim}
\begin{gathered}
\lim_{\hbar\to0}\hbar\left\langle n',l'\right|r^{j}\partial_{r}\left|n,l\right\rangle =L\left\{\frac{2E}{K\sqrt{1-e^2}}\left\langle n',l'\right|r^{j}\left|n,l\right\rangle +\left\langle n',l'\right|r^{j-1}\left|n,l\right\rangle
+\frac{2}{\sqrt{1-e^2}}\left(\frac{p}{1-e^2}\right)^{j-1}\left(-\eta\right)^{-\Delta n-\Delta l}\left(\frac{\eta e}{2}\right)^{j+1}\right.\\
\left.\times\sum_{m=0}^{\infty}L_{m+\Delta l+\Delta n }^{j+1-m-\Delta n}\left(\frac{\eta e\Delta n}{2}\right)L_{m}^{j-m-\Delta l}\left(-\frac{\eta e \Delta n}{2}\right)\eta^{-2m} \right\}\,.
% \times\sum_{k_{i}=0}^{\infty}\frac{\left(-\frac{\Delta n}{2}\right)^{k_{i}}\left(\eta e\right)^{k_{i}}}{k_{i}!}\sum_{l_i=0}^{k_{i}}\left(-1\right)^{l_i}\left(\begin{array}{c}
% k_{i}\\
% l_i
% \end{array}\right)\left(\begin{array}{c}
% \Delta l+j+k_{i}-l_i\\
% \Delta n+\Delta l
% \end{array}\right)\\
% \left.\times\,\,_{2}F_{1}\left(\Delta l-j-l_i-1,\Delta n+l_i-k_{i}-j,\Delta n+\Delta l+1;\eta^{-2}\right).
\end{gathered}
\end{equation}

Otherwise, the resulting expression is less compact.

§.§.§ Bessel Matrix Elements

Since the prefactors in (<ref>)-(<ref>) are finite in the classical limit, we can commute the $\hbar\rightarrow 0$ limit past them and get

\begin{equation}\label{eq:rad2}
\lim_{\hbar\rightarrow 0}\left\langle n',l'\right|r^{-1}j_{l_\gamma}\left(\omega_{\Delta n}\,r\right)\left|n,l\right\rangle=\frac{\omega_{\Delta n}}{2l_\gamma+1}~\left\{\,\lim_{\hbar\rightarrow 0}\left\langle n',l'\right|j_{l_\gamma+1}\left(\omega_{\Delta n} \,r\right)\left|n,l\right\rangle+\lim_{\hbar\rightarrow 0}\left\langle n',l'\right|j_{l_\gamma-1}\left(\omega_{\Delta n} \,r\right)\left|n,l\right\rangle\,\right\}\,,
\end{equation}

as well as

\begin{eqnarray}\label{eq:rad1}
\lim_{\hbar\rightarrow 0}\left\langle n',l'\right|j_{l_\gamma}\left(\omega_{\Delta n} \,r\right)\left|n,l\right\rangle&=&2^{l_\gamma}\sum_{j=0}^{\infty}\frac{\left(-1\right)^{j}\left(j+l_\gamma\right)!}{j!\left(2j+2l_\gamma+1\right)!}\,\omega_{\Delta n}^{2j+l_\gamma} \lim_{\hbar\rightarrow 0}\left\langle n',l'\right|r^{2j+l_\gamma}\left|n,l\right\rangle
\end{eqnarray}

\begin{eqnarray}\label{eq:rad3}
\lim_{\hbar\rightarrow 0}\left\langle n',l'\right|j_{l_\gamma}\left(\omega_{\Delta n} \,r\right)\partial_r\left|n,l\right\rangle&=&2^{l_\gamma}\sum_{j=0}^{\infty}\frac{\left(-1\right)^{j}\left(j+l_\gamma\right)!}{j!\left(2j+2l_\gamma+1\right)!}\,\omega_{\Delta n}^{2j+l_\gamma} \lim_{\hbar\rightarrow 0}\left\langle n',l'\right|r^{2j+l_\gamma}\partial_r\left|n,l\right\rangle\,.
\end{eqnarray}

§ 0PA ENERGY LOSS FROM SPONTANEOUS EMISSION

In section <ref>, we got the following expression for the 0PA EM energy loss:

\begin{eqnarray}\label{eq:dE4app}
\left\langle\frac{dE}{dt}\right\rangle_{\alpha}&=&\frac{iq^2}{\mu^2}\,\lim_{\hbar\rightarrow 0}\sum_{\Delta n,\Delta l,\Delta m}\omega_{\Delta n}\,\sum_{l_{\gamma}=0}^\infty\sum_{m_{\gamma}=-l_{\gamma}}^{l_\gamma}\,\left\{ \,i\omega_{\Delta n}\,\left|\left\langle n' l' m'\right| j_{l_{\gamma}}(\omega_{\Delta n} r)\,Y_{l_\gamma}^{m_\gamma}\,p^i\left|nlm\right\rangle\right|^2-\right.\nonumber\\[5pt]
&&~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\left.\mu\,\left\langle  n' l' m'\right|\partial_i\left[j_{l_{\gamma}}(\omega_{\Delta n} r)\,Y_{l_\gamma}^{m_\gamma}\right]\,p^i\left|nlm\right\rangle^*\,\left\langle n'l'm'\right|j_{l_{\gamma}}(\omega_{\Delta n} r)\,Y_{l_\gamma}^{m_\gamma}\left|nlm\right\rangle\right\}\,,\nonumber\\
\end{eqnarray}

where $\omega_{\Delta n}=\Upsilon\Delta n$. Here we show that this expression is the classical limit of the quantum expression for spontaneous emission. To see this, we first look at the expression

\begin{eqnarray}\label{eq:dE4app1}
\sum_{l_{\gamma}=0}^\infty\sum_{m_{\gamma}=-l_{\gamma}}^{l_\gamma}\,\omega_{\Delta n}\,\left|\left\langle n' l' m'\right| j_{l_{\gamma}}(\omega_{\Delta n} r)\,Y_{l_\gamma}^{m_\gamma}\,p^i\left|nlm\right\rangle\right|^2=\frac{1}{(4\pi)^2}\int\,\frac{d^3k}{|k|}\,\delta\left(|k|-\omega_{\Delta n}\right)\,\left|\left\langle n' l' m'\right| e^{-i\vec{k}\cdot\vec{p}}\,p^i\left|nlm\right\rangle\right|^2\,.\nonumber\\
\end{eqnarray}

The equality here is a consequence of the partial wave decomposition of the $e^{-i\vec{k}\cdot\vec{p}}$ plane wave, followed by integration over $\vec{k}$. Similarly, we have

\begin{eqnarray}
&&\mu\,\left\langle  n' l' m'\right|\partial_i\left[j_{l_{\gamma}}(\omega_{\Delta n} r)\,Y_{l_\gamma}^{m_\gamma}\right]\,p^i\left|nlm\right\rangle^*\,\left\langle n'l'm'\right|j_{l_{\gamma}}(\omega_{\Delta n} r)\,Y_{l_\gamma}^{m_\gamma}\left|nlm\right\rangle\nonumber\\[5pt]
&&=\frac{\mu}{(4\pi)^2}\int\,\frac{d^3k}{|k|}\,\frac{1}{\omega_{\Delta n}}\delta\left(|k|-\omega_{\Delta n}\right)\,{\left\langle n' l' m'\right| \partial_i e^{-i\vec{k}\cdot\vec{p}}\,p^i\left|nlm\right\rangle}^*\left\langle n' l' m'\right| e^{-i\vec{k}\cdot\vec{p}}\left|nlm\right\rangle\nonumber\\[5pt]
&&=-\frac{i}{(4\pi)^2}\frac{\mu}{\mu+\hbar\omega_{\Delta n}}\int\,\frac{d^3k}{|k|}\,\delta\left(|k|-\omega_{\Delta n}\right)\,\left|\left\langle n' l' m'\right| e^{-i\vec{k}\cdot\vec{p}}\hat{k}\cdot \vec{p
\end{eqnarray}

Putting it all together, we have

\begin{eqnarray}\label{eq:dE4app2}
\left\langle\frac{dE}{dt}\right\rangle_{\alpha}&=&-\lim_{\hbar\rightarrow 0}\sum_{\Delta n,\Delta l,\Delta m}\,\frac{q^2\omega_{\Delta n}}{(4\pi)^2\mu^2}\,\int\,d^3k\,\delta(|k|-\omega_{\Delta n})\left\{\left|\left\langle n' l' m'\right| e^{-i\vec{k}\cdot\vec{p}}\,p^i\left|nlm\right\rangle\right|^2-\left|\left\langle n' l' m'\right| e^{-i\vec{k}\cdot\vec{p}}\hat{k}\cdot \vec{p}\left|nlm\right\rangle\right|^2\right\}\nonumber\\[5pt]
&=&-\,\lim_{\hbar\rightarrow 0}\sum_{\Delta n>0,\Delta l,\Delta m}\,\frac{q^2\,\left(E_{n}-E_{n'}\right)}{8\pi^2\mu^2\hbar}\,\int\,d^3k\,\delta(|k|-\omega_{\Delta n})\,\sum_{\sigma=\pm}\,\left|\left\langle n' l' m'\right| e^{-i\vec{k}\cdot\vec{p}}\,\left(\vec{\varepsilon}^{~\sigma}_{\vec{k}}\cdot\vec{p}\right)\left|nlm\right\rangle\right|^2\,,
\end{eqnarray}

where $\vec{\varepsilon}^{~\sigma}_{\vec{k}}$ are the usual EM transverse polarization vectors. This is the classical limit of the famous spontaneous emission formula.

§ AUXILIARY CALCULATIONS

Here we derive for completeness various asymptotic formulae  relevant for the classical limit of the radial matrix elements. These include asymptotic forms of the Pochhammer symbol, the Gauss hypergeometric function ${}_2F_1$, and the Appell function $F_2$.

§.§ Pochhammer Symbol

As a warm-up, we calculate the classical limit of the Pochhammer Symbol,

\begin{eqnarray}\label{eq:Poclim}
\lim_{\hbar \rightarrow 0}\hbar^m(A\hbar^{-1})_m=A^m\,,
\end{eqnarray}

where $(a)_m$ is the Pochhammer symbol.

§.§ Gauss Hypergeometric

Now we move on to calculate the classical limit of our ubiquitous Gauss Hypergeometric functions. Consider the following hypergeometric function

\begin{eqnarray}
\end{eqnarray}

where $i,j,k$ are positive integers and $\Delta\ll 1$. We wish to express it in the form

\begin{eqnarray}\label{eq:fexp}
\end{eqnarray}

To do this, we Taylor expand in $\Delta$, and use the identity

\begin{eqnarray}
\frac{d^{l}\,{}_{2}F_{1}(x)}{dx^{l}}=\frac{(-i)_l(j)_l}{(k)_l}{}_2F_1\left[-i+l,j+l,k+l,x\right]\,.
\end{eqnarray}

After Taylor expanding, we have

\begin{eqnarray}\label{eq:prechu}
\end{eqnarray}

where in the last line, we used a $z\rightarrow 1/z$ Euler transformation (in the polynomial case). We now simplify this further using the Chu-Vandermonde identity,

\begin{eqnarray}\label{eq:chu}
\end{eqnarray}

so that $f(\Delta)=\sum_l\,\frac{(-i)_l(j)_l}{l!\,(-i+j-k+1)_l}(-\Delta)^l$.

Now we are ready to take the classical limit of the hypergeometric function. More specifically, we consider the case in which $(i,j)=\hbar^{-1}(I,J)$, $k=\hbar^{-1}(-I+J)-\Delta k$ and $\Delta =-\delta \hbar^2$, with $\hbar\rightarrow 0$. In this case we have

\begin{eqnarray}\label{eq:fd}
\lim_{\hbar\rightarrow 0}\,f(\Delta )&=&\lim_{\hbar\rightarrow 0}\,\sum_l\,\frac{(-I\hbar^{-1})_l(J\hbar^{-1})_l}{l!\,(\Delta k+1)_l}(\hbar^2 \delta)^l=\sum_l\,\frac{\left(-I J \delta\right)^l}{l!\,(\Delta k+1)_l}=(\Delta k)!{\left(\frac{2}{\kappa}\right)}^{\Delta k}J_{\Delta k}(\kappa)\equiv f_{\Delta k}\,,
\end{eqnarray}

where $\kappa=2\sqrt{IJ\delta}$. Here we used the classical limit of the Pochhammer symbol in the middle equality. Finally, we can express ${}_2F_1$ evaluated at $1$ using gamma functions and get:

\begin{eqnarray}\label{eq:fexpfin}
\lim_{\hbar\rightarrow 0}\,{}_2F_1\left[-i,j,k,1+\Delta\right]&=&\lim_{\hbar\rightarrow 0}\,\frac{(k-j)_i}{(k)_i}\,f_{\Delta k}=\lim_{\hbar\rightarrow 0}\,(-1)^i\frac{\Gamma(j-k+1)\Gamma(k+i)}{\Gamma(j-k-i+1)\Gamma(k)}\,f_{\Delta k}\,.
\end{eqnarray}

§.§ Appell $F_2$

We are looking to expand

\begin{equation}
\end{equation}

as a power series in $\Delta$. For this purpose we verify by explicit calculation that

\begin{eqnarray}
\frac{\partial^{r+s}}{\partial x^r\partial y^s}F_{2}\left(a,-b_1,-b_2,c_1,c_2;x,y\right)=
\frac{(a)_{r+s}(-b_1)_r(-b_2)_s}{(c_1)_r(c_2)_s}\,F_{2}\left(a+r+s,-b_1+r,-b_2+s,c_1+r,c_2+s;x,y\right)\,.
\end{eqnarray}

Using this expression for our Taylor expansion in $\Delta$, we have

\begin{eqnarray}\label{eq:F2Del}
\sum_{r,s=0}^\infty\,\frac{(a)_{r+s}(-b_1)_r(-b_2)_s}{(c_1)_r(c_2)_s\Gamma(r+1)\Gamma(s+1)}\times\nonumber\\[5pt]
\end{eqnarray}

To compute $F_2$ at $x=1,y=1$ we use the representation of $F_2$ as the sum of products of Gauss ${}_2F_1$:

\begin{eqnarray}
\end{eqnarray}

In the cases of interest in this work, $A-C_1$ and $A-C_2$ can not be both negative simultaneously. Without loss of generality, assume that $A-C_1\ge0$.  After a little massaging we get

\begin{equation}
\begin{gathered}
\times\sum_{m=0}^{A-C_{1}+B_{2}-B_{1}}\left(\begin{array}{c}
\end{array}\right)\left(\begin{array}{c}
\end{array}\right)\frac{\left(B_{1}+C_{1}+m-1\right)!}{\left(B_{1}+C_{1}-A+m\right)!}.
\end{gathered}
\end{equation}

Using this expression in (<ref>), we finally get

\begin{eqnarray}\label{eq:F2fin}
\sum_{r,s=0}^\infty\,\frac{(a)_{r+s}(-b_1)_r(-b_2)_s}{(c_1)_r(c_2)_s\Gamma(r+1)\Gamma(s+1)}\frac{\left(-1\right)^{b_{2}-b_{1}+r-s}(b_{1}-r)!(b_{2}-s)!}{\left(c_{1}+r\right)_{b_{1}-r}\left(c_{2}+s\right)_{b_{2}-s}\left(a+r+s-1\right)!}\times\nonumber\\[5pt]
\end{array}\right)\left(\begin{array}{c}
\end{array}\right)\frac{\left(b_{1}+c_{1}+m-1\right)!}{\left(b_{1}+c_{1}-a-r-s+m\right)!}(-\Delta)^r\Delta^s\,.
\end{eqnarray}

This expression might seem cumbersome, but it simplifies greatly upon taking the classical limit.