Transition of EMRIs through resonance: corrections to higher order in the on-resonance flux modification

TL;DR

This paper extends models of EMRI resonances by calculating higher-order corrections to the on-resonance flux, providing a more precise understanding of the gravitational self-force effects during resonance in extreme mass ratio inspirals.

Contribution

It introduces higher-order corrections to existing resonance models of EMRIs and validates the instantaneous frequency approach for analyzing resonance effects.

Findings

Extended the Gair, Bender, and Yunes model to third order in flux modification.

Developed a perturbative solution algorithm for non-linear resonance equations.

Discussed the mathematical properties and scope of the perturbative approach.

Abstract

Extreme mass ratio in-spirals (EMRIs) are candidate events for gravitational wave detection in the millihertz range (by detectors like LISA and eLISA). These events involve a stellar-mass black hole, or a similar compact object, descending in the gravitational field of a supermassive black hole, eventually merging with it. Properties of the in-spiralling trajectory away from resonance are well known and have been studied extensively, however little is known about the behaviour of these binary systems at resonance, when the radial and lateral frequencies of the orbit become commensurate. We describe the two existing models, the instantaneous frequency approach used by Gair, Bender, and Yunes, and the standard two timescales approach implemented by Flanagan and Hinderer. In both cases, the exact treatment depends on the modelling of the gravitational self-force, which is currently not available. We extend the results in Gair, Bender and Yunes to higher order in the on-resonance flux modification, and argue that the instantaneous frequency approach is also a valid treatment of the resonance problem. The non-linear differential equations which arise in treating resonances are interesting from a mathematical view point. We present our algorithm for perturbative solutions and the results to third order in the infinitesimal parameter, and discuss the scope of this approach.

Equations

$$a_{\alpha }=\frac{{\mathrm{d}}^2{x}_{\alpha }}{\mathrm{d}{\tau }^2}=ϵF_{\alpha }(\mathbf{q},\mathbf{J}),$$

$$\frac{{\omega }_r}{{\omega }_{\theta }}=\frac{n_{\theta }}{n_r},\{n_r,n_{\theta }\}\subset \mathbb{Z}.$$

$$\frac{\mathrm{d}{J}_{\mu }}{\mathrm{d}\tau }=ϵ{\mathcal{F}}_{\mu }(\mathbf{q},\mathbf{J},\mathbf{F})+\mathcal{O}\left({ϵ}^2\right),$$

$$\frac{\mathrm{d}{q}_{\mu }}{\mathrm{d}\tau }={\omega }_{\mu }+ϵ{\mathcal{H}}_{\mu }(\mathbf{q},\mathbf{\omega })+\mathcal{O}\left({ϵ}^2\right).$$

$$\frac{\mathrm{d}{\omega }_{\mu }}{\mathrm{d}\tau }=ϵ{\mathcal{G}}_{\mu }(\mathbf{q},\mathbf{\omega })+\mathcal{O}\left({ϵ}^2\right).$$

$$\displaystyle\begin{split}\frac{\mathop{}\!\mathrm{d}\omega_{\mu}}{\mathop{}\!\mathrm{d}\tau}=\epsilon\sum_{l,\,m}&\left(A_{lm}\cos\!\Big{[}\!\left(l\omega_{r}+m\omega_{\theta}\right)\tau\Big{]}\vphantom{\Big{[}\Big{]}^{\!\!\!2\,}}\right.\\ &\left.\vphantom{\Big{[}\Big{]}^{\!\!\!2\,}}\!\!\!+B_{lm}\sin\!\Big{[}\!\left(l\omega_{r}+m\omega_{\theta}\right)\tau\Big{]}\right)\!+\mathcal{O}\!\left(\epsilon^{2}\right).\end{split}$$

$$\displaystyle\begin{split}\frac{\mathop{}\!\mathrm{d}^{2}\phi_{\mu}}{\mathop{}\!\mathrm{d}\tau^{2}}=\epsilon\sum_{l,\,m}&\,\Big{(}C_{lm}\cos\left(l\phi_{r}+m\phi_{\theta}\right)\Big{.}\\ &\Big{.}\!\!\!+D_{lm}\sin\left(l\phi_{r}+m\phi_{\theta}\right)\!\Big{)}\!+\mathcal{O}\!\left(\epsilon^{2}\right)\end{split}$$

$$\displaystyle\begin{split}\frac{\mathop{}\!\mathrm{d}\bar{\omega}}{\mathop{}\!\mathrm{d}\tau}=\epsilon\sum_{l,\,m}&\left(A_{lm}\cos\!\Big{[}\!\left(l\bar{\omega}+m\omega_{r}\right)\tau\Big{]}\vphantom{\Big{[}\Big{]}^{\!\!\!2\,}}\right.\\ &\!\!\!\left.\vphantom{\Big{[}\Big{]}^{\!\!\!2\,}}+B_{lm}\sin\!\Big{[}\!\left(l\bar{\omega}+m\omega_{r}\right)\tau\Big{]}\right)\!+\mathcal{O}\!\left(\epsilon^{2}\right).\end{split}$$

$$\displaystyle\begin{split}\frac{\mathop{}\!\mathrm{d}^{2}\bar{\phi}}{\mathop{}\!\mathrm{d}\tau^{2}}=\epsilon\sum_{l,\,m}&\,\Big{(}C_{lm}\cos\left(l\bar{\phi}+m\phi_{r}\right)\Big{.}\\ &\!\!\!\Big{.}+D_{lm}\sin\left(l\bar{\phi}+m\phi_{r}\right)\!\Big{)}\!+\mathcal{O}\!\left(\epsilon^{2}\right).\end{split}$$

$$\displaystyle\hat{\omega}=\frac{\bar{\omega}}{\sqrt{\epsilon}},\;\;\;\hat{x}=\sqrt{\epsilon}\,\tau\;\;\;\text{in eq.~{}(\ref{eq:freq_general2})}$$

$$\displaystyle\hat{\phi}=\bar{\phi},\;\;\;\hat{\phi}_{r}=\sqrt{\epsilon}\,\phi_{r},\;\;\;\hat{x}=\sqrt{\epsilon}\,\tau\;\;\;\text{in eq.~{}(\ref{eq:phase_general2})}.$$

$$\displaystyle\begin{split}\frac{\mathop{}\!\mathrm{d}\hat{\omega}}{\mathop{}\!\mathrm{d}\hat{x}}=\sum_{l,\,m}&\left(A_{lm}\cos\!\Bigg{(}\!l\hat{x}\hat{\omega}+\frac{m\omega_{r}\hat{x}}{\sqrt{\epsilon}}\!\Bigg{)}\vphantom{\Bigg{[}\Bigg{]}^{\!\!\!2\,}}\right.\\ &\left.\vphantom{\Bigg{[}\Bigg{]}^{\!\!\!2\,}}\!\!+B_{lm}\sin\!\Bigg{(}\!l\hat{x}\hat{\omega}+\frac{m\omega_{r}\hat{x}}{\sqrt{\epsilon}}\!\Bigg{)}\!\!\right)\!+\mathcal{O}\!\left(\epsilon\right),\end{split}$$

$$\displaystyle\begin{split}\frac{\mathop{}\!\mathrm{d}^{2}\hat{\phi}}{\mathop{}\!\mathrm{d}\hat{x}^{2}}=\sum_{l,\,m}&\left(C_{lm}\cos\!\left(\!l\hat{\phi}+\frac{m\hat{\phi}_{r}}{\sqrt{\epsilon}}\!\right)\vphantom{\Bigg{[}\Bigg{]}^{\!\!\!2\,}}\right.\\ &\left.\vphantom{\Bigg{[}\Bigg{]}^{\!\!\!2\,}}\!\!\!+D_{lm}\sin\!\Bigg{(}\!l\hat{\phi}+\frac{m\hat{\phi}_{r}}{\sqrt{\epsilon}}\!\Bigg{)}\!\!\right)\!+\mathcal{O}\!\left(\epsilon\right).\end{split}$$

$$\frac{\mathrm{d}\hat{\omega }}{\mathrm{d}\hat{x}}$$

$$\frac{{\mathrm{d}}^2\hat{\varphi }}{\mathrm{d}{\hat{x}}^2}$$

$$\frac{\mathrm{d}\hat{\omega }}{\mathrm{d}\hat{x}}$$

$$\frac{{\mathrm{d}}^2\hat{\varphi }}{\mathrm{d}{\hat{x}}^2}$$

$$\begin{array}{rl}& {\mathcal{A}}_n=\sqrt{A_n^2+B_n^2},{\mathcal{B}}_n=\sqrt{C_n^2+D_n^2},\\ & \delta {\chi }_n=\arctan \left(\frac{B_n}{A_n}\right)\text{and}\delta {\psi }_n=\arctan \left(\frac{D_n}{C_n}\right).\end{array}$$

$$\frac{\mathrm{d}\hat{\omega }}{\mathrm{d}\hat{x}}$$

$$\frac{{\mathrm{d}}^2\hat{\varphi }}{\mathrm{d}{\hat{x}}^2}$$

$$\displaystyle\begin{split}\omega&=\big{(}\mathcal{A}_{0}\cos(\delta\chi_{0})\big{)}^{-\nicefrac{{1}}{{2}}}\>\hat{\omega}\\ x_{\omega}&=\big{(}\mathcal{A}_{0}\cos(\delta\chi_{0})\big{)}^{\nicefrac{{1}}{{2}}}\>\hat{x}\end{split}\quad\quad$$

$$\displaystyle\begin{split}x_{\phi}&=\big{(}\mathcal{B}_{0}\cos(\delta\psi_{0})\big{)}^{\nicefrac{{1}}{{2}}}\>\hat{x}\\ \phi&=\hat{\phi}\end{split}\quad\quad$$

$${\omega }^{\prime }$$

$${\varphi }^{\prime \prime }$$

$$\frac{\mathrm{d}{J}_{\mu }}{\mathrm{d}\tau }$$

$$\frac{\mathrm{d}{\omega }_{\mu }}{\mathrm{d}\tau }$$

$$\frac{\mathrm{d}{q}_{\mu }}{\mathrm{d}\tau }={\omega }_{\mu }\left(\mathbf{J}\right)+ϵ{\mathcal{K}}_{\mu }(\mathbf{q},\mathbf{J})+\mathcal{O}\left({ϵ}^2\right).$$

$$\bar{\nu }\left(x_{\nu },ϵ\right)$$

$$\bar{q}\left(x_{\nu },ϵ\right)$$

$$\frac{\mathrm{d}{\bar{\nu }}_1}{\mathrm{d}{x}_{\nu }}$$

$$\frac{\mathrm{d}{\bar{q}}_0}{\mathrm{d}{x}_{\nu }}$$

$${\omega }^{\prime }$$

$${\omega }^{\prime }(x)=1+\kappa \cos \left(x\omega \right),$$

$$\omega (x)=\sum _{n=0}^{\infty }{\kappa }^n{\varpi }_n(x).$$

Full text

I Introduction

Gravitational waves are one of the most fascinating predictions of general relativity. The theoretical underpinning for them was first presented by Einstein (1918). In 1974, Taylor and Hulse discovered indirect evidence for their existence when they observed the first pulsar in a binary system, in which the object separation was gradually decreasing, and were the first to suggest that gravitational wave effects are at play (Hulse and Taylor, 1975; Taylor and Weisberg, 1982). The pair were awarded the Nobel Prize in Physics in 1993 for their discovery – until recently the only evidence for the existence of gravitational waves. Nearly a century after Einstein’s paper, gravitational waves were directly observed for the first time in September 2015 by the two advanced ligo detectors, thus culminating a decades-long efforts to detect gravitational radiation (Abbott et al., 2016a). Currently, several experiments, including Advanced ligo and Advanced Virgo, which are searching for waves in the range ${10}^{2}\text{\,}\mathrm{Hz}{10}^{4}\text{\,}\mathrm{Hz}$ (The LIGO Scientific Collaboration et al., 2015; Acernese et al., 2015), are being deployed to search for further signals. The first joint run is scheduled to take place in late 2016 or early 2017 (Abbott et al., 2016b) and it is expected that a significant number of events will be observed in the near future (Abadie et al., 2010; The LIGO Scientific Collaboration et al., 2016). The European Space Agency has selected gravitational wave detection from space as the science theme to be addressed by the L3 large satellite mission, currently scheduled for launch in 2034. This space-based detector, called Laser Interferometer Space Antenna (lisa), will operate in the millihertz frequency range. The mission design is still being formulated, but it will be finalised by the end of this decade (Amaro-Seoane et al., 2012). The need for multiple approaches to detection (currently there are more than a dozen projects, both ground- and space-based, see Moore et al. (2015) for an overview) is driven by the broad spectrum of gravitational wave frequencies (the typical frequency of a gravitational wave system is given by $f\sim(GM/L^{3})^{1/2}$, where $M$ is the mass of the system, and $L$ is its length scale), namely in the range ${10}^{-16}\text{\,}\mathrm{Hz}{10}^{4}\text{\,}\mathrm{Hz}$, produced by a variety of potential gravitational radiation sources (Thorne, 1995).

Systems which are expected to produce gravitational waves are those that are accelerating and exhibit asymmetry in their configuration: for instance a binary (Aasi et al., 2014), a supernova (Gossan et al., 2016), or a spinning star with a pronounced bulge. Among the range of binary celestial systems, we focus on the case of a compact object (of several solar masses, usually a stellar-mass black hole, or a neutron star) gradually spiralling towards a supermassive black hole (typically between $10^{5}$ and $10^{12}$ solar masses, although the ones of interest to us are between $10^{5}$ and $10^{7}$ $M_{\odot}$, since these generate gravitational waves at frequencies in the range of sensitivity of lisa (Amaro-Seoane et al., 2012; Babak et al., 2015)). Since the mass $m$ of the compact object is in the range of $10-10^{2}M_{\odot}$, the ratio of the masses is typically $\epsilon=m/M\sim 10^{-6}\text{ to }10^{-4}$. Termed extreme-mass-ratio inspirals (emris), these binaries are expected to radiate gravitational waves with frequency in the range ${10}^{-4}\text{\,}\mathrm{Hz}{10}^{-2}\text{\,}\mathrm{Hz}$. This fact places them in a position to be likely candidates for detection by the lisa mission. emris reside near the centres of galaxies (where supermassive black holes are situated (Kormendy and Richstone, 1995)) and clusters of stars and stellar remnants are available for capture by the black hole. Plenty of scientific effort has been invested into modelling the evolution of emris (Barack and Cutler, 2004; Gair and Glampedakis, 2006; Drasco and Hughes, 2006; Chua and Gair, 2015), computing their waveforms (Babak et al., 2007), predicting the rates and timescales of inspiral (Glampedakis et al., 2002; Hopman, 2009), and estimating the parameters of the orbit (Barack and Cutler, 2004; Arun et al., 2009). Applications of gravitational wave observations of emris to inference on astrophysical populations (Gair et al., 2010), cosmology (MacLeod and Hogan, 2008) and fundamental physics (Gair et al., 2013) have also been explored in detail.

The mass of the compact object in an emri is typically at least 4 orders of magnitude smaller than the mass of the supermassive black hole. This allows us to treat the presence of the compact body in the background space-time of the central object as a perturbation. According to the uniqueness theorem for black holes, they can be uniquely characterised by three parameters - mass, spin and charge. Since electromagnetic repulsion is much stronger than gravitational attraction (about 40 orders of magnitude), we can assume that astrophysical black holes have negligible net charge. Hence, we consider the Kerr metric, where black holes are described only by their mass and spin. The Kerr space-time is well-known, with a complete set of integrals of motion, and the geodesics can be characterised by three constants (Carter (1968) and Appendix B). For a point mass, the trajectory would be that of a geodesic around the black hole, but while the point-mass notion is convenient mathematically, it does not apply in the context of the current problem. Therefore, we need to address the finite size and mass of the compact object. The deviation from geodesic motion can be quantified by allowing for the gravitational self-force on the compact object, (Nakamura and Sasaki, 1982; Wardell et al., 2012), also known as the gravitational radiation reaction force, since it is responsible for the gravitational radiation emission from the object. There is a radiative component of the self-force that leads to evolution of the orbital frequencies, as well as a conservative component that changes the instantaneous values of the orbital frequencies. The trajectory can be modelled as a sequence of geodesics, each one with lower energy than the previous one, and therefore (on average) closer to the central object. This method is called the osculating element approach (Pound and Poisson, 2008; Gair et al., 2011; Sopuerta and Yunes, 2012). Individual orbits are well approximated as geodesics, however the defining parameters of these geodesics are continuously changing, and this is what shapes the inspiral orbit. This effect will be visible over a large number of orbits – typically over a timescale of $\mathcal{O}\!\left(\epsilon^{-1}\right)$.

As the test body progresses along its in-spiral, at certain points two of the three frequencies of motion will become rational multiples of each other, putting the system into resonance (Flanagan and Hinderer, 2012). While in general the behaviour of the inspiral is well-studied and modelled (Barack, 2009), during resonance the system exhibits unusual behaviour, accompanied by changes in the orbital parameters (Flanagan and Hinderer, 2012). The evolution of emris during resonances, and their effect on gravitational wave emission is currently an interesting topic in gravitational wave astronomy (Berry et al., 2016).

There are two approaches in the literature that treat the behaviour of extreme-mass-ratio inspirals on resonance. On one hand is the instantaneous value treatment by (Gair et al., 2010), with manifestations in both frequency and phase formalisms, and on the other is the approach by (Flanagan and Hinderer, 2012), which separates the timescales associated with the in-spiral into small and large. We will describe both methods here, but focus on extending the former approach to higher order, as we anticipate the predictions of the latter to be consistent with the phase formalism in the former.

I.1 The Gravitational Self-force

The geodesic principle of General Relativity dictates that free point particles, of zero mass, traverse time-like geodesics of the space-time. Unfortunately, while treating bodies in General Relativity as point masses is mathematically convenient, in reality this approach is far from accurate, and introduces immense difficulties to the problem we are considering. As the size $R$ of a body is scaled down to 0, its mass $M$ should decrease accordingly, scaling to 0 at the same time as the size of the body (Gralla and Wald, 2008), otherwise a black hole solution would form before the point-limit of the body is reached. We can account for the finite size and mass of the body if we consider the effect of these on the deviation from geodesic motion. At leading order, we ignore the presence of the space-time generated by the compact object in the background space-time, and in this idealised system, the test body follows a time-like geodesic, whose shape can be accurately computed (see Appendix C and Figure C for details). The next order accounts for the non-infinitesimal size and mass of the object, and this causes a deviation from the geodesic motion described above. We quantify this phenomenon by introducing the gravitational self-force, responsible for the inward spiral of the compact object (Nakamura and Sasaki, 1982; Flanagan and Hinderer, 2012). This force is at the same time responsible for dissipating the energy made available by the inspiral motion in the form of gravitational radiation.

A proper treatment of the gravitational self-force is not the objective of this paper. The reader may wish to consult Poisson et al. (2011), as well as the more concise review Poisson (2011). Further recommended reading is Nakamura and Sasaki (1982) and Wald (2011), or the more mathematical approach described in Gralla and Wald (2008).

In the context of this article, it is worth knowing that since the self-force is the perturbation in the mass parameter to the geodesic, its leading order dependence on the mass of the test body is linear. Further, the self-force will depend on the shape of the orbit, which for the Kerr space-time is specified by three constants of the motion, for example the energy $\mathcal{E}$, the 3rd component of the angular momentum $\mathcal{L}_{z}$, and Carter’s constant, $\mathcal{Q}$. Furthermore, the individual components of this self-force will depend on the position of the test body along the orbit, specified, for example, by the action-angle coordinates $q_{\alpha}$, where the index $\alpha\in\{t,r,\theta,\phi\}$ (this implies there are 4 separate action-angle coordinates, with 4 corresponding components of the self-force. In fact, only 3 of these force components are independent, as can be verified by differentiating the condition $g_{\alpha\beta}u^{\alpha}u^{\beta}=1$ with respect to the proper time $\tau$ to obtain $u^{\alpha}a_{\alpha}=0$ (Gair et al., 2011)). With this in mind, we can assume that the linear order gravitational self-force can be written in the form (Poisson, 2011; Gair, Yunes, and Bender, 2012)

[TABLE]

where $\epsilon=m/M$ is the ratio of the compact object mass to the mass of the supermassive black hole, and the vector $\mathbf{J}$ contains the three constants of motion defining the shape of the orbit. The values of $\mathbf{q}$ and $\mathbf{J}$ in general depend on the way the orbit evolved before the moment at which they are computed. Thus, it can be regarded that the instantaneous self-force is in fact given by an integral over the past history of the orbit.

I.2 Orbital resonances

We already reasoned why the gravitational self-force depends on the shape of the orbit, i.e. the values of the three constants of motion in the Kerr metric. These quantities can be mapped to the three fundamental frequencies of orbital motion $\{\omega_{r},\omega_{\theta},\omega_{\phi}\}$ (Schmidt, 2002). Therefore, instead of eq. (1), we can write the expression for the self-force as $\epsilon\,\mathbf{F}(\mathbf{q},\bm{\omega})$. The fundamental frequencies are another valid way to label orbits of the Kerr space-time, and can be used in lieu of the geodesic constants of motion. It should be noted that the reverse substitution, from $\bm{\omega}$ to $\mathbf{J}$ is not always possible, on the account of the isofrequency pairs exhibited by the proper time frequencies (Warburton et al., 2013). For the discussion in this paper, we can assume local invertibility, however it should be noted that this is not universal.

The three fundamental frequencies evolve independently according to their equations of motion. However, at certain points along the trajectory, two of the frequencies will become commensurate, and the system will temporarily exhibit a qualitatively different behaviour, a resonance, with a reduced number of degrees of freedom. Since the Kerr metric (Appendix A and Figure 23) does not depend explicitly on the azimuthal angle $\phi$, the corresponding fundamental frequency $\omega_{\phi}$ is not relevant for any resonant motion. While the azimuthal frequency is trivial, the radial and the polar frequencies will have independently varying values along the in-spiral. Resonance occurs when these two frequencies become commensurate with each other, or, in more formal terms, when their ratio is a rational number (Ruangsri and Hughes, 2014):

[TABLE]

During resonance, the apparent effect of the self-force on the orbit is modified, in particular there is a fractional change of order $\mathcal{O}\!\left(\epsilon^{0}\right)$ to the dynamics of the in-spiral – the evolution rate of parameters and the in-spiral rate (Flanagan and Hinderer, 2012). A typical emri passes through several resonance points before the two objects merge. Obviously, this complicates the numerical analysis of the in-spiral, since one has to develop a formalism to account for the occurrence of resonances, and to quantify the effects these have on the orbital dynamics (Flanagan and Hinderer, 2012).

In the current paper, we will first describe the two approaches to resonance modelling that have appeared in the literature, which were suggested by (Hinderer and Flanagan, 2008) and (Gair, Yunes, and Bender, 2012). We then describe what the differences in the model assumptions are, explain why there is a preferred case of physical relevance to resonances in emris, and argue that full knowledge of the self-force is required to robustly identify the most relevant model. We will then extend the results of (Gair, Yunes, and Bender, 2012) to obtain complete results describing the transition through resonance in that model, at leading order in the mass-ratio.

II Instantaneous frequency approach

This method is described in detail in Gair, Yunes, and Bender (2012), although here we will give an outline of the approach to aid further discussions. As explained above, individual orbits of the compact object around the black hole can be thought of as tri-periodic geodesics with fundamental frequencies $\{\omega_{r},\omega_{\theta},\omega_{\phi}\}$. This encourages us to employ the “osculating element” formalism, where the in-spiral is modelled as a sequence of geodesics, with the constants of motion continuously evolving as the orbit approaches the black hole, and furthermore, the effect of the compact body on the trajectory can be treated as a perturbation (Pound and Poisson, 2008). Since these constants of the motion do not change on a given geodesic, at leading order, and only evolve (continuously) as the object spirals into the supermassive black hole, we can argue that, to leading order, this effect is driven by the self-force on the object. Hence, we can write the evolution equation for $J_{\mu}=(\mathcal{E},\mathcal{L}_{z},\mathcal{Q})$ as follows:

[TABLE]

where $\mathcal{F}_{\mu}$ is a function of the coordinates, the generalised coordinates, and the self-force, which describes the rate of change of orbital parameters as the object moves along the inspiral. We argued that the self-force is itself a function of the coordinates and the fundamental frequencies (see Section I.2), hence we can write $\mathcal{F}_{\mu}=\mathcal{F}_{\mu}(\mathbf{q},\bm{\omega})$. Similarly, the action-angle coordinates of the system evolve proportionally to the fundamental frequencies, with a correction owing to the self-force:

[TABLE]

We will not concern ourselves with the precise form of the functions $\mathcal{F}_{\mu}$ and $\mathcal{H}_{\mu}$, but rather will point out that they can be expanded in tri-harmonic Fourier series as discussed above, which is valid without any loss of generality. $\omega_{\mu}$ is the instantaneous value of the frequency at this moment along the trajectory. This is calculated according to eq. (3), and the value is continuously updated to reflect the motion of the body along its inspiral and the effect of the self-force. Without the self-force, the frequencies will, of course, remain constant, since the trajectory would be that of a geodesic around the black hole, unaffected by gravitational wave emission.

Before carrying out the Fourier mode expansion, we recall that a geodesic in Kerr space-time is uniquely specified by three constants. Instead of the orbital parameters $\mathcal{E}$, $\mathcal{L}_{z}$, and $\mathcal{Q}$, we could use the set of fundamental frequencies – the radial frequency, $\omega_{r}$, the polar frequency, $\omega_{\theta}$, and the azimuthal frequency, $\omega_{\phi}$, to specify the geodesics, as these are constant along a given geodesic. Since these can be mapped to the aforementioned integrals of the motion (Schmidt, 2002), we can use these maps to argue that the fundamental frequencies evolve along the trajectory according to a differential equation of the same form as eq. (3):

[TABLE]

The two orbital frequencies relevant to the occurrence of resonances are the radial and the polar (lateral) frequencies. As noted previously, the third frequency $\omega_{\phi}$, that of azimuthal motion, does not play a role in the occurrence of resonances due to the axisymmetric nature of the Kerr space-time. We use this property to write a bi-periodic expansion for the function of the self-force in Fourier modes of the two frequencies $\{\omega_{r},\omega_{\theta}\}$:

[TABLE]

We can regard this equation as an adiabatic evolution equation – at leading order the terms on the right-hand side are constant in time since they are taken along the same geodesic, while the differential on the left-hand side expresses their evolution along the trajectory of the compact object towards the black hole. Of course, the quantities on the right-hand side are instantaneous, and their values are constantly updated through the evolution equations. The sum is over the non-negative integers $l$ and $m$. Since we have written the time-dependence explicitly in the Fourier series, the coefficients of the expansion depend smoothly on the constant parameters of the geodesic. Therefore, we can expand the coefficients in orders of the mass ratio, and neglect all terms sub-leading in $\epsilon$.

Alternatively, instead of using the $\{\omega_{\mu}\}$ as the basis for the Fourier expansion, we could rewrite the evolution equation in terms of the phase angle variables $\{\phi_{\mu}\}$. Phase variables are related to the fundamental frequencies through $\dot{\phi}_{\mu}=\omega_{\mu}$, and since on a geodesic the frequencies remain constant, $\phi_{\mu}=\omega_{\mu}\tau$. Changing variables from frequencies to phase angles in eq. (6), we arrive at a second-order differential equation:

[TABLE]

Geodesic resonances occur when the two frequencies become commensurate, or, more precisely, when their ratio is a rational number (eq. (2)). Therefore, near resonance the term $\left(l_{0}\,\omega_{r}-m_{0}\,\omega_{\theta}\right)$ will approach 0 (for some $l_{0}$ and $m_{0}$ and their integer multiples), and the sum will contain both slowly-oscillating resonance terms, and rapidly oscillating off-resonance terms. Let us now restructure eq. (6) to change variables to $\bar{\omega}=l_{0}\,\omega_{r}-m_{0}\,\omega_{\theta}$, which is small near a resonance. Now the arguments of the oscillating modes on the right-hand side can be written as a linear combination of $\bar{\omega}$ and $\{\omega_{\mu}\}$ (we can use either of them, but choose $\omega_{r}$ for definitiveness):

[TABLE]

Similarly, we can rewrite the phase equation in terms of $\bar{\phi}$ and $\phi_{r}$, where $\bar{\phi}=l_{0}\,\phi_{r}+m_{0}\,\phi_{\theta}$ vanishes on resonance:

[TABLE]

Let us change variables to absorb the factor of the mass ratio appearing on the right-hand sides, namely:

[TABLE]

and

[TABLE]

Using this substitution, eqs. (8) and (9) become

[TABLE]

For any non-zero value of the index $m$, corrections to the frequency or phase equations are of the form $\cos(\gamma x+\delta)$, which oscillate rapidly owing to the $\gamma\sim\epsilon^{-\nicefrac{{1}}{{2}}}\sim 10^{3}$ factor both on resonance and away from it (and therefore will average out to zero over the timescales under consideration — this is illustrated in Figure 1). Hence, the correction to the first order solution of the phase equation is of the form $(1/\gamma^{2})\cos(\gamma x+\delta)$ and that of the frequency equation is $(1/\gamma)\sin(\gamma x+\delta)$. These are terms of order $\epsilon$ and $\epsilon^{\nicefrac{{1}}{{2}}}$, respectively. Hence, we can perform an expansion to establish that they are corrections of higher order in the mass ratio and hence can be ignored at the same level as the other higher order in mass ratio terms that we have omitted. We can therefore consider only terms with $m=0$, and relabel the remaining modes by $n$, obtaining

[TABLE]

We combine the two oscillating modes together, and this introduces a phase to the oscillating term. This could sometimes be set to 0 by suitable choice of initial coordinates, however it is not always possible, since the resonance phase depends on the prior evolution of the binary system. To preserve generality of the problem, we will include it in our equations. The presence of phase does not complicate the method of solution, only brings in additional terms, as we will demonstrate later in Section IV.

[TABLE]

where

[TABLE]

Finally, we split the first term from the sum, $n=0$, which is constant (with respect to time or phase, respectively):

[TABLE]

and make a suitable, final change of variables:

[TABLE]

and

[TABLE]

to bring the equations to their most simplified form:

[TABLE]

where $\kappa_{n}\!=(\mathcal{A}_{n}/\mathcal{A}_{0})\sec(\delta\chi_{0})$ and $\lambda_{n}\!=(\mathcal{B}_{n}/\mathcal{B}_{0})\sec(\delta\psi_{0})$ are in principle known constants, dependent on the orbital parameters at resonance, prime denotes a derivative with respect to the rescaled time $x$ and the quantities $\{\delta\chi_{i}\}$ and $\{\delta\phi_{j}\}$ are constants. The sums over $n$ run from 1 to $\infty$, however we will later consider finite number of modes on the right-hand side. In order to achieve convergence when solving these equations, we need to make the assumption that the series $\{\kappa_{n}\}$ and $\{\lambda_{n}\}$ follow a power law of the form $n^{-\alpha}$, with $\alpha\gtrsim 2$, this is discussed in detail in Section IV.1. In this paper, we are attempting to find solutions to eqs. (19a) and (19b) in increasing order of complexity. It is worth noting that while one of them is a first-order equation, the other is of second order. They derive from equivalent representations of the geodesics, however the first and second order equations are no longer equivalent, because they make different assumptions about the expansion of the self-force.

Before proposing the respective solutions, we will focus on presenting an alternative approach to describing the behaviour of an emri near resonance.

III Two-timescale approach

An alternative approach to describing the transition of emris through resonances was first suggested by (Hinderer and Flanagan, 2008), and is often called the “two-timescale method” since it differentiates between processes which occur on short (orbital) timescales, and those happening on long (orbit evolution) timescales. A good overview is offered by Kevorkian (1987).

Here we will present this method as described by van de Meent (2014), while it is also advisable that the reader consults Berry (2013), and the original article by Hinderer and Flanagan (2008). Again, we consider the motion of a compact object spiralling towards a supermassive black hole, and describe the evolution in terms of the perturbative parameter, the mass ratio $\epsilon$. At leading order, we recover the geodesic equations of motion, as before. At next order, we need to introduce corrections due to the self-force acting on the compact object, in analogy with eq. (3) in the instantaneous frequency approach:

[TABLE]

While the equations in Berry (2013) are presented in terms of the Mino time $\lambda$ (Mino, 2003), we can map $\lambda$ to the proper time $\tau$ through a function involving the coordinates of the compact object. We can absorb this dependence into the functions $\mathcal{I}_{\mu}$ on right-hand side of eq. (20) and work in terms of the proper time. As before, we can exchange the integrals of the motion $\mathbf{J}$ in favour of the fundamental frequencies $\bm{\omega}$, since mappings between these exist (Schmidt, 2002), albeit not invertable over the entire domain. Further we can implicitly write the self-force as a function of the coordinates and their corresponding frequencies (cf. eq. (5)) to obtain:

[TABLE]

Similarly, the action-angle variables evolve as per eq. (4):

[TABLE]

We construct the resonant frequency $\bar{\nu}=l_{0}\,\omega_{r}-m_{0}\,\omega_{\theta}$, and the corresponding action-angle variable combination $\bar{q}=l_{0}\,q_{r}-m_{0}\,q_{\theta}$. $\bar{\nu}$ is the analogous variable to $\bar{\omega}$ in Section II, but we use a different notation here to differentiate between the two methods. Again, we ignore the third fundamental frequency $\omega_{\phi}$, since the motion is axisymmetric. Let $\bar{J}$ denote the action conjugate to $\bar{q}$, which can be mapped to $\bar{\nu}$. Further, let us rescale the proper time to write $x_{\nu}=\sqrt{\epsilon}\,\tau$. We now expand the frequency and coordinate variables in powers of $\epsilon$:

[TABLE]

Plugging these power laws into the (averaged) evolution equations (21) and (22) yields the resonance equations:

[TABLE]

This relays the leading terms in the expansions (23). Further members of the series $\{\bar{\nu}_{n}\}$ and $\{\bar{q}_{n}\}$ can be found iteratively by substituting the results from eqs. (24) back into the resonance evolution equations.

In Flanagan and Hinderer (2012), eq. (23a) was further simplified by expanding $\bar{\nu}_{1}$ as a series in $\sqrt{\epsilon}\,x_{\nu}$. The leading term is then $\nu_{1}(x_{\nu}-x_{0})$ for a constant $\nu_{1}$. Proceeding with the two timescale expansion then yields a result accurate to first order in the resonant modification, $\{\lambda_{n}\}$, which agrees with the result obtained from the asymptotic Fourier series approach in Gair, Yunes, and Bender (2012). Corrections at higher order in $\lambda_{n}$, but leading order in mass ratio, cannot be found in this way, since the assumption of a linearly evolving frequency is only valid close to the resonance point. The two-timescale approach can be extended to higher order by using a less restrictive ansatz for $\bar{\nu}_{1}$, which is given implicitly by eqs (23) and (25) in van de Meent (2014). Explicit solutions to higher order in the resonant modification at leading order in mass ratio have not been computed in the two-timescale approach. However, we anticipate that these will agree with the phase version of the asymptotic Fourier series approach. The purpose of this paper is to extend previous calculations to higher order in resonant modification and to include multiple modes in the resonance. To do this we will use the asymptotic Fourier series approach and this will be the focus for the remainder of this paper.

IV Solutions of the resonance equations

We are now in a position to solve the instantaneous frequency approach equations that we derived earlier in Section II. Evolution of the action-angle coordinates according to the instantaneous frequency of the system yields the 2 non-linear differential eqs. (19a) and (19b). Due to the constant additive term on their right-hand sides, these equations cannot be solved exactly, and therefore we need to come up with an alternative approach. In this Section we propose solutions to these equations in an increasing order of complexity. Solutions for the single-mode versions of eqs. (19) are provided in Gair, Yunes, and Bender (2012), we will start by re-deriving those to depict how the order-by-order solution algorithm works. Following this, we relax the assumptions which have been considered in Gair, Yunes, and Bender (2012) in order to build up solutions to the generalised source equations.

IV.1 Regime and domain of the variables

Before attempting to solve the resonance equations, we need to specify their domains and any special limits we are taking. The time parameters $x_{\{\omega,\,\phi,\,\nu\}}$ range from $-\infty$ to $+\infty$, however we are interested in the vicinity of $x=0$ as this is when resonance occurs, hence where we have specified the initial conditions. In this notation, $-\infty$ corresponds to times long before resonance, and $+\infty$ to times when the resonance has had an effect on the orbit. Therefore, in solving these equations, we will consider definite integrals of the form $\int_{0}^{x}\mathop{}\!\mathrm{d}\bar{x}$, where $x$ can take large positive or negative values. Since the equations are not coupled, and each rescaled time variable appears in the differential equation for its corresponding quantity, we will drop the subscripts on $x$, and will keep track of them by the subject of the equation we are solving.

More attention needs to be devoted to the fractional resonance modifications, $\kappa_{i}$ and $\lambda_{j}$, which are themselves functions of the orbital parameters. As discussed in Gair, Yunes, and Bender (2012), if these modifications are of order unity, there can be a qualitative difference in the on-resonance behaviour, with the possibility of sustained resonances. These were explored in more detail in van de Meent (2014). A complete study of the size of on-resonance flux modifications has not yet been carried out. Some insight on the likely size of these corrections can be obtained by studying solutions to the Teukolsky equation, which provide a decomposition of the off-resonance energy flux into different Fourier modes. Such results can be found in Drasco and Hughes (2006). Section I C of that paper indicates that the size of the terms in this decomposition initially follows a power law of the form $\sim n^{-\alpha}$ before eventually falling off exponentially (see Figure 2 and Table 1 for an illustration different rates of fall off). From Fig. 3 of that paper it appears that $\alpha\gtrsim 2$ for all $n$, but as a lower limit in our calculations we will use the value of $\alpha=5/2$. Drasco and Hughes (2006) also indicate that the size of the first term in the series is less than unity, and in fact of order $\kappa_{1}\sim\lambda_{1}\sim 0.1$. This implies that the resonance is likely to be “weak”, but the Teukolsky results do not directly give the size of the on-resonance flux change. The size of the resonance effect was assessed more carefully in flanaganhughesruangsri2012, who computed the flux change for a number of different resonances for a few systems and found that the resonant flux change was only $\sim 1\%$ in all cases. This also suggests that $\kappa_{1}$ and $\lambda_{1}$ are small, typically $\sim 0.01$. Based on this evidence it seems unlikely that strong resonances will occur in practice and so in this article we focus on the weak resonance case with $\kappa_{1}\ll 1$, which will allow us to perform an expansion in $\kappa_{1}$. However, we should note that flanaganhughesruangsri2012 did not thoroughly explore the resonance parameter space and so it is still possible that there are regions of parameter space where strong sustained resonances may exist. If this is possible, it would be very interesting, as discussed in detail in van de Meent (2014).

Regarding the other relevant parameters, the mass-ratio $\epsilon$ is typically in the range $\sim{10^{-6}}-\sim{10^{-3}}$. The quantity $\epsilon^{1/2}$ is therefore $\mathcal{O}\!\left(10^{-3}-10^{-1.5}\right)$. The phases on resonance, $\delta\chi_{i}$ and $\delta\psi_{j}$, are constants in the range $[0,2\pi)$, which depend only on the coefficients of the Fourier series in eqs. (6) and (7), and the index $n$.

IV.2 Frequency resonance equation

Consider the generalised frequency resonance equation:

[TABLE]

where prime denotes differentiation with respect to the rescaled time parameter $x_{\omega}$. Before attempting to solve this equation in the general case, we consider only the single oscillatory mode ($n=1$, besides the $n=0$ mode) on the right-hand side, and as agreed, drop the subscript on $x$. The resulting problem is not much simpler to solve, however it will provide useful insight about the method:

[TABLE]

where $\kappa\equiv\kappa_{0}$. Inspired by the numerical solution of eq. (25), shown in Fig. (3), we split the solution in two parts: the first for the region where the oscillations are still growing, and the second where the oscillations gradually decay. After we establish appropriate solutions to eq. (25) in each of these regions, we will come up with a way to match them on the boundary region.

Let us initially look for a solution in the $\kappa|x|\ll 1$ region. That is, find an expansion of the frequency $\omega$ as a function of $x$ in the region immediately around $x=0$, in powers of the fractional resonance modification $\kappa$:

[TABLE]

We remind ourselves that the parameter $\kappa\ll 1$, as per Section IV.1. In order to find explicit forms of the functions $\{\varpi_{n}\}$, we substitute the expansion (26) in eq. (25) and proceed to match terms on both sides. The results we now derive are plotted in Figure 4. At lowest order, we disregard all terms of $\mathcal{O}\!\left(\kappa\right)$ and higher:

[TABLE]

where it is important to note that the integration has been carried out between 0 and $x$. Further terms in the series (26) will be reduced by one or more factors of $\kappa\ll 1$, hence we can treat them as a perturbation to $\varpi_{0}(x)$:

[TABLE]

At the next higher order we find the equation for $\varpi_{1}(x)$:

[TABLE]

integrating the right-hand side yields the solution:

[TABLE]

where $C(\bullet)$ is the Fresnel cosine integral (see Appendix D and Figure 25 for details). At $\mathcal{O}\!\left(\kappa^{2}\right)$, we find:

[TABLE]

Solving this differential equation (which also involves a simple integration by parts) provides the solution:

[TABLE]

It is obvious that with increasing order the functional form of $\varpi_{n}(x)$ become more and more involved. It will become apparent later on why we do not need solutions at $\mathcal{O}\!\left(\kappa^{4}\right)$ and beyond. However, at third order, we find the following differential equation for the functional $\varpi_{3}(x)$.

[TABLE]

We divide $\varpi_{3}$ into two parts, $\varpi_{3}(x)=\varpi_{3a}(x)+\varpi_{3b}(x)$, defined via the two differential equations:

[TABLE]

Unfortunately, the term $\varpi_{3a}(x)$ above is the only integral in our analysis which cannot be performed analytically. However, we can proceed by expanding the Fresnel term on the right hand side for small $x$ and integrating the resulting terms. We then define $\varpi_{3a_{\rm approx}}(x)$ to be the terms arising from this integration, given analytically by eq. (34). This has the property $\varpi_{3a}^{\prime}(x)-\varpi_{3a_{\rm approx}}^{\prime}(x)\sim\mathcal{O}\!\left(x^{9}\right)$ and hence $\varpi_{3a}(x)-\varpi_{3a_{\rm approx}}(x)\sim\mathcal{O}\!\left(x^{10}\right)$. Finally, we need to verify that this is a valid approximation; we can do this by tabulating the function

[TABLE]

for values of $x$ between $-1$ and 1, and comparing this with our (approximate) analytical result for the integral, eq. (34). The integrand, both functions, and the residuals of the approximation are shown in Fig. 5.

The term in eq. (32b) can be integrated analytically, and the result is given by expression (35). Gathering the above results together, we can write down an approximate solution to the Frequency equation, which is given in eq. (36). We note that this is valid for $\kappa|x|\ll 1$ up to $\mathcal{O}\!\left(\kappa^{2}\right)$, while the result for $\varpi_{3}(x)$ is only valid when $|x|\ll 1$. Later on, as the need arises, we will find a solution for $\varpi_{3}(x)$ which is valid in other regimes.

[TABLE]

[TABLE]

Before we continue with our analysis, let us consider the asymptotic value of this solution in the limit of it’s validity, i.e. for $1\ll|x|\ll\kappa^{-1}$. The non-oscillating part of solution (36) can be written as

[TABLE]

We see that any constant term in the solution for $\omega(x)$ is at least $\mathcal{O}\!\left(\kappa\right)$, rather than $\mathcal{O}\!\left(\kappa^{0}\right)$. To confirm the validity of solution (36) we compute the quantities

[TABLE]

where $\omega^{(n)}(x)$ is eq. (36) restricted only to order $n$ (e.g. $\omega^{(1)}(x)=\varpi_{0}(x)+\kappa\varpi_{1}(x)$). If our solution is correct, then we would find that the fractional error, computed by the above formula is $\mathcal{O}\!\left(\kappa^{n+1}\right)$. This scaling is evident in the plot of these functions, shown in Figure 6.

Let us now seek a solution in the region of decaying oscillations, $|x|\gg 1$. At lowest order, the solution is

[TABLE]

At higher order, we attempt an ansatz, composed of a constant, a linear term, and a decaying oscillating term. We allow for three undetermined constants, $a_{1},b_{1}$ and $c_{1}$, as suggested by Gair, Yunes, and Bender (2012):

[TABLE]

where the subscript “1” on the slope correction term indicates that it pertains to the solution $\omega_{1}(x)$, as the value of $n$ will change for higher-order solutions. When we substitute this ansatz into the model eq. (25):

[TABLE]

we ignore $\mathcal{O}\!\left(x^{-1}\right)$ terms and beyond, and use our knowledge of the behaviour of each separate term to reason about the values of the undetermined constants in eq. (42) below. Using standard trigonometric identities, we can expand the second term on the right-hand side, as shown explicitly in eq. (43) below.

[TABLE]

At this point we introduce a shorthand notation for the oscillating terms. Please refer to Appendix E (and the plots in Figure 26) for details.

[TABLE]

Using these, we cast our last equation (43) in the form:

[TABLE]

Comparing terms on both sides of eq. (45), we notice that in order for this identity to hold, it is necessary that $c_{1}\sim\mathcal{O}\!\left(\kappa\right)$. We remind ourselves of the assumption that $\kappa\ll 1$ which justifies the small angle approximation we used above. Since the range of the sine function is between $-1$ and 1, and $c_{1}\sim\mathcal{O}\!\left(\kappa\right)\ll 1$, we can safely assume that $c_{1}\sin\!\left[ax+\left(1+b\right)x^{2}\vphantom{()^{\!\!\!2\,}}\right]\ll 1$ for all $x$ and expand both trigonometric terms in eq. (45):

[TABLE]

As we only retain the leading order terms in these expansions, we must disregard any $\mathcal{O}\!\left(\kappa^{3}\right)$ or higher terms in our further calculations. Since we exclude them at this stage, such terms can no longer appear in the final solution for $\omega_{1}(x)$. We match the coefficients of the two $\mathcal{O}\!\left(\kappa\right)$ terms on either side of the relation (46): $2\left(1+b\right)\,c_{1}\sim\kappa$. Therefore:

[TABLE]

By using the double-angle trigonometric identities on the remaining term on the right-hand side of (46), we find

[TABLE]

Combining this with result (47), we establish that $b_{1}\!\sim\mathcal{O}\!\left(\kappa^{2}\right)$, i.e., 0 to linear order. This result follows solely from the linear independence of the $\{C_{n}\}$ terms. While they are not orthogonal (this can be easily verified by analytical integration), they are linearly independent, and therefore the coefficients on both sides need to match for eq. (48) to hold true. Similarly, we obtain the revised result $c_{1}\sim\kappa/2+\mathcal{O}\!\left(\kappa^{3}\right)$. The additive constant $a\equiv\omega(0)$ can be determined by comparing the limit $|x|\gg 1$ of solution (36) to the ansatz (40) in the regime $\kappa|x|\ll 1$:

[TABLE]

where on each side there are terms which we have ignored in the respective analysis. We find that a=\pm\big{(}\sqrt{2\pi}/4\big{)}\kappa to linear order, and we can write down the leading-order solution (50) to the frequency equation (for a single oscillatory mode and no initial phase). Note that here, and for the rest of this article, whenever the symbols $\nicefrac{{\pm}}{{\mp}}$ are used, the upper sign is valid for $x\geq 0$, and the lower for $x\leq 0$. The function in question is thus defined piecewise, but we will use this notation throughout for brevity.

[TABLE]

Once we have established the suitability of this ansatz, we propose a further trial solution by adding an extra term to eq. (40). It accounts for the terms in eq. (48) we ignored when considering only first-order corrections:

[TABLE]

The full ansatz can now be written as:

[TABLE]

We substitute this into the source equation (25) and carry out a similar procedure to establish the new values of the constants, as detailed below. On the left-hand side, we ignore terms of $\mathcal{O}\!\left(x^{-1}\right)$ in the derivative:

[TABLE]

while on the right-hand side, we need to unpack the argument of the cosine function. Previously we established that $c_{1}\sim\mathcal{O}\!\left(\kappa\right)$, and since we expect that the new terms originating from $\omega_{2}(x)$ will match the remaining term in eq. (48), we expect $c_{2}\sim\mathcal{O}\!\left(\kappa^{2}\right)$. This justifies the further use of the small-angle approximation for the latter two terms. We point out that we need to ignore non-constant terms of $\mathcal{O}\!\left(\kappa^{3}\right)$ in this expansion, as they will not be matched by a suitable term in the ansatz.

[TABLE]

Matching the coefficients of terms in eqs. (53) and (54), we obtain the coupled equations:

[TABLE]

These imply that $c_{1}\sim\kappa/2+\mathcal{O}\!\left(\kappa^{3}\right)$ and $b_{2}\!\sim-\kappa^{2}/4+\mathcal{O}\!\left(\kappa^{4}\right)$. The last equation above determines the value of the newly introduced constant $c_{2}$:

[TABLE]

To find the constant $a$ at second order, we need to match this solution to the $\mathcal{O}\!\left(\kappa^{2}\right)$ part of solution (36) in the buffer zone where both solutions are valid, as per eq. (57). Expanding (36) for $x\gg 1$ and the above ansatz for $\kappa|x|\ll 1$ (both to order $\mathcal{O}\!\left(x^{-1}\right)$ inclusive), we find that all functional term match on both sides. The constant term is the value of $a$ at second order, \big{(}\!\mp\sqrt{\pi}/16\big{)}\kappa^{2}.

Combining these results together, we establish the improved solution $\omega_{2}(x)$, given by eq. (58).

[TABLE]

[TABLE]

At the next order, the task of finding appropriate terms becomes more complicated, since we need to balance the terms in eq. (34). At this order there is a logarithmic trend, as well as an oscillatory part, hence we propose the following ansatz which reflects the nature of the solution:

[TABLE]

The details of calculating the expressions for these constants are burdensome, hence we have left them in Appendix I. The important message here is that we have still managed to find a set of terms which satisfy the governing equation (25), without conflicting the previous solutions (i.e. without producing new terms at $\mathcal{O}\!\left(\kappa,x^{-1}\right)$).

To determine the new value of the constant $a$, we again need to turn to our solution near the origin. In particular, here we need an expansion of $\varpi_{3}(x)$ for large $x$. However, we used a small-$x$ approximation to obtain $\varpi_{3a_{\rm asymp}}(x)$ earlier, which is not valid in the buffer zone where both solutions (36) and (59) are valid. We need to go back to the original problem (32a), and instead of expanding in the small-$x$ regime, we expand for $|x|\gg 1$ (which is valid when $\kappa|x|\ll 1$). To do this, we expand the Fresnel term on the right hand side of eq. (32a) for large $|x|$, truncate the series at $\mathcal{O}\!\left(x^{-5}\right)$ and then compute the antiderivative of this series, the function $\varpi_{3a_{\rm asymp}}$ given in eq. (62). This procedure ensures that $\varpi_{3a}(x)-\varpi_{3a_{\rm asymp}}(x)\sim\alpha+\mathcal{O}\!\left(x^{-4}\right)$. The constant $\alpha$ is the asymptotic value of $\varpi_{3a}(x)-\varpi_{3a_{\rm asymp}}(x)$, which may be found by integrating $\varpi_{3a}(x)-\varpi_{3a_{\rm asymp}}(x)$ from [math] to infinity

[TABLE]

The quantity $x_{c}$ is arbitrary, but by choosing a value $x_{c}\sim 1/k_{0}\sim 10$ we may bound the value of the second integral by $10^{-5}$. Numerical evaluation of the terms on the first line of eq. (60) then give $\alpha=\mp\,0.00991$. We deduce that, asymptotically, we may write $\varpi_{3}(x)=\alpha+\varpi_{3a_{\rm asymp}}(x)+\varpi_{3b}(x)$, with higher-order corrections which scale as $x^{-4}$.

Combining these resulting terms with the ones coming from eq. (35), we arrive at the left-hand side of eq. (63). The right-hand side of eq. (63) is comprised of the terms in an expansion of the ansatz (59) for small $\kappa|x|$.

The terms match nicely as we would expect, apart from the coefficients of \sin\!\big{(}x^{2}\big{)}/x and \cos\!\big{(}x^{2}\big{)}/x. The reason for this is that to determine these terms accurately we need to extend the asymptotic solution to the next order in $1/x$. At that order we get a term of the form $C_{1}(x)/x^{3}$ in the solution, and it is straightforward to deduce that, at leading order, the coefficient in front of this term is $-\,\kappa/4$. This additional term contributes terms at order $\kappa^{3}$ on the right hand side of eq. (63) and those terms are \pi\cos\!\big{(}x^{2}\big{)}/(64x) and -\sin\!\big{(}x^{2}\big{)}/(16x), which are exactly what is needed to match the terms from the $\kappa|x|\ll 1$ solution on the left hand side of eq. (63).

We can now read off the third-order piece of the coefficient $a$ from eq. (63)

[TABLE]

The third-order solution $\omega_{3}(x)$ is given in eq. (64). This is the most accurate solution to the single-mode, zero-phase frequency equation that we will present. The lack of necessity for higher order iterations will be explained when discussing the multi-mode frequency equation.

Once we have obtained these solutions, it is worthwhile to investigate their properties visually. To confirm their validity, Figure 7 demonstrates that these solutions are correct at their respective order by computing the quantities $\varkappa_{n}(x)$ defined in eq. (38), and showing graphically that they have the expected magnitudes. Please note we have shown only the region $x>0$, with the region $x<0$ being a mirror image of this one.

Figure 8 shows how the intercept $a$, the slope $(1-b)$, and the coefficients of the oscillating modes $\{c_{n}\}$ and $\{d_{n}\}$ vary with the parameter $\kappa$.

[TABLE]

[TABLE]

[TABLE]

The solution can be made more accurate, by adding terms at higher orders, of the form $c_{n}S_{n}(x)/x$ and $d_{j}S_{n}(x)/x^{2}$ to the ansatz (52). Below we see the general case of having $m$ oscillatory terms in the trial solution $\omega_{m}(x)$, where $m$ could be either a finite number or infinity:

[TABLE]

We shall discuss later (after we solve the multimode version of the phase equation) why there is no practical need for solutions beyond third order. If we seek to find lower-order corrections in $\mathcal{O}\!\left(x^{-p}\right)$, we add terms of the form:

[TABLE]

However, instead of following through with this procedure, we look at extensions of eq. (25) to find the general solution of eq. (19a). Initially, we allow for multiple oscillation modes, each with a different amplitude (parameter $\kappa_{i}$), and following this, we include a non-vanishing phase in our calculations.

In the following we will use bold quantities, e.g., $\bm{\omega}(x)$, to indicate values of the coefficients in the multiple mode case. We note that these are not vectors, but the scalars that appear in the multiple mode solution. To solve for the evolution of the frequency resonances,

[TABLE]

we can again look at the problem from two scales: small and large $|x|$. Near the origin, we can write down a series expansion similar to eq. (26):

[TABLE]

Similarly to before, we can derive solutions for $n=0,1$:

[TABLE]

Further solutions are not going to be pursued, as higher orders of $\kappa_{n}$ (for small $n$) would be of the same size as $\kappa_{n}$ for some larger $n$. Hence, $\bm{\omega}(x)$ near the origin can be most completely written down as

[TABLE]

For $|x|\gg 1$, we propose an extended ansatz of the form:

[TABLE]

with derivative

[TABLE]

We substitute this into the governing eq. (67) and proceed to match the relevant coefficients on both sides of this relation. It is important to explain here that the leading order dependence of each coefficient $\bm{c}_{n}$ will be the same as if considering its corresponding mode on its own (see Appendix II for derivation), with higher-order corrections due to cross-terms emerging from the cosine term (see Appendix III), and therefore $\bm{c}_{n}\sim\mathcal{O}\!\left(\kappa_{n}/n\right)\sim\mathcal{O}\!\left(\kappa_{1}/n^{\alpha+1}\right)$. Thus, each coefficient $\bm{c}_{n}$ is smaller than the previous one in the series, and we can expand the trigonometric terms on the right-hand side, to obtain the form (73). Evidently, to be able to match the constant term $\bm{b}$ on the left-hand side, we need to find a constant term on the right-hand side too. Such a constant could only arise from the second term on the right-hand side, in the case $n=m$, and yields eq. (74). The remaining terms must combine to form an expression for $\bm{c}_{n}$, and we match them using the linear independency of the functions $\{C_{n}(x)\}$, yielding eq. (75). On the right-hand side of this equation, we omitted a pre-factor of $(1+\bm{b})^{-1}\sim 1+\mathcal{O}\!\left(\kappa_{1}^{2}\right)$, since we know that $\bm{b}\sim\mathcal{O}\!\left(\kappa_{1}^{2}\right)$.

$\displaystyle\bm{b}+2(1+\bm{b})\sum_{n}n\,\bm{c}_{n}C_{n}(x)\sim\sum_{n}\kappa_{n}\,C_{n}(x)-\frac{1}{2}\sum_{n,\,m}n\,\kappa_{n}\,\bm{c}_{m}C_{n-m}(x)+\frac{1}{2}\sum_{n,\,m}n\,\kappa_{n}\,\bm{c}_{m}C_{n+m}(x)$

(73)

$\displaystyle\bm{b}\sim-\,\frac{1}{2}\sum_{n}n\,\kappa_{n}\,\bm{c}_{n}$

(74)

$\displaystyle\bm{c}_{n}\sim\frac{\kappa_{n}}{2n}-\frac{1}{4n}\sum_{m,\,p}m\,\kappa_{m}\,\bm{c}_{p}\,\delta_{|m-p|,\,n}+\frac{1}{4n}\sum_{m,\,p}m\,\kappa_{m}\,\bm{c}_{p}\,\delta_{(m+p),\,n}$

(75)

We now proceed to demonstrate that the second and third terms in the relation (75) are sub-leading compared to the first term, and therefore can be ignored in the final expression. To do this, we remind ourselves of our assumption that $\kappa_{n}\sim\kappa_{1}/n^{\alpha}$, with $\alpha\gtrsim 2$ and substitute this into our expression for $\bm{c}_{n}$. From the first term we can see that $\bm{c}_{n}\sim\mathcal{O}\!\left(\kappa_{n}/n\right)$. Using this, the formula for $\bm{c}_{n}$ can be written down solely in terms of $\kappa_{1}$. After applying the action of the Krönecker delta operators, the second and third terms of eq. (75) become

[TABLE]

and

[TABLE]

We need to inspect and bound each of these terms, in order to show they are sub-leading compared to the first term in eq. (75), $\kappa_{1}/n^{\alpha+1}$. In eq. (76), we split off the $p=1$ term from the sum, and bound the remaining terms from above by an integral over the index $p$:

[TABLE]

We re-scale by writing $p=n\tilde{p}$, and the integral becomes:

[TABLE]

We split the integration region in two intervals: from $(1/n)$ to $1$, and from $1$ to $\infty$. The first of the resulting integrals can be bounded as follows

[TABLE]

while the second one can be solved and the solution expressed in terms of the hypergeometric function ${}_{2}F_{1}$, independently of the integration variable $\tilde{p}$. We find that the term $\Xi_{n}(\kappa_{1})$ is bounded above by ($\xi_{i}$ are constants):

[TABLE]

Considering the term in eq. (77), we can find a bound using a similar argument, albeit the reasoning is slightly longer, as the sum has definite bounds on both ends, and the integrand is not monotonic for the interval of integration. We strip off the first and last terms of the sum: $(n-1)^{-\alpha+1}$ and $(n+1)^{-\alpha-1}$, respectively. The remainder can be bound from above by an integral,

and a suitable substitution can be applied to obtain:

[TABLE]

We proceed to split the integral for some $n>1$, say 2:

[TABLE]

Each of the two integrals can be bounded as follows:

[TABLE]

Evaluating these and putting everything back together, we conclude that a bound on the term (77) is given by

[TABLE]

where $\upsilon_{i}$ are constants. We have found that the latter two terms in eq. (75) are of the order of $\kappa_{1}^{2}/n^{\alpha}$. We still need to determine how this compares to the magnitude of the first term in this result, namely $\kappa_{1}/n^{\alpha+1}$:

[TABLE]

We can therefore say that for values of $n$ smaller than $\kappa_{1}^{-1}$ the second and third terms are much smaller than the first and can therefore be ignored. For values of $n$ above this limit all terms are of comparable magnitude. However, this $\bar{n}$ can be related to the order of the parameters $\{\kappa_{i}\}$:

[TABLE]

and so ignoring these higher terms leads to an error smaller than $\kappa_{1}^{3}$. Therefore, we have discovered that we are free to ignore the terms discussed above, and in general ignore the interference between the magnitude of $\kappa_{1}$ and the large values of the index $n$, as long as we do not concern ourselves with corrections beyond third order. We may write:

[TABLE]

As before, the constant $\bm{a}$ in eq. (71) can be determined if we compare our multimode ansatz to solution (70) in the buffer zone where both solutions are valid. This means expanding (70) for large $x$, and comparing with eq. (71) (see eq. (91). This suggests that the value of $\bm{a}$ is:

[TABLE]

Our final solution with this ansatz is given by eq. (92).

To test the validity of this solution, we devise a test similar to eq. (38), but include the presence of multiple modes in the definition of the quantifier:

[TABLE]

Here, $\bm{\omega}(x)$ is given by eq. (92), $\bm{\omega}_{\text{num}}$ is the respective numerical solution, and the subscript $m$ signifies that all equations and solutions have been limited to the first $m$ terms only. In Figures 10 and 27 we demonstrate graphically the validity of our solution for $m=\{2,5,10,20\}$ and in the intervals $|x|\leq 1$ and $1\leq x\leq 20$, respectively.

We could try to extend this approach further, by naively including $S_{2n}(x)$ terms in the ansatz. These terms, however, are redundant, as for $n$ ranging from 1 to $+\infty$, the set $\{2n\}$ is a subset of $\{n\}$, and therefore the latter terms in the sum can be absorbed into the former one. Thus, eq. (71) is the only ansatz we can use for the equation with an infinite number of modes. Finding higher-order solutions is also complicated by the fact that taking powers of the infinite sum introduces an increasing number of cross-terms.

[TABLE]

[TABLE]

Instead, we return to the single-mode version of the frequency equation and include a non-zero phase $\delta\chi$:

[TABLE]

To find the solution for small $|x|$, we can expand the cosine term on the right-hand side of the equation, and proceed in an analogous way to the zero-phase case:

[TABLE]

The zeroth-order solution is $\varpi_{0}(x)=x$, while higher-order solutions will differ from our previous results. The first-order function can be found from

[TABLE]

After integrating, we arrive at the following result:

[TABLE]

The next-order function obeys the differential equation

[TABLE]

and its solution is given in eq. (108). The third-order function, $\varpi_{3}(x)$ deserves a little more attention, as was the case for the zero-phase equation. Part of the derivative is again not integrable analytically, but could be approximated for the cases $|x|\ll 1$ and $|x|\gg 1$, so we split the solution in two parts, given in eqs. (98) and (99).

[TABLE]

and

[TABLE]

To find a solution for $\varpi_{3a}(x)$, we proceed in an analogous way to the $\delta\chi=0$ case, and take the appropriate limits of interest of the Fresnel function. For $|x|\ll 1$, the expansion is a polynomial, finding an analytical approximation is straightforward, and the result is given by eq. (109). Hence the third-order solution near the origin is

[TABLE]

We also calculate the asymptotic expression for $\varpi_{3a}(x)$, by expanding the Fresnel function for $|x|\gg 1$, the solution is given by eq. (110). The additive constant accounting for the value of the integral at the origin is

[TABLE]

The asymptotic expression for $\varpi_{3}(x)$ can be written as

[TABLE]

The validity of these solutions can be demonstrated in an analogous way to the zero-phase case, by computing:

[TABLE]

where $\omega^{(n)}(x)$ is the solution to eq. (93) for $|x|\ll 1$, and $\omega_{\text{num}}$ is the respective numerical solution. Plots of these functions for $n=0,\dots,3$ are given in Figure 13.

No further calculations are needed to find the constants $\{b_{i},c_{1},c_{2},c_{3},d_{1},d_{2},l\}$, since we can use modified ansätze which include a phase term in their oscillating factors:

[TABLE]

To establish the constants $\{a_{i}\}$, we match the ansatz against the new (asymptotic) solution:

[TABLE]

where the last term is the combination of terms (102) and $\alpha$ is given by eq. (101). The three separate orders of the constant $a$ as a function of the phase difference $\delta\chi$ are:

[TABLE]

The constant $a$ itself is given by

[TABLE]

Plots of these functions are given in Figure 12.

[TABLE]

[TABLE]

Solutions $\varpi_{0}(x)$, $\varpi_{1}(x)$, $\varpi_{2}(x)$ and $\varpi_{3}(x)$ are plotted in Fig. (9), while a more rigorous derivation of the principle we use here is presented in Appendix IV. This paves the way for finding a solution to the general equation with an infinite number of modes.

Before we proceed with this, we check graphically the validity of solution (112). Fig. 14 depicts the fractional errors $\varkappa_{4}(x)$ for the usual values of $\kappa$, and for 4 different values of the initial phase offset.

Finally, we consider the multi-mode frequency equation with a non-vanishing initial phase:

[TABLE]

where the array of constants $\{\delta\chi_{n}\}$ is known and can in principle be computed from the constants of motion through eq. (16). Before proceeding, we elect to write the equation out in the following more revealing way

[TABLE]

For small $|x|$, we find consecutively:

[TABLE]

and

[TABLE]

As before, higher-order solutions are not necessary, so we can write the solution for $|x|\ll 1$ as

[TABLE]

Next we choose an elaborate ansatz of the form:

[TABLE]

where we note that we neglect terms of $\mathcal{O}\!\left(x^{-2}\right)$, which correspond to oscillations which die out quickly with $|x|$. The derivative of this ansatz is given by

[TABLE]

This needs to match the right-hand side of eq. (19a). We perform analysis similar to our earlier discussion on the zero-phase multi-mode version to find the solution (121). Note that the slope correction is still zero at leading order, and we discover that the constants $\{\delta\chi_{n}\}$ only appear in the argument of the oscillation modes. The constant $\bm{a}$ is found by comparing the above ansatz with the $|x|\gg 1$ limit of solution (116).

[TABLE]

Validating these solutions can be done in the same manner as earlier, by computing the quantities analogous to eq. (38) and eq. (90), but this time incorporating the phase variable into the checks. These are depicted in Figures 15 and 28. For these figures, we have assumed that the phases $\{\delta\chi_{n}\}$ cycle the set $\{0,\pi/2,\pi,3\pi/2\}$,

[TABLE]

As before, it is very challenging to correctly determine higher-order solutions in terms of the parameters $\{\kappa_{n}\}$ as higher order corrections include contributions from cross-terms between modes, which significantly complicate the analysis. Instead, we will focus on the more interesting phase resonance equation, where this problem is avoided, and higher-order solutions can still be found in the multi-mode case.

[TABLE]

This concludes our discussion of the frequency resonance equation. We have produced a series of solutions with increasing generality, and have provided graphical evidence of their validity. We note that some results have been intentionally left to Appendix F, since they are simple extensions of less general cases presented within this Section. We can now turn our attention to the phase resonance equation and its solutions.

IV.3 Phase resonance equation

The phase-resonance equation is derived in Section II and makes use of the action-angle variables $\{q_{\mu}\}$, instead of the fundamental frequencies $\{\omega_{\mu}\}$ to write down the Fourier modes of the self-force function. In the case of a multitude of modes in this expansion, as well as a non-vanishing initial phase offset, the equation takes the form

[TABLE]

Similarly to our approach to the frequency resonance equation, we examine the case in which we ignore all higher-order modes beyond the first one ($n=1$), and assume a zero initial phase, so eq. (19b) simplifies to

[TABLE]

The equation is again non-linear due to the additive constant, and an analytic solution cannot be provided in the general case. Instead we employ an iterative approach to establish a series of solutions in increasing order of the mode parameter $\lambda\equiv\lambda_{1}\ll 1$. At first, we ignore the term involving $\lambda$, and seek a solution to the equation:

[TABLE]

which straightforwardly yields the zeroth-order solution to the phase resonance equation (this and higher order solutions are denote by PH) in the current regime

[TABLE]

Further solutions are obtained by substituting the current solution on the right-hand side of the source equation (122), while on the left-hand side we include this solution plus the next-order correction. In the general case, we would need to solve the equation (for $n\geq 1$)

[TABLE]

where $\phi_{n}(x)$ is the complete $n^{\text{th}}$-order solution. We apply this iterative technique to find the first-order correction: we use the function $\phi_{0}(x)$ as the source for the next order solution $\phi_{1}(x)=\phi_{0}(x)+\phi_{(1)}(x)$. That means we write $\phi_{0}(x)$ as the argument on the right-hand side, while on the left-hand side we solve for $\phi_{1}(x)$, a combination of the zeroth-order $\phi_{0}(x)$ and the first-order correction $\phi_{(1)}(x)$. Throughout our current discussion, we shall use $\phi_{n}(x)$ to denote the complete solution to order $n$, while $\phi_{(n)}(x)$ would denote the $n^{\text{th}}$-order correction to $\phi_{0}(x)$, i.e. $\phi_{n}(x)=\phi_{0}(x)+\phi_{(n)}(x)$. Using this, we write:

[TABLE]

From this, we find a differential equation for $\phi_{(1)}(x)$:

[TABLE]

The details of finding this solution are described in Appendix I, here we only quote the final result:

[TABLE]

Combining this with $\phi_{0}(x)$, we find the complete first-order-in-$\lambda$ solution to eq. (122), given by (LABEL:eq:phase_solution_1). It is worth explaining that the solution has been rearranged into a form which is more useful for our discussion later on. More specifically, we have grouped together terms which oscillate around the value $\phi(x)=0$ – both (C\!\left(\mathchoice{{\hbox to4.99995pt{\hss\leavevmode\hbox{\set@color\raisebox{0.0pt}{\scalebox{2.0}{\raisebox{0.0pt}{\displaystyle\cdot}}}}\hss}}}{{\hbox to4.99995pt{\hss\leavevmode\hbox{\set@color\raisebox{0.0pt}{\scalebox{2.0}{\raisebox{0.0pt}{\textstyle\cdot}}}}\hss}}}{{\hbox to4.16661pt{\hss\leavevmode\hbox{\set@color\raisebox{0.0pt}{\scalebox{2.0}{\raisebox{0.0pt}{\scriptstyle\cdot}}}}\hss}}}{{\hbox to3.61105pt{\hss\leavevmode\hbox{\set@color\raisebox{0.0pt}{\scalebox{2.0}{\raisebox{0.0pt}{\scriptscriptstyle\cdot}}}}\hss}}}\right)\mp\nicefrac{{1}}{{2}}) and \sin(\mathchoice{{\hbox to4.99995pt{\hss\leavevmode\hbox{\set@color\raisebox{0.0pt}{\scalebox{2.0}{\raisebox{0.0pt}{\displaystyle\cdot}}}}\hss}}}{{\hbox to4.99995pt{\hss\leavevmode\hbox{\set@color\raisebox{0.0pt}{\scalebox{2.0}{\raisebox{0.0pt}{\textstyle\cdot}}}}\hss}}}{{\hbox to4.16661pt{\hss\leavevmode\hbox{\set@color\raisebox{0.0pt}{\scalebox{2.0}{\raisebox{0.0pt}{\scriptstyle\cdot}}}}\hss}}}{{\hbox to3.61105pt{\hss\leavevmode\hbox{\set@color\raisebox{0.0pt}{\scalebox{2.0}{\raisebox{0.0pt}{\scriptscriptstyle\cdot}}}}\hss}}}) are such functions (see Appendix D for details about the former). Moreover, we have collected the terms proportional to $x^{2}$ and $x$ into a single quadratic term (by completing the square, which inevitably introduces a further constant, however it is of second order in the mode parameter $\lambda$, and not relevant here). Hence the solution is comprised of a modified (shifted along the $x$-axis) quadratic term and a term $\rho_{1}$ which oscillates around $\phi=0$, and whose frequency of oscillations increases with the value of $x$.

A number of results in this section contain the symbols “$\pm$” and “$\mp$”. Similarly to the discussion of the solutions of the frequency resonance equation, these do not signify two alternative equations, but rather that the lower signs are valid for $x\leq 0$, while the upper signs hold for $x\geq 0$, both part of the same solution.

It is important to note that for large values of $x$, the oscillating terms become small, and the approximate slope of $\phi_{1}(x)$ is given by $\sqrt{\pi}\,\lambda/2$. We will see later how this slope evolves when we add higher-order-in-$\lambda$ terms, but for now let us denote $\sigma_{1}(\lambda)=\lambda/2$, which is the shift in the position of the leading parabolic term.

Applying the same iterative approach, we use the combined function $\phi_{1}(x)$ as the source for the second-order asymptotic solution $\phi_{2}(x)=\phi_{0}(x)+\phi_{(2)}(x)$, where the second-order correction $\phi_{(2)}(x)$ is comprised both of first- and second- order terms. This notation is convenient since $\phi_{0}^{\prime\prime}(x)=1$ cancels with the constant term on the right-hand side, and upon integrating twice, we obtain an expression for $\phi_{(2)}(x)$, the complete solution to second order in $\lambda$, only without the zeroth-order quadratic term. On the right-hand side, we substitute the complete first-order solution, $\phi_{1}(x)=\phi_{0}(x)+\phi_{(1)}(x)$,

[TABLE]

In fact, from eqs. (127) and (129) we can infer the general form of this differential equation (for $n>1$):

[TABLE]

The details involved in solving eq. (129) are presented in Appendix I, together with the full-form solution. Here, we omit these details, and simply quote the final combined result $\phi_{2}(x)=\phi_{0}(x)+\phi_{(2)}(x)$ in the form of a modified quadratic term, plus two terms, both of them $\mathcal{O}\!\left(\lambda^{2}\right)$ that oscillate and remain small for all $x$.

[TABLE]

We note some interesting features of this solution. Most notably, the shift in the quadratic $x$ term has changed from $\sigma_{1}(\lambda)=\lambda/2$ to $\sigma_{2}(\lambda)=\lambda/2-\sqrt{2}\,\lambda^{2}/4$ (and the slope for large $|x|$ has changed from $\sqrt{\pi}\,\sigma_{1}(\lambda)$ to $\sqrt{\pi}\,\sigma_{2}(\lambda)$). The term which is proportional to $\lambda$ has the same functional form as the corresponding term in solution (LABEL:eq:phase_solution_1), however we see that in (131) each instance of $x/\!\sqrt{\pi}$ has been modified by exactly $\pm\sigma_{1}(\lambda)$. We also have a new class of terms, all of which also oscillate around the $x$-axis, and they are proportional to $\lambda^{2}$. We will investigate the numerical performance of these solutions after we explore the next-order solution.

The solutions presented in Gair, Yunes, and Bender (2012) do not continue further than this order. Since our goal is to solve the multi-mode equation (19b), as a form of sanity check, we would like to present the solution to third order in the single-mode parameter $\lambda$, i.e. the solution to the following generating equation:

[TABLE]

The solution follows via the same technique as the previous two, and again it is not practical to present the details of the derivation here. Instead, the reader is invited to consult Appendix I for the technicalities of the calculation as well as the full solution without abbreviations. With every successive iteration we obtain new terms, and we have grouped them together in a similar fashion to before. Here we have terms which are proportional to $\lambda^{3}$, as well as modified terms, originating from the lower-order solutions. The result (133) is the complete solution to third order, for a single mode with amplitude $\lambda$, with vanishing initial phase $\delta\psi=0$.

[TABLE]

In (133), the $x$-shift of the quadratic term is now

[TABLE]

and consequently, the slope of the graph for large $x$ is given by $\sqrt{\pi}\,\sigma_{3}(\lambda)$. In fact Gair et al. (2012) numerically computed the change in the slope of the phase solution at $\lambda^{3}$ order:

[TABLE]

As a form of consistency check, we can directly compare this with the (analytically derived) factor multiplying the $\lambda^{3}$ term in the coefficient of $x$ in eq. (133) and confirm that both numerical values agree:

[TABLE]

Moreover, we can notice a parallel between the terms proportional to $\lambda$ in each of (LABEL:eq:phase_solution_1) — (133): all three brackets have the same dependence on $x$, meaning that for $\lambda=0$, they would be equal. However, in each of them, each instance of $x/\!\sqrt{\pi}$ is modified by a different function of the mode parameter $\lambda$. Compared to (LABEL:eq:phase_solution_1), each $x/\!\sqrt{\pi}$ in (131) is shifted by $\pm\,\sigma_{1}(\lambda)$, and in (133), it is shifted by $\pm\,\sigma_{2}(\lambda)$. Similar evolution can be noticed in the terms proportional to $\lambda^{2}$ in (131) and (133). In the latter, each instance of a term involving $x$ is shifted by $\sigma_{1}(\lambda)$ compared to the former. While it is hard to predict what terms at higher order would look like, because the $x$-dependence cannot be inferred from the first three solutions, we can nevertheless expect that at fourth order, the term proportional to $\lambda^{3}$ would look similar to its counterpart in (133), with each $x/\!\sqrt{\pi}$ term modified by $\pm\,\sigma_{1}(\lambda)$, and so on. It is evident that the functions $\sigma_{n}(\lambda)$ play an important role in the phase resonance solutions, and at least in the evolution of the slope of the graph for large $|x|$. In Figure 16, we plot these functions and present what they look like for different values of $\lambda$.

Let us investigate the properties of solutions (124) — (133): each of them contains a quadratic term of the form $\pi/2\,(x/\!\sqrt{\pi}\pm\sigma_{i})^{2}$. Since the magnitudes of the rest of the terms in each solution scale with $\lambda^{n}$, the functions look similar when plotted together, as can be seen in Figure 31. It is hard to distinguish them in this plot, unless we zoom in closely on a region. For example, if we focus on the region around $x=0$ we see that while $\phi_{0}(x)=x^{2}/2$ goes through the origin, the rest of the solutions deviate slightly from zero. This is due to the fact that a perturbative solution of order $n$ is only approximate, and is correct up to a remainder of $\mathcal{O}\!\left(\lambda^{n+1}\right)$. Indeed, if we consider each solution $\phi_{\{1,2,3\}}(x)$ for $x=0$:

[TABLE]

we find that they are indeed of the expected order in each case (see Figure 17 for a plot of these functions). These expressions have further higher order corrections coming from sub-leading order terms in the expansion of Fresnel cosine integrals (see eq. (175b) in Appendix D). Since these solutions are hard to examine graphically, we need to employ other techniques to compare the accuracy of the solutions. Specifically, we compare them with a purely numerical solution of the phase equation.

A simple plot like the lower one in Fig. 31 shows that our highest-order solution, eq. (133), does indeed compare well — within $\mathcal{O}\!\left(\lambda^{4}\right)$ — to a purely numerical solution (Wolfram Research, Inc., ). We need to quantify the comparison between this numerical solution and our analytical solutions as the mode parameter $\lambda$ varies.

It is important to demonstrate that each of the corrections $\phi_{(n)}(x)$ appropriately captures all of the behaviour of the solution up to order $\mathcal{O}\!\left(\lambda^{n}\right)$, i.e., that the difference between the full numerical solution $\phi_{\text{num}}(x)$ and the correction is of order $\mathcal{O}\!\left(\lambda^{n+1}\right)$ at most. To do this, we employ a technique that allows us to solve numerically for that difference (instead of naively subtracting each analytical solution from the numerical map). Consider eq. (122) which is our initial source equation:

[TABLE]

As per our recipe, we would substitute the highest-known-order solution $\phi_{n}(x)$ on the right-hand side, and solve the equation for the next-order solution $\phi_{n+1}(x)$. In reality, the solution we are looking for appears on either side of the equation, and furthermore the functions $\phi(x)$ also contain the orders beyond $n$ which we have so far been ignoring in eqs. (125) and (130):

[TABLE]

Here, $\phi_{0}(x)=x^{2}/2$ and $\phi_{(n)}(x)$ is the correction to $\phi_{0}$. $\vartheta_{n+1}(x)$ is the remainder of the complete solution which is not included into $\phi_{n}(x)$. We can easily find a generic equation satisfied by the functions $\vartheta_{n+1}(x)$ by using expression (137) in eq. (122) and making use of eq. (122). In general, $\vartheta_{n+1}(x)$ satisfies the differential equation

[TABLE]

For $n=0$ we can find directly from the source equation

[TABLE]

By solving for $\vartheta_{n+1}(x)$ and confirming that the result is $\mathcal{O}\!\left(\lambda^{n+1}\right)\forall x$ we verify that $\phi_{n}(x)$ capture the complete $\mathcal{O}\!\left(\lambda^{n}\right)$ behaviour. The plots for $\vartheta_{n}(x)/x$ (to remove the expected linear-in-$x$ general trend) are shown in Figure 18, where we have performed the computation for 4 distinct values of $\lambda$: $0.1,0.05,0.01,$ and $0.005$.

The corrections $\phi_{(n)}(x)$ themselves are plotted (for $\lambda^{(1)}=0.1$) in Figure 32, found in the Appendix. Similarly, plots of the oscillating functions $\rho_{1}(x)$, $\rho_{2}(x)$, and $\rho_{3}(x)$ can be found in Figures 33, 34, and 35, respectively, with other relevant discussions in Appendix II.

Before moving on to the case of multiple modes, it is intriguing to look into several other toy models. So far we have only looked at the resonance equation with a single mode $n=1$. We now consider the same equation, but with the sole mode $n=r$ on the right-hand side. The leading-order solution is still (124).

[TABLE]

The explicit solution of the above equation is redundant, since it can be obtained from eq. (122) by setting

[TABLE]

and therefore its corrections could be obtained from eqs. (LABEL:eq:phase_solution_1) – (133) by applying these transformations. The explicit form of the solutions is presented in Appendix III. In Figure 19 we show these solutions (at third order, i.e. complete) for different values of the mode number $r$. The modified $\sigma$-terms are presented in Figure 36, while the numerical validation of these solutions can be found in Figure 37.

Another toy model which deserves attention is the case of two distinct modes on the right-hand side:

[TABLE]

The leading-order solution is unchanged, while the higher-order corrections are presented in Appendix IV. We note that at first order, the correction is just the sum of the corrections due to each of the separate oscillating modes on the right-hand side. This, however, is not the case at higher order, where cross-terms appear. It is important to note that these cross-terms are not sub-leading to the other quadratic terms in the solution. This is in contrast to the frequency-resonance equation, where these cross-terms were sub-leading in this respect (see Appendix III). In Figure 38 we can examine the $\sigma$-terms associated with this equation, while the validity of the solutions is presented in Figure 39.

Our objective in this paper is to extend this approach beyond the scope which has been developed thus far. So far we have explored the solution to third order in the mode parameter. However, we are still limited to the case of a single mode, and we are assuming that the initial phase is $\phi(0)=0$. We now relax the first of these restrictions and allow for an infinite number of oscillating modes, each with a distinct parameter $\lambda_{i}$. This is achieved through allowing for multiple modes in the Fourier expansion, thus lifting the assumption that the contribution derives primarily from the first Fourier mode.

We derived the phase resonance equation in Section II. In the regime of an unrestricted number of oscillating modes, still assuming that the initial phase is 0, the form of the phase resonance equation is given by

[TABLE]

Similar to the previous equation under consideration, this is a second order non-linear differential equation, and therefore finding an exact solution is not possible, and a perturbative solution quickly grows to become heavily involved. For brevity, here we only present the final solutions, while details about the solutions and relevant discussion are collected in Appendix G. At zero’th order in the mode parameters $\{\lambda_{n}\}$, we ignore any oscillating terms on the right-hand side, and the resulting differential equation is identical to that in the single-mode case, eq. (123); the solution is obtained by integrating twice:

[TABLE]

Building onto this result, we construct the equivalent to eq. (125), which in the multi-mode case is now given by

[TABLE]

We will solve this consecutively for $m=0,1$ and 2, and we will obtain the respective solutions to the phase resonance equation for the multi-mode case. At first instance, we use $\bm{\phi}_{0}(x)=x^{2}/2$ as the argument of the trigonometric function on the right-hand side, and investigate the solution to this equation to first order in $\{\lambda_{n}\}$, which will yield the first-order correction $\bm{\phi}_{(1)}(x)$:

[TABLE]

The details of finding the first-order correction can be found in Appendix V. Combining $\bm{\phi}_{(1)}(x)$ with (144), we establish the complete first-order solution.

[TABLE]

Here $\bm{\rho}_{1}^{n}(x)$ is an oscillating term proportional to $|\lambda|$, equivalent to $\rho_{1}(x)$, and given explicitly in Appendix I. Compared to (144), eq. (147) has a modified quadratic term and a part which oscillates around 0 for all $x$.

We apply the same technique to find the second-order solution, namely $\bm{\phi}_{(2)}(x)$, which contains terms that are quadratic in the parameters $\{\lambda_{n}\}$. Similarly to previous iterations, we use the current highest-order solution $\bm{\phi}_{1}(x)$ to establish a differential equation for $\bm{\phi}_{(2)}(x)$. The differential equation we want to solve is given by

[TABLE]

We should remark again that the solution of this equation would yield $\bm{\phi}_{(2)}(x)$, the second-order correction to the quadratic solution. It will contain terms which are first-order in the mode parameters, as well as other terms which are of second order. This is an artefact of the method we use for asymptotically establishing higher-order corrections to our source equation. In practice, we need to solve the equation

[TABLE]

Then, by adding the zero’th order part $\bm{\phi}_{0}(x)$ to the solution of eq. (149) we will obtain the complete solution to second order in $\{\lambda_{n}\}$. Solving this equation is rather involved, therefore it is better to leave the practicalities of obtaining the solution to Appendix V. Here we only present the final solution:

[TABLE]

While the solution to second order is already quite involved, we continued our calculations and derived the solution to third order, so we can draw a comparison with the single-mode result eq. (133). The generating equation involves the new correction $\bm{\phi}_{(3)}(x)$ on the left-hand side, and the current solution, $\bm{\phi}_{2}(x)$ as the argument on the right-hand side:

[TABLE]

The solution is rather long, and hence we have placed it, together with explanations, in Appendix V. The correction at third order is given by eq. (249).

After having derived these solutions, it is important to verify their validity numerically, using the same method as in the single-mode case. The results from this verification are presented in Figure 20, where we have shown that the equation is valid for several different values of the total number of modes. To depict the effect a larger number of modes has on the solution, we have used the same initial value of the resonant flux modification $\lambda$.

So far we have solved the single- and multi-mode equations up to third order in the parameters $\{\lambda_{n}\}$. It is important to carefully consider these solutions we have obtained thus far, and evaluate the need for further iterations. Let us investigate the governing equation

[TABLE]

and consider those terms for which the small-angle approximation for the trigonometric functions breaks down:

[TABLE]

This occurs when $n\,\lambda_{1}\sim 1$. Therefore, we need to impose a restriction on the upper bound of $n$:

[TABLE]

Although the cross-term corrections here are of the same magnitude as the single-mode solutions, the bound is comparable to that in the frequency equation (cf. eq. (86)). We invoke the coefficient scaling law which we discussed in Section IV.1 to link this to values of $\lambda_{1}$:

[TABLE]

for the suggested values of $2<\alpha<3$ (Drasco and Hughes, 2006). We can therefore argue that beyond that order, we would not be able to distinguish our solutions, valid in the regime where the small-angle approximation can be applied, from effects which arise from the violation of this assumption. For this reason, higher-order corrections cannot be easily calculated without restricting the total number of modes $n$. Finally, proceeding beyond third order is pointless since for the range of values we are investigating, $\lambda^{3}\sim\epsilon^{\nicefrac{{1}}{{2}}}\sim\mathcal{O}\!\left(10^{-3}\right)$, hence in calculating further corrections we need to account for terms arising from higher-order mass-ratio corrections. In this article, we only investigate solutions up to third order inclusive, and disregard higher-order corrections.

We now return to the original differential equation (19b), and again consider a single mode, however this time we allow for a non-vanishing initial phase $\delta\psi\neq 0$

[TABLE]

The zero’th-order solution, when we ignore the term proportional to $\lambda$ on the right-hand side, is again given by $\varphi_{0}(x)=x^{2}/2$. To establish first- and higher-order corrections, we proceed according to the same procedure as before, solving equations of the form:

[TABLE]

Using this, we find the first-order solution below.

[TABLE]

We notice that the new slope $\varsigma_{1}(\lambda)$ involves the constant $\delta\psi$, and that the oscillating term $\varrho_{1}$ has a different form, now involving a Fresnel sine term too. They reduce to their zero-phase counterparts for $\delta\psi=0$. Following our established algorithm, we can easily derive the second- and third-order corrections to $\varphi_{0}(x)$. Details of these derivations are left to Appendix V, here we only present the final results.

[TABLE]

It is important to investigate how the slope of these solutions evolves with $\delta\psi$. At each order we encounter a different function of the initial phase. Figure 21 shows how these functions behave over the range of $\delta\psi$. Furthermore, we need to verify that these solutions are valid for different values of the initial phase. We computed the quantities $\vartheta_{n}(x)$ for different values of the initial phase (keeping $\lambda$ constant) using the same method as before, and show the results in Figure 22.

Finally, we approach the phase resonance equation in the most general regime: including an infinite number of modes each with a different non-zero phase:

[TABLE]

The details of the solution are too lengthy to be presented in enough detail here and can be found in Appendix VII: the first-order correction is given by eq. (255) and the second-order correction by eq. (256). Upon closer examination these bear the same features of our previously establishes solutions. Finally, we have established the third-order solution to this equation, $\bm{\varphi}_{3}(x)$. It is important to note that solution (160) is the most general solution to the phase resonance equation that we could find. It allows for multiple modes and a non-vanishing initial phase, and treats perturbations up to and including third order. The third-order slope function $\bm{\varsigma}_{3}$ is given by eq. (257) in Appendix VII. Furthermore, each of the three corrections are given by eqs. (258), (259), and (260).

[TABLE]

This is the final and most important result of this article. While it is not easily usable in theoretical calculations, it could find applications in modelling resonances in new kludge waveforms, as this is not a feature currently implemented. We have provided a solution to the highest possible order of perturbation, which is important for meaningful numerical simulations. Furthermore, we have also provided solutions to simplified versions of this problem, which could be employed in various theoretical calculations regarding emris.

V Discussion and conclusions

In this paper we have described in Sections II and III the assumptions and conditions which form the basis of two different methods for analysis of resonances in emris. The two-timescale approach, described in van de Meent (2014) presumes that the functional dependence of the orbital characteristics can be expanded in power series of the mass-ratio and time from resonance. The instantaneous frequency model, as per Gair et al. (2012), relies on the instant values of the on-resonance frequencies as a basis for asymptotic analysis, and this is a more complicated function of the (rescaled) time.

After describing the basic approach, Section IV described in detail the solution to the frequency and phase resonance equations derived in the course of describing the instantaneous frequency approach. Some of these solutions were previously published in Gair et al. (2012). The difference between the frequency and phase resonances lies in the assumption about the phasing of the self-force. It is more natural that the phase will be determined by the instantaneous orbital phase (as in the phase treatment) than the orbital phase at some fiducial $\tau=0$ (as in the frequency investigation). Fully-detailed self-force calculations would be needed to make an informed choice between the two, but it is likely that the phase approach is physically relevant for the resonance problem. The frequency resonance solutions exhibits interesting effects, in particular the “memory” effect on the slope of the graph, which would be interesting if it really occurred in any physical situation. Here we extended these results, by including higher-order terms and a generalised form of the solutions. We have derived solutions for the case of multiple oscillating modes, given by eq. (92). Working further, we have relaxed the assumption of a vanishing initial phase, and have presented the single-mode solution in eq. (112). Finally, combining these results we presented the complete frequency solution, eq. (121).

As mentioned above, the phase resonance equation is probably more relevant physically, as it describes the evolution of the phase coordinates near resonance, and its solutions can be related to the two-timescale method results. It is also interesting from a mathematical view point, since it is a second-order non-linear differential equation which can be approached using asymptotic methods. We extended the previously established sequence of results by deriving the third-order solution for the case of a single oscillating mode, eq. (133). We presented further toy models, discussed in detail in Appendices III and IV before demonstrating an order-by-order solution to the phase resonance equation with an arbitrary (possibly infinite) number of modes, eqs. (147), (150), and (249). Following this we presented the solution to the phase equation when the alternative initial assumption is relaxed: namely, the starting phase is allowed to vary. The first- and second-order solutions are given by eqs. (157) and (158). Finally, we have derived the solution in the most general case, allowing for an arbitrary (possibly infinite) number of modes, each with an individual initial phase. This solution was given in eq. (160). In the case that there is only a single mode, this calculation agrees with the result given in Gair et al. (2012), to leading order in the parameter $\lambda$, which is equivalent to the order of $G_{1}(\mathbf{J})$ in the two-timescale analysis.

Both the instantaneous frequency approach and the two-timescale analysis are valid treatments of the emri problem — though their foundational assumptions appear different, they both have the same computational complexity and arrive at equivalent conclusions. Previous work on the two-timescale method has only been presented to the lowest order in the mass ratio $\epsilon$ and the resonance flux change parameters $\{\kappa,\lambda\}$. As has been discovered, in certain scenarios higher-order corrections in the latter can be more important than those to the former, hence it is vital to investigate the form of these corrections. This is what the current paper accomplishes.

Further work on solving the phase-resonance equation would use these equations in numerical waveform models. These could be utilised in implementing a computationally cheap approximate (kludge) model of emri waveforms which include resonance effects (Barack and Cutler, 2004; Babak et al., 2007). At the moment such models lack the necessary features to predict the change in frequency or phase on resonance, which is precisely what the results in this article offer. Therefore, this should be pursued in a future publication.

Electronic resources

Many of the results in this article are available in a Mathematica notebook, which can be found at:

http://resonances.deyanmihaylov.net/

Acknowledgements

Deyan P Mihaylov is supported by the stfc. He is thankful to Dr Christopher J Moore and to Alvin Chua for the engaging conversations about general relativity, cosmology, and gravitational wave astronomy.

Appendix A Kerr space-time

Since the background we consider in this article is exclusively the Kerr space-time, it is worth providing a summary of the metric for a rotating non-charged black hole. Alternatively, a charged and spinning black hole is described by the Kerr-Newman metric, while a non-rotating charged black hole solution is given by the Reissner-Nordstöm background. The Kerr metric, in Boyer-Lindquist coordinates (Boyer and Lindquist, 1982) $(t,r,\theta,\phi)$, is given by:

[TABLE]

and $a=\mathcal{L}_{z}/M$ is the specific black hole spin parameter. This metric is singular at the points where $\theta=\{0,\pi\}$, however these are the usual coordinate singularities associated with the spherical polar system, remedied by the introduction of a second coordinate system. Further, we encounter coordinate singularities at the points where $\Delta(r)=0$, i.e., for values of the radial coordinate $r_{\pm}=M\pm\sqrt{M^{2}+a^{2}}$. To demonstrate that these are merely coordinate singularities (and not an inherent curvature singularity of the metric), it is customary to make the change to Kerr coordinates $(v,r,\theta,\chi)$, given by (Visser, 2007):

[TABLE]

This shows that the metric could be analytically continued through the surface $r_{\pm}$ to the region $0<r<r_{+}$:

[TABLE]

The surface $r=r_{+}=M+\sqrt{M^{2}-a^{2}}$ is a null hypersurface with normal 4-vector given by

[TABLE]

and the region $r\leq r_{+}$ is (part of) the black hole region. $\Omega_{\mathrm{}}$ is the angular velocity of the black hole with respect to an observer at infinity. There are several important surfaces in the Kerr space-time. Firstly, the curvature singularity $r=0$ is in fact a ring singularity, which can be explained as follows. A rotating black hole would not be spherically symmetric, as matter would bulge around the equator, the effect quantified by the spin parameter $a$. Since a point singularity cannot support angular momentum, the lowest-dimensional shape required is a ring. This is evidenced by substituting $t\equiv const.$, $r=0$, $\theta=\nicefrac{{\pi}}{{2}}$ (i.e., the singularity associated with $\Sigma(r,\theta)=0$) in eq. (161a), obtaining $\mathop{}\!\mathrm{d}s^{2}=a^{2}\mathop{}\!\mathrm{d}\phi^{2}$, i.e., a ring of radius $a$. Another important surface is the event horizon $r_{+}=M+\sqrt{M^{2}-a^{2}}$. To understand the final feature of interest, we consider the norm of the Killing vector field $(\partial/\partial t)^{a}$ in Boyer-Lindquist coordinates:

[TABLE]

and notice that it is space-like ($(\partial/\partial t)^{2}<0$) in the following region outside the event horizon $r=r_{+}$:

[TABLE]

This region is called the ergosphere, and its outside surface is named the ergosurface, which intersects the event horizon at the poles $\theta=0,\pi$. These surfaces are illustrated on the Kerr black-hole schematic on Fig. 23.

Appendix B Geodesic constants

One immediately notices that the Kerr metric does not depend on the time coordinate, $t$, or the azimuthal angle, $\phi$. Therefore, using the Euler-Lagrange equations, it can be shown that the conjugate momenta $p_{\mu}=\partial\mathcal{L}/\partial\dot{q}^{\mu}$ associated to those coordinates are constant with respect to the affine parameter of the geodesic:

[TABLE]

where in general relativity, the Lagrangian is formally given by $\mathcal{L}(q,\dot{q})=\nicefrac{{1}}{{2}}\,g_{\mu\nu}\,\dot{q}^{\mu}\dot{q}^{\nu}$. Thus we have identified 2 quantities which are constant along geodesics of the Kerr space-time. These are related to the energy of the orbit $\mathcal{E}$, and the axial component of the angular momentum $\mathcal{L}_{z}$, respectively:

[TABLE]

For equatorial geodesics ($\theta=\nicefrac{{\pi}}{{2}}$), these two constants are sufficient to uniquely determine a given geodesic. However, for general geodesics with 3 spatial degrees of freedom, we can identify a third constant of the motion, which is not associated with a linear symmetry of the background space-time.

Apart from the two Killing vectors associated with the Kerr metric, there exists a second-order Killing tensor for this background, $Q_{\mu\nu}$. This Killing tensor gives rise to another conserved quantity, which was originally discovered in 1968 as a consequence of the separable Hamilton-Jacobi equation. This quantity is called the Carter’s constant, and can be written as (Carter, 1968):

[TABLE]

Appendix C Geodesics of the Kerr metric

Fig. 24: A geodesics of the Kerr metric, and several projections, evolved for a period of $2\times 10^{5}$ s $\approx 5.56$ hours. The orbit has eccentricity $e=0.6$ and semi-latus rectum $10M$, while the central black hole has mass $M=10^{6}M_{\odot}$ and spin $0.9M^{2}$. All distances are in units of the black-hole mass $M$. The angle of inclination is $\iota=45^{\circ}$, and the plane of inclination precesses.

The equatorial geodesics of the Kerr metric can be characterised by the constants $\mathcal{E}$ and $\mathcal{L}_{z}$, combined with the Hamiltonian, which is an additional integral of motion. General 3-dimensional geodesics, however, require the introduction of a further integral of motion, Carter’s constant, discussed in Appendix B. The equations of motion for a massless particle, as per Misner, Thorne, and Wheeler (1973), are given by

[TABLE]

In the general case, these equation cannot be solved analytically to find a closed form of the shape of the geodesic. However, it is possible to integrate the equations numerically to find an approximate shape for a Kerr geodesic, and one possible solutions is shown in Fig. C. The Kerr geodesics are not closed, but rather fill the entire parameter space $r_{\text{{min}}}<r<r_{\text{{max}}}$, $\theta_{\text{{min}}}<\theta<\theta_{\text{{max}}}$, $0<\phi<2\pi$ (which has the shape of a torus) before returning to the starting $(r,\theta,\phi)$. Moreover, the plane of inclination of the orbit is not stationary, but precesses around the black hole. In the case of a resonance, when two of the frequencies become commensurate, this configuration breaks down and the geodesic becomes a (simpler) closed curve around the black hole, defined only by 2 parameters.

Appendix D Fresnel integrals

In the course of solving the phase resonance equations, as well as in the presentation of the frequency resonance solution, we encounter integrals of trigonometric functions with quadratic arguments:

[TABLE]

These are known as the Fresnel integrals, after the French engineer and physicist Augustin-Jean Fresnel, who used them in optics. Here, we have adopted the notation where the argument of the integrand is premultiplied by $\nicefrac{{\pi}}{{2}}$. Please note that this definition is equivalent to the one used in Gair et al. (2010):

[TABLE]

These integrals are not directly solvable analytically, and therefore other methods for integration must be used. One approach to finding the two integrals in (172) is by expanding the integrands in converging Taylor series:

[TABLE]

Often, a symbol like the ones appearing on the left-hand sides of eqs. (174) are used in lieu of the actual series solutions. Hence, it is useful to know the first few terms of the series for each of the integrals:

[TABLE]

The integrands and the integral functions are plotted in Fig. 25. We make note of several features: at $x=0$, $S(x)$ starts flat, and progresses as $\sim x^{3}$, while $C(x)$ is linear. For positive or negative infinity, both functions tend to $\pm\nicefrac{{1}}{{2}}$, and can be represented by an asymptotic expansion for large absolute values of $x$:

[TABLE]

Appendix E Shorthand trigonometric notation

In eqs. (44) we define the following to make the discussion regarding the frequency mode equation easier to read:

\displaystyle S_{n}(x)\equiv\sin\!\left[n\Big{(}ax+\left(1+b\right)x^{2}\Big{)}\vphantom{\Big{(}\Big{)}^{\!\!\!2\,}}\!\right]

(44a)

\displaystyle C_{n}(x)\equiv\cos\!\left[n\Big{(}ax+\left(1+b\right)x^{2}\Big{)}\vphantom{\Big{(}\Big{)}^{\!\!\!2\,}}\!\right].

(44b)

for some arbitrary constants $a$ and $b$.

Please note that these are different from the Fresnel integrals, denoted by the same letters, but without subscripts, discussed in Appendix D. These functions are plotted for several representative values of $n$ in Fig. 26. It is worth noting several identities which we make use of in the article, but do not derive explicitly. These follow directly from the product-to-sum formulae and the double angle formulae in trigonometry.

[TABLE]

We only consider positive values of the subscripts, with the following rules for negative values: $S_{-n}(x)=-S_{n}(x)$ and $C_{-n}(x)=C_{n}(x)$. It is important to note that these functions are linearly independent, i.e., they form a basis in functional space. This is most easily verified by considering the Wronskian for the sequence of functions

[TABLE]

For $\{f_{0}(x),f_{1}(x),f_{2}(x)\}$, the Wronskian is $-(a+2(1+b)x)^{3}$, while for the first 5 functions in the series, it is $72(a+2(1+b)x)^{10}$, both of which do not vanish in the general case. Since their Wronskian does not vanish, we can conclude that these functions are linearly independent over $\mathbb{R}\,\backslash\{0\}$:

[TABLE]

It should be noted that these functions are not orthogonal. First, we show how to integrate a single factor of $S_{n}(x)$:

[TABLE]

or similarly, a single factor of $C_{n}(x)$:

[TABLE]

Here, $S(x)$ and $C(x)$ are the Fresnel integrals, whose properties and limits are discussed in Appendix D. Using these results, we can easily demonstrate that these functions do not form an orthogonal basis:

[TABLE]

which does not vanish in general. Similar expressions can be obtained for integrals of $S_{n}(x)\,C_{m}(x)$ and $C_{n}(x)\,C_{m}(x)$.

Appendix F Solutions to the Frequency resonance equation

I Third-order solution

In the body of the article we derived in detail the solution to the single-mode frequency resonance equation (25) up to second order in the parameter $\kappa$. Here we extend the same solution method to derive the behaviour of $\omega(x)$ to $\mathcal{O}\!\left(\kappa^{3}\right)$. Building onto our previous ansätze, eqs. (40) and (52), we suggest the third-order function

[TABLE]

where $S_{3}(x)=\sin\!\left[3\!\left(ax+(1+b)\,x^{2}\right)\right]$ and we assume that $c_{3}\sim\mathcal{O}\!\left(\kappa^{3}\right)$, since it needs to match terms of $\mathcal{O}\!\left(\kappa^{3}\right)$ in eq. (54). We substitute this proposed solution into our governing equation (25), and following rearrangement of both sides using familiar techniques, but this time keeping terms $\mathcal{O}\!\left(x^{-1}\right)$, we find

[TABLE]

Matching terms on both sides yields a system of coupled equations, which determines the coefficients:

[TABLE]

These yield the following expressions for the undetermined coefficients of the anzats, eq. (59):

[TABLE]

and the solution to third order in $\kappa$ can be written as

[TABLE]

The value of the initial value parameter $a_{3}$ is determined in the main body of the article.

II Frequency resonance equation with $\mathbf{n\neq 1}$

In Section IV.2 we solved the single-mode frequency resonance equation to obtain the solution (50). Instead of including the first oscillatory mode $n=1$, we now solve this equation for $n=r$:

[TABLE]

We propose a version of the ansatz (40), with a modified oscillating term, reflecting the modified governing equation:

[TABLE]

Upon substitution and carrying out the derivative and the Taylor expansion (valid since $r\kappa_{r}\sim 1/r^{\alpha-1}$), we obtain:

[TABLE]

Matching the respective coefficients in the manner already familiar from above, we obtain the following expressions:

[TABLE]

We have shown that $b\sim\mathcal{O}\!\left(\kappa_{r}^{2}\right)$ and ${c_{1}}^{(r)}\sim\mathcal{O}\!\left(\kappa_{r}/r\right)$ for the $r^{\mathrm{th}}$ oscillating mode. Note that since (Drasco and Hughes, 2006) propose a value of $\alpha\gtrsim 2$, the series ${c_{1}}^{(n)}$ converge, and our use of the Taylor expansion is valid.

III Frequency resonance equation with two distinct modes

If we investigate the same equation with two distinct modes on the right-hand side, say $n=\{r,s\}$, with $r<s$:

[TABLE]

then the solution would include two leading-order terms associated with each of the two separate oscillating terms, with corrections at sub-leading order due to cross-terms between these. Specifically, we find ${c_{1}}^{(r)}\sim\kappa_{r}/2r(1+b)$ and ${c_{1}}^{(s)}\sim\kappa_{s}/2s(1+b)$, which implies that the coefficient $b$ is given by the expression

[TABLE]

and consequently

[TABLE]

where each expression also includes sub-leading terms of order $\mathcal{O}\!\left(\kappa_{r}\kappa_{s}\right)$. This proves the important result that for the frequency equation with 2 distinct modes, all mixed terms (which depend on both $\kappa_{r}$ and $\kappa_{s}$) are of sub-leading order compared to the quadratic contributions originating from the individual modes as found earlier.

IV Frequency equation with a non-vanishing phase term

Let us consider the frequency equation with a single oscillating mode but with a non-vanishing phase term $\delta\chi$:

[TABLE]

and let us expand the cosine term on the right-hand side to obtain the following differential equation:

[TABLE]

We may naively try the ansatz which previously proved suitable for the zero-phase version of the equation, given by eq. (40), however it will not work this time, as there will be unmatched terms left on the right-hand side. Instead, we suggest the following alternative ansatz, including an extra term proportional to $C_{1}(x)$

[TABLE]

Upon differentiating and substituting on both sides of eq. (93), we obtain a more involved relation than before:

[TABLE]

Making the informed assumption that $b\sim\mathcal{O}\!\left(\kappa\right)$, we find the constants $c_{1}$ and $d_{1}$ by matching terms on both sides:

[TABLE]

which is different from the previously found value of $c_{1}$. These yield the same value of the slope correction $b$ as before:

[TABLE]

If we substitute these into our proposed ansatz, we see that we recover the proposed alternative ansatz:

[TABLE]

where the sine factor in the last term is identical to $\mathcal{S}_{1}(x)$, as defined in eq. (104a).

Appendix G Solutions to the phase resonance equation

In Section IV.3, we explain our algorithm for finding iterative solutions to the phase resonance equation, and present the solutions to third order in the parameter $\lambda$. However, we omit the details of finding these solutions, since these are often too tedious and irrelevant to our discussion. Here, we present these solutions in more detail. The zeroth-order equation (123) is straightforward to solve, by integrating both sides twice.

I Phase equation with a single mode and zero phase

We source the first-order correction from the zeroth-order solution, and arrive at eq. (127). Integrating both sides between [math] and $x$, making use of the Fresnel integrals, with initial condition $\phi^{\prime}(0)=0$, we find:

[TABLE]

Integrating once again, using the method of integration by parts, we obtain an initial form of the first-order correction:

[TABLE]

While this is indeed a valid result, we prefer to re-arrange it in a form which is more useful for our tasks. If we look at the graph of this result for $\lambda=0.1$ (Fig. 32), we can notice several features. First of all, since there are no constants dependant on $\lambda$ in eq. (197), it vanishes at $x=0$. Furthermore, we immediately notice that the solution is a linear function which is modified by a decaying oscillating function. Therefore, we split the solution into two parts: a linear part, whose coefficient is a function of $\lambda$, and another part, a collection of periodic and quasi-periodic functions (trigonometric functions and Fresnel integrals) that collectively oscillates around $\phi=0$, and whose magnitude is small compared to the linear term. Since the graph is symmetric relative to the vertical axis, it is necessary to write down solutions for both cases $x>0$ and $x<0$. The Fresnel integrals oscillate around the value $\pm\nicefrac{{1}}{{2}}$ (instead of 0) for large $\pm\,x$, so we subtract this value from each of them in the rapidly oscillating terms. Finally, we note of an emerging pattern: the argument of the Fresnel integral is $x/\sqrt{\pi}$, therefore we express the linear term and the remaining factors of $x$ in terms of this scaled variable. The upper sign is for positive $x$, while the lower sign is for $x<0$.

[TABLE]

We note that this result is strictly first order in the mode parameter $\lambda$. To construct the complete solution of the phase resonance equation to first order in $\lambda$, eq. (LABEL:eq:phase_solution_1), we combine eq. (198) with the zeroth-order function (124):

[TABLE]

We have combined the quadratic and linear terms into a single quadratic term, by implicitly adding the constant term $(\pi\lambda^{2}/8)$, which is nevertheless ignored at linear order in the mode parameter, thus obtaining the form (LABEL:eq:phase_solution_1):

[TABLE]

Finally, we make note of the functional form of the oscillating term, which will be investigated in detail later together with the functions incurred at higher order (let $X=x/\!\sqrt{\pi}$ and see Appendix II for plots and discussion):

[TABLE]

The above result allows us to solve the equation at the next higher order, which is discussed below.

Finding the solution up to second order in $\lambda$ involves substituting $\phi_{1}(x)=\phi_{0}(x)+\phi_{(1)}(x)$ on the right-hand side of the source equation (122), and solving for $\phi_{(2)}(x)$. Explicitly, eq. (129) is (we use the variant (197) of $\phi_{(1)}(x)$):

[TABLE]

We write the last term in the argument of the cosine as \lambda\,\rho_{1}\big{(}x/\!\sqrt{\pi}\big{)} and invoke our knowledge of trigonometric identities to write the right-hand side as:

[TABLE]

The function \rho_{1}\big{(}x/\!\sqrt{\pi}\big{)} oscillates around 0, and further, is pre-multiplied by the coefficient $\lambda\ll 1$, hence we can assume that it is always small (compared to the quadratic and linear $x$-terms, see Figure 30 for comparison of the relative magnitude of the terms), and use this fact to substitute the trigonometric functions involving $\rho_{1}$ for Taylor series to leading order in $\lambda$. We ignore all terms of $\mathcal{O}\!\left(\lambda^{3}\right)$ and higher, since they are not included at the current order.

[TABLE]

We carry out the same procedure with the term multiplied by $\lambda^{2}$, and neglect resulting terms proportional to $\lambda^{3}$:

[TABLE]

Thus, we are left with the integral below, which can be easily performed using standard techniques, integration by parts and utilising the Fresnel integrals notation. We keep in mind to throw out terms of $\mathcal{O}\!\left(\lambda^{3}\right)$ and higher.

[TABLE]

We need to integrate this further once again, keeping track of the initial condition, $\phi(0)=0$. We obtain the following expression for the correction, where again the upper sign is valid for $x>0$ and the lower for $x<0$:

[TABLE]

We have obtained an expression which gives the correction to the quadratic term to second order in $\lambda$, which can be seen plotted in Figure 32. Again, we seek to express this in a modified form where the linear term is explicit, and all the oscillating terms are grouped by order of the parameter $\lambda$, which allows for better understanding of the solution.

[TABLE]

We note that some of these terms, if expanded, would be $\mathcal{O}\!\left(\lambda^{3}\right)$ or higher, but we leave them in for presentational purposes. We simplify the coefficient of the linear term (the slope of the derivative $\phi_{2}^{\;\prime}(x)$), by writing the Fresnel integral in power series, and retaining only its linear term (the next higher-order term is of fifth order):

[TABLE]

When we combine this with the zeroth-order function $\phi_{0}(x)=x^{2}/2$, we find the complete second-order solution, given by eq. (131). Specifically, we can complete the square with the quadratic, linear, and constant terms, to obtain

[TABLE]

In deriving this solution, we have used familiar integration techniques, and the definition of the Fresnel integrals, however it is useful to quote one of the “tricks” for the benefit of the reader who wishes to replicate these results:

[TABLE]

It is important to pay attention to the two oscillating terms in eq. (208): the first one, has the same functional form as the oscillating term in eq. (198), and can in fact be related to $\rho_{1}(X)$ in the following way:

[TABLE]

The second oscillating term has not been encountered before, and will be investigated in detail in Appendix II:

[TABLE]

As we can see, the second-order solution $\phi_{2}(x)=\phi_{0}(x)+\phi_{(2)}(x)$ is already heavily involved. Extending the results in (Gair et al., 2010), we use this function to source the third-order solution, see below for details.

Substituting the expression (133) for $\phi_{2}(x)$ into eq. (132) gives the equation we need to solve at third order:

[TABLE]

It is important to realise that $\lambda\,\rho_{1}(\bullet)$ and $\lambda^{2}\rho_{2}(\bullet)$ are small in magnitude compared to the other terms, as is obvious when these are plotted on the same graph, see Figure 30. Hence we can treat these as a perturbation to the rest of the argument. This becomes obvious when the functions are plotted on the same graph, see Figure 30. This prompts us to expand the cosine term using the sum-to-product formulae, and expand the terms containing the perturbation:

[TABLE]

We have combined both third-order terms together, since it makes the integration more straightforward. Integrating the expression (219), subject to the boundary condition $\phi^{\prime}(0)=0$ gives the derivative of the third-order correction:

[TABLE]

Integrating this again, using a combination of integration by parts, substitution, double-angle formulae and other trigonometric identities, finally applying the current boundary condition $\phi(0)=0$, we arrive at an expression for the third-order correction to the solution (124). Below is given an initial form of the solution, and as before the upper sign holds for positive $x$, while the lower sign is valid for negative $x$ (see Figure 32 for a plot of this function):

[TABLE]

As before, we re-arrange this form of the solution to obtain groups of terms which oscillate around [math] (both for negative and positive $x$) and beside them, terms which are either linear or constant with $x$, but are otherwise functions of $\lambda$.

[TABLE]

A useful trick that we use in evaluating these integrals is the following: if $f\sim\mathcal{O}\!\left(\lambda\right)$ is a polynomial function of the mode parameter, we can establish which terms are relevant at each order of the solution:

[TABLE]

The linear and constant terms can be simplified by expanding the sine function and the Fresnel integral (see eq. (175b)):

[TABLE]

Combining eq. (222) with $\phi_{(0)}(x)=x^{2}/2$, we obtain the complete solution to third order in $\lambda$, eq. (133). In particular, we have combined the quadratic, linear, and constant terms into a single quadratic term:

[TABLE]

Finally, we pay attention to the form of the oscillating terms. If we compare eq. (222) to eqs. (198) and (208), we can immediately notice how the oscillating terms proportional to $\lambda$ and $\lambda^{2}$ are related to previous results:

[TABLE]

We are in a position to establish the functional form of the third-order oscillating term, $\rho_{(3)}(X)$:

[TABLE]

II Supplementary plots and remarks

The corrections have been plotted in Figure 32 for value of $\lambda=0.1$. Interesting is their behaviour around the origin: even though we have imposed the condition $\phi(0)=0$ at each order of the solution, it is important to note that they have corrections of the next-higher order of $\lambda$: $\phi_{(1)}(x)$ is accurate up to $\mathcal{O}\!\left(\lambda^{2}\right)$, $\phi_{(2)}(x)$ to $\mathcal{O}\!\left(\lambda^{3}\right)$, and so on. Therefore, at the origin $x=0$, $\phi_{(1)}(x)$ satisfies the initial condition and goes through (0,0), since it does not contain any constant terms in $\lambda$ (cf. eq. (198)). The next-order correction, $\phi_{(2)}(x)$, however, does not go through the origin: we can see this both from eq. (208) and from Figure 32. Note that this is not only due to the constant term $\pi\lambda^{2}/8$, but also because the oscillating term proportional to $\lambda$ in eq. (208) does not go through the origin. Similarly, $|\phi_{(3)}(0)|\sim\mathcal{O}\!\left(\lambda^{4}\right)$, and the third-order correction is also valid.

As we can see, certain function emerge in the arguments of the phase functions: the coefficients of the linear $x$ term in each of the phase corrections can be written as $\pm\sqrt{\pi}\,\sigma_{\{1,2,3\}}(\lambda)$, where the functions are given by:

$\displaystyle\sigma_{1}(\lambda)=\frac{1}{2}\,\lambda,$

(228)

$\displaystyle\sigma_{2}(\lambda)=\frac{1}{2}\,\lambda-\frac{\sqrt{2}}{4}\,\lambda^{2},$

(229)

$\displaystyle\sigma_{3}(\lambda)=\frac{1}{2}\,\lambda-\frac{\sqrt{2}}{4}\,\lambda^{2}+\frac{3\sqrt{3}-1}{16}\,\lambda^{3}.$

(134)

These have been plotted (for $\lambda=0.1$) in Fig. 16, and discussed in more detail in Section IV.3.

The last issue we will discuss in this section is the functional form of the oscillating terms $\rho_{\{1,2,3\}}(X)$, given by eqs. (200), (213), and (227). These represent the functions which arise as the oscillating terms at first, second, and third orders, respectively. Two of them, $\rho_{1}(X)$ and $\rho_{3}(X)$, oscillate around the abscissa, and are depicted in Figure 33 and Figure 35. The magnitude of the oscillations decreases with $|X|$, while their frequency increases, which means that they can safely be treated as perturbations to terms which grow quadratically or linearly with $X$. The function $\rho_{2}(X)$ differs slightly from the others in that it does not oscillate around the abscissa, but rather around the curve $1/4\pi X^{2}$, as can be seen from analytically expanding the function for large values of $|X|$ (i.e., $|X|>1$):

[TABLE]

The details about $\rho_{2}(X)$ are plotted on Figure 34, and we can see that despite the dissimilarities to the other two functions, it still decreases in magnitude with increasing $|X|$, and can be treated as a perturbation to the other terms.

III Phase resonance equation with $\mathbf{n\neq 1}$

If we consider the source equation (19b) in the single-mode version again, however this time instead of keeping the mode with $n=1$, we include only the $r^{\text{th}}$ mode, then we need to solve the following differential equation:

[TABLE]

The algorithm for solution is identical to the one presented above for the case with $r=1$, hence we only quote the results at each order. The zeroth-order correction is left unchanged, given by eq. (124). The first-order correction is:

[TABLE]

In the second-order correction we begin to notice the evolution of the terms from the original solution eq. (208):

[TABLE]

We can already notice some ways in which this result differs from the second-order correction with $r=1$, for instance the change of variables from $x/\!\sqrt{\pi}$ to $\sqrt{r}\,x/\!\sqrt{\pi}$, and the factor of $1/r$ added to the multipliers of the oscillating terms. Before we explain how this arises from the governing equation, we present the third-order result:

[TABLE]

These results can be obtained from the ones involving $\lambda_{1}\equiv\lambda$ by transforming eq. (122) with a change of variables:

[TABLE]

which yields eq. (231), hence the solutions to this equation can be found from eq. (122) by applying a change of variables. It is interesting to see how the number of the mode $r$ affects the graphical behaviour of the solution. In Figure 19 we have plotted the effect of the mode number $r$ on the correction ${\phi_{3}}^{(r)}(x)$, and in Figure 36 we see the effect of $r$ on the slope of the solution. The validity of solutions (232) – (234) is confirmed graphically in Figure 37.

IV Phase resonance equation with two distinct modes

As an interlude to solving the multi-mode version of the phase resonance equation, we present the solution for two modes on the right-hand side, labelled by integers $r$ and $s$. Without loss of generality, we can let let $r<s$.

[TABLE]

The zeroth-order solution remains unchanged, $\phi_{0}^{(r,s)}=x^{2}/2$. The first-order correction can be written in the form:

[TABLE]

i.e., the two separate solutions do not mix at first order, and effectively $\phi_{(1)}^{(r,s)}(x)=\phi_{(1)}^{(r)}(x)+\phi_{(1)}^{(s)}(x)$.

At next order, however, mixing occurs, since we need to substitute the above expression (236) in both terms on the right-hand side of eq. (142). The second-order correction to the original solution can be found to be a combination of (modified) single-mode corrections plus a new mixing term at second order.

[TABLE]

Here, $\bar{l}$ denotes the element of the set $\{r,s\}$ which is currently not $l$. It is best to interpret these results graphically.

Finally, we present the third-order solution for a source equation with two distinct modes.

[TABLE]

V Multi-mode phase equation

We know by designation that the correction $\bm{\phi}_{(1)}(x)$ will be of first order in the $\{\lambda_{n}\}$. Therefore, we ignore any terms of second or higher order when solving the differential equation (146), as they belong to higher-order corrections:

[TABLE]

Substituting for $\bm{\phi}_{0}(x)$ from eq. (144), we obtain the differential equation whose solution is the first-order correction:

[TABLE]

We integrate once to find the first derivative (we have to apply the appropriate boundary condition $\bm{\phi}^{\prime}(0)=0$), keeping in mind that the integral of a sum is simply a sum of the integrals. We obtain the following result:

[TABLE]

Integrating this once again, and keeping track of the boundary condition $\bm{\phi}(0)=0$ yields the first-order correction:

[TABLE]

As we can see, in exactly the same way as solution (236), this first-order result exhibits no mixing between the solutions for different values of $n$, which means it can be expressed as a sum over many single-mode solutions of the form (232). Combining this with the zeroth-order solution, eq. (144), we can write down the complete solution of the phase equation to first order in the case of an infinite number of Fourier-like coefficients $\{\lambda_{n}\}$. As before, it is beneficial to combine the terms proportional to $x^{2}$ and $x$ into a single term, again retaining only terms linear in $\{\lambda_{n}\}$:

[TABLE]

The second-order correction to eq. (143) requires us to solve eq. (149). The argument of the trigonometric function on the right-hand side is given by the first-order correction, eq. (242) combined with the zeroth-order solution:

[TABLE]

Even though the process is more involved than before, the method of solution here is the same: consecutively integrating the right-hand side twice, while observing the boundary conditions. After the first integration, we obtain

[TABLE]

While performing this integration, we encounter some unusual integrals, hence we quote the results here:

[TABLE]

Similarly, when we integrate the expression (244), we encounter the following integral, whose solution is made possible by the fact that we are summing over both indices $n$ and $m$, and hence we can re-label them interchangeably.

[TABLE]

Using this, the second-order correction is found to be as follows (remember that we are using $\bm{\phi}^{\prime}(0)=\bm{\phi}(0)=0$).

[TABLE]

It is important to note that this is consistent with our previous results. For $n=m=1$, we retrieve the second-order single-mode solution, eq. (208) while if we choose $n,m=\{r,s\}$, we find the second-order solution for an equation with 2 distinct modes, namely the voluminous eq. (IV). Moreover, we recognise that the oscillating terms here have similar form to the ones in the single-mode solution, hence we can introduce the notation $\{\bm{\sigma}_{n}\}$ and $\{\bm{\rho}_{n}^{m}(x)\}$.

To obtain the third-order correction, we need to solve the following multi-mode second-order differential equation:

[TABLE]

Following through with two rounds of integration yields the following solution for the third-order correction:

[TABLE]

As with the second-order result, this one reduces to the appropriate previous expressions for $n=m=p=1$, or for $n=m=p=\{r,s\}$. Moreover, it is easy to recognise the generalised form of the terms proportional to $\lambda^{n}$ and how they reduce to their single-mode counterparts. A more detailed discussion about the features and validity of these solutions for different mode counts is presented in the main text of this article.

VI Phase equation with a single mode and non-zero phase

Having considered the multi-mode zero-phase equation, we revert back to the single mode case and introduce a non-zero phase $\delta\psi$:

[TABLE]

The first-order solution is $\varphi_{0}=x^{2}/2$, and we substitute it into eq. (156) to obtain the differential equation for $\varphi_{(1)}$:

[TABLE]

The method of solution is identical to expand the trigonometric function on the right-hand side:

[TABLE]

We now proceed to integrate this exactly as before, only this time there are twice as many terms that need to be integrated. The final solution can be presented in a fashion similar to before, however this time we have $\delta\psi$-dependence both in the slope and in the oscillating terms.

[TABLE]

This solution can be used again in the source equation to find the second-order correction:

[TABLE]

We can apply the same technique to find the third-order solution. The third-order correction is given below.

[TABLE]

Details about these solutions are discussed in the main text of the article. We can now address the final problem, namely that of a phase equation with many oscillating modes, each with an independent initial phase.

VII Multi-mode phase equation with non-zero phase

The solutions to the phase resonance equation with an infinite number of modes and non-zero phases

[TABLE]

can be found using the established algorithm. Note this is the most complete solution to the phase resonance equation. The first-order correction to the zeroth-order solutions is given by the expression below.

[TABLE]

Similarly, we can find the second-order correction, given here.

[TABLE]

Finally, we get to the point where we can derive the most complete solution to the problem we are considering in this article. The third-order solution in this case is given by eq. (160).

[TABLE]

In this solution, $\bm{\varsigma}_{3}$ is the third-order slope function, dependent on the resonance flux parameters $\{\lambda_{n}\}$ and the mode numbers $\{n\}$ (this solutions does not require the mode numbers to be consecutive). The lower-order functions $\bm{\varsigma}_{1}$ and $\bm{\varsigma}_{2}$ are given as part of the corrections $\bm{\varphi}_{(1)}$ and $\bm{\varphi}_{(2)}$, respectively.

[TABLE]

The three oscillating functions $\{\bm{\varrho}_{i}\}$ are given below.

[TABLE]

[TABLE]

and

[TABLE]