1 Introduction↩︎

It has now been a decade since the first direct detection of a GW signal by the LIGO [1] that marks the beginning of GW astronomy. Over the past ten years, we are able to understand much more about astrophysics, such as the population properties of stellar-mass compact binary system [2], fundamental physics, such as BH mechanics [3], and more, by analyzing the gravitational waveforms observed by ground-based detectors such as LIGO [4], Virgo [5] and KAGRA [6].

The next big leap in the field of GW astronomy would be the commission of space-based GW detectors such as the LISA [7]. These space-based detectors target a much lower frequency band—in the millihertz range—compared to those ground-based ones. As a result, they are sensitive to different astrophysical sources of GWs. One such sources is the EMRIs [8], which are the gravitational radiation emitted by massive BHs when perturbed by smaller bodies, such as stars and BHs that are much lighter, orbiting around them.

In contrast to the gravitational waveforms coming from the merger of a stellar-mass compact binary system that we can observe at ground-based detectors effectively infinitely far away, the GW signals coming from EMRIs can remain detectable for months to years instead of just mere seconds [9], [10]. Therefore, it is crucial for us to be able to compute these EMRI waveforms accurately, such that we can compare these theoretical predictions with observations and extract properties about the sources of those EMRIs.

1.1 Primer in solving for gravitational waveforms from Kerr black holes using Green’s functions↩︎

The gravitational waveforms that we observe at spatial infinity contain two modes, namely the plus polarization \(h_{+}\) and the cross polarization \(h_{\times}\), respectively. They are encoded in the perturbed Weyl scalar \(\psi_4\) as [11] \[\frac{1}{2}\dfrac{\partial^2}{\partial t^2} \left( h_{+} - ih_{\times} \right) = \psi_4(r \to \infty).\] In his seminal work, Teukolsky showed that the equation governing the linear perturbation to the scalar \(\psi_{4}\), which is a PDE, can be solved using separation of variables and a FT [12]. Schematically, in the Kinnersley tetrad, this decomposition can be written as \[\rho^{-4}\psi_4(t, r, \theta, \varphi) = \sum_{\ell m \omega} R_{\ell m \omega}(r) {}_{-2}S_{\ell m \omega}(\theta, \varphi) e^{- i\omega t},\] where \(\rho=-(r-ia\cos\theta)^{-1}\), and \((t,r,\theta,\varphi)\) are the BL coordinates. Throughout this paper, we use geometric units where \(c = G = M = 1\).

For the angular sector \((\theta, \varphi)\), the solutions are known as the SWSHs \({}_{-2}S_{\ell m \omega}(\theta, \varphi) = {}_{-2}S_{\ell m \omega}(\theta)e^{i m \varphi}\). We refer readers to the Appendix A of Ref. [13] for more details.¹ As for the radial sector, the solutions \(R_{\ell m \omega}(r)\) are governed by an ODE, aptly referred to as the radial Teukolsky equation in literature, given by \[\label{Eq46TeukolskyRadialEquation} \left[\Delta^2\frac{d}{dr}\left(\frac{1}{\Delta}\frac{d}{dr}\right)-V_{\rm T}(r)\right]R_{\ell m\omega}(r)= - \mathcal{T}_{\ell m\omega}(r),\tag{1}\] with a potential \(V_{\rm T}\) given by \[\label{Eq46TeukolskyEffectivePotential} V_{\rm T}(r)= -\dfrac{K^2 + 4i(r-1)K}{\Delta} + 8i\omega r + \lambda_{\ell m \omega},\tag{2}\] where \(\Delta = (r - r_+)(r - r_-)\), \(r_{\pm} = 1 \pm \sqrt{1 - a^2}\), \(K = (r^2 + a^2)\omega - ma\), and \(\lambda_{\ell m\omega}\) is the separation constant from the angular sector. For the sake of simplicity, we will drop the \(\ell m \omega\) subscript when there is no risk of confusion hereinafter.

Conceptually, the inhomogeneous radial Teukolsky equation in Eq. 1 can be solved using the Green’s function method. We start with two linearly independent homogeneous solutions that satisfy one of the two boundary conditions for the inhomogeneous solution \(R^{\rm inhomo}\) that we want, respectively. In this case, we impose the boundary conditions that the solution is purely ingoing at the horizon and purely outgoing at spatial infinity, and the corresponding homogeneous solutions are denoted by \(R^{\rm in}(r)\) and \(R^{\rm up}(r)\), respectively. Specifically, these solutions have the following asymptotic forms \[\begin{align} \tag{3} R^{\mathrm{in}}(r) & = \begin{cases} B^{\mathrm{trans}}_{\mathrm{T}} \Delta^{2} e^{-i \kappa r_*}, & r \to r_+ \\ B^{\mathrm{inc}}_{\mathrm{T}} \dfrac{e^{-i\omega r_*}}{r} + B^{\mathrm{ref}}_{\mathrm{T}} r^3 e^{i\omega r_*}, & r \to \infty \\ \end{cases}, \\ \tag{4} R^{\mathrm{up}}(r) & = \begin{cases} C^{\mathrm{ref}}_{\mathrm{T}} \Delta^{2} e^{-i \kappa r_*} + C^{\mathrm{inc}}_{\mathrm{T}} e^{i \kappa r_*}, & r \to r_+ \\ C^{\mathrm{trans}}_{\mathrm{T}} r^3 e^{i\omega r_*}, & r \to \infty \end{cases}, \end{align}\] where \(B^{\rm trans,inc, ref}_{\rm T}\) and \(C^{\rm trans,inc, ref}_{\rm T}\) are the transmission, incidence, and reflection coefficients for the \(R^{\rm in}\) and \(R^{\rm up}\) solutions, respectively, and \(\kappa = \omega-ma/(2r_{+})\).

With these two solutions, we can construct a Green’s function \(G_{\rm T}(r, \tilde{r})\) as \[G_{\rm T}(r, \tilde{r}) = \begin{cases} \dfrac{1}{W_R} R^{\rm up}(r)R^{\rm in}(\tilde{r}), & r > \tilde{r} \\ \dfrac{1}{W_R} R^{\rm in}(r)R^{\rm up}(\tilde{r}), & r < \tilde{r} \end{cases},\] where \(W_R\) is the scaled Wronskian for the two Teukolsky solutions given by \[W_R = \dfrac{1}{\Delta} \left( R^{\rm in} \dfrac{dR^{\rm up}}{dr} - R^{\rm up} \dfrac{dR^{\rm in}}{dr} \right).\] The inhomogeneous solution \(R^{\rm inhomo}(r)\) is then given by \[\label{eq:inhomo95R95conv95integral95full} \begin{align} R^{\rm inhomo}(r) = & \frac{R^{\rm up}(r)}{W_R} \int_{r_+}^{r} d\tilde{r}\, \dfrac{R^{\rm in}(\tilde{r}) \mathcal{T}(\tilde{r})}{\Delta^{2}(\tilde{r})} \\ & + \frac{R^{\rm in}(r)}{W_R} \int_{r}^{\infty} d\tilde{r} \, \dfrac{R^{\rm up}(\tilde{r}) \mathcal{T}(\tilde{r})}{\Delta^{2}(\tilde{r})} \end{align}.\tag{5}\] As \(r \to \infty\), the solution becomes \[\label{eq:inhomo95R95conv95integral} R^{\rm inhomo}_{\ell m \omega}(r \to \infty) = \underbrace{\frac{1}{2 i \omega B^{\mathrm{inc}}_{\mathrm{T}}} \int_{r_+}^{\infty} d\tilde{r}\, \dfrac{R^{\rm in}(\tilde{r}) \mathcal{T}(\tilde{r})}{\Delta^{2}(\tilde{r})}}_{ Z^{\infty}_{\ell m \omega} } r^3 e^{i\omega r_*},\tag{6}\] where we have substituted an expression for \(W_{R}\), and we can see that \(R^{\rm inhomo}_{\ell m \omega}\) indeed satisfies the purely outgoing boundary condition at spatial infinity.

Computationally, one will need to first solve the homogeneous radial Teukolsky equation to get \(R^{\rm in, up}(r)\) and the Wronskian \(W_{R}\) before performing the convolution integral over the source term. This can be done using various methods, such as the MST method [14]–[16], the SN formalism [13], [17], [18], and recently a method based on analytical series expansion [19]. When the source term \(\mathcal{T}\) is compact (i.e., nonvanishing only at finite values of \(r\)), for instance, the bound motion of a test particle orbiting around a BH for EMRI waveform modeling, Eq. 5 is perfectly fine for these kind of calculations (see, for example, Ref. [20]).

However, Eq. 5 is no longer suitable for numerical computations when \(\mathcal{T}\) is extended and does not decay fast enough. In these cases, the integral for \(Z^{\infty}_{\ell m \omega}\) is divergent.² A prototypical example of this scenario would be the radial infall of a test particle from infinity towards a BH. In fact, this was the very motivation that led to the development of the SN formalism [17], [18], which gives a well-behaved source term and a convergent convolution integral even when the source is extended, on top of providing an efficient numerical scheme to solve the homogeneous Teukolsky equation in Eq. 1 .

1.2 This work↩︎

In this paper, we revamp the SN formalism for the driven case to take full advantages of the formalism for computing gravitational radiation from Kerr BHs. Specifically, we give a derivation of the source term for the inhomogeneous SN equation that is valid for any equation of motion (i.e., not necessary a geodesic), and introduce a new scheme for solving gravitational waveforms that bypasses the additional integration that was required to obtain the appropriate source term, which is a common criticism of the formalism.

This paper is organized as follows. In Sec. 2.1, we first review the basics of the SN formalism. Then in Sec. 2.2, we describe our new scheme for solving the inhomogeneous SN equation using integration by parts, followed by the recipes for computing gravitational waveforms and fluxes at infinity using the SN formalism in Sec. 2.3. We present in Sec. 3.1 and Sec. 3.2 our results for bound and unbound orbits, respectively. Finally, we discuss some applications, limitations, and extensions of this work in Sec. 4.

2 Sasaki-Nakamura formalism↩︎

Here we review the basics of the SN formalism for the sake of completeness. In essence, the formalism introduces a new variable \(X_{\ell m \omega}\) in place of \(R_{\ell m \omega}\) used in the Teukolsky formalism. This new variable is constructed such that the ODE it satisfies has a short-ranged potential and a source term that gives a convergent integral when using the Green’s function method. We refer readers to Refs. [17], [18] for the detailed construction of the variable.

2.1 Basics of the Sasaki-Nakamura formalism↩︎

The variable \(X_{\ell m \omega}\) satisfies the SN equation, which is given by \[\label{Eq46Inhomogeneous95SN} \left[\frac{d^2}{d{r_*}^2}- \mathcal{F}_{\ell m \omega}\frac{d}{dr_*}- \mathcal{U}_{\ell m \omega}\right]X_{\ell m\omega}(r_*)=\mathcal{S}_{\ell m\omega}(r),\tag{7}\] where \[\begin{align} \mathcal{F}(r)&= \frac{\eta'}{\eta}\frac{\Delta}{r^2+a^2},\\ \mathcal{U}(r)&=\frac{\Delta U_1}{(r^2+a^2)^2}+G^2+\frac{\Delta G'}{r^2+a^2}-\mathcal{F}G,\\ G(r)&=-\frac{2(r-1)}{r^2+a^2}+\frac{r\Delta}{(r^2+a^2)^2},\\ U_1(r)&=V_{\rm T} +\frac{\Delta^2}{\beta}\left[\left(2\alpha+\frac{\beta'}{\Delta}\right)'-\frac{\eta'}{\eta}\left(\alpha+\frac{\beta'}{\Delta}\right)\right], \end{align}\] with \[\begin{align} \alpha &=3i K'+\lambda+\frac{6\Delta}{r^2}-i\frac{K\beta}{\Delta^2},\\ \beta &=\Delta\left(-2i K+\Delta'-\frac{4\Delta}{r}\right), \\ \eta &=c_0+\frac{c_1}{r}+\frac{c_2}{r^2}+\frac{c_3}{r^3}+\frac{c_4}{r^4}, \end{align}\] and \[\begin{align} c_0 & = -12i\omega + \lambda(2 + \lambda) -12 a\omega \left( a\omega - m \right), \\ c_1 & = 8iam\lambda + 8ia^2 \omega( 3 - \lambda), \\ c_2 & = -24 i a \left( a \omega - m\right) + 12 a^2 \left[ 1 - 2 \left( a \omega - m \right)^2 \right], \\ c_3 & = 24 i a^3 \left( a \omega - m \right) - 24a^2, \\ c_4 & = 12a^4. \end{align}\] Moreover, the ODE is written with respect to the tortoise coordinate \(r_{*}\) given by \[\label{Eq46Tortoise} \begin{align} r_{*}(r)& = \int^{r} \frac{\tilde{r}^2+a^2}{\Delta}d\tilde{r},\\ & = r+\frac{2r_{+}}{r_{+}-r_{-}}\ln\frac{r-r_{+}}{2}-\frac{2r_{-}}{r_{+}-r_{-}}\ln\frac{r-r_{-}}{2}. \end{align}\tag{8}\] The Teukolsky variable \(R_{\ell m\omega}(r)\) in Eq. 1 can be constructed from the SN variable \(X_{\ell m\omega}\left(r_{*}(r)\right)\) in Eq. 7 using \[\label{Eq46SNtransformationRtoX} R_{\ell m\omega}(r)=\Lambda^{-1}\left[X_{\ell m\omega}(r_{*}(r))\right]+\frac{\left(r^2+a^2\right)^{3/2}}{\eta}\mathcal{S}_{\ell m\omega},\tag{9}\] where \(\Lambda^{-1}\) is the differential operator for the inverse SN transformation defined as \[\Lambda^{-1}\left[X_{\ell m\omega}\right]=\frac{1}{\eta}\left[\frac{\alpha\Delta+\beta'}{\sqrt{r^2+a^2}}X_{\ell m\omega}-\frac{\beta}{\Delta}\left(\frac{\Delta X_{\ell m\omega}}{\sqrt{r^2+a^2}}\right)'\right].\] Comparing Eq. 9 with the homogeneous case, e.g., in Ref. [13], we see that there is now an additional contribution coming from the source term \(\mathcal{S}_{\ell m \omega}\). By construction, \(\mathcal{S}_{\ell m \omega}\) decays fast enough as \(r \to \infty\) such that it does not contribute to Eq. 9 when evaluated at spatial infinity. Importantly, this implies that \[\label{eq:inhomo95X95to95inhomo95R95asym95relation} R_{\ell m \omega}(r\to\infty) = \lim_{r \to \infty} \Lambda^{-1}\left[X_{\ell m\omega}(r_{*}(r))\right].\tag{10}\]

The derivation of the expression relating the SN source term \(\mathcal{S}_{\ell m \omega}\) with the Teukolsky source term \(\mathcal{T}_{\ell m \omega}\) can be found in Appendix 6. Here, we just state the result, which is \[\label{Eq46Source95S} \mathcal{S}_{\ell m\omega}=\frac{\eta\Delta\mathcal{W}}{(r^2+a^2)^{3/2}r^2}\exp\left(-i\int^r\frac{K}{\Delta}d\tilde{r}\right),\tag{11}\] where we define an auxiliary function \(\mathcal{W}(r)\) related to the Teukolsky source term \(\mathcal{T}_{\ell m\omega}\) by \[\label{Eq46d2W} \frac{d^2 \mathcal{W}}{dr^2}=-\frac{r^2}{\Delta^2}\mathcal{T}_{\ell m\omega}(r)\exp\left(i\int^r\frac{K}{\Delta}d\tilde{r}\right).\tag{12}\] One in principle needs to integrate the ODE in Eq. 12 to obtain the SN source term for Eq. 7 . The Teukolsky source term \(\mathcal{T}_{\ell m \omega}\) itself for a point particle is given by [22] \[\label{Eq46T95definition} \begin{align} \mathcal{T}_{\ell m\omega}(r)=\;&\mu\int_\gamma d\tau\;e^{i\omega t(\tau)-i m\varphi(\tau)} \\ &\Delta^2 \left\{\left(A_{nn0}+A_{n\bar{m}0}+A_{\bar{m}\bar{m}0}\right)\delta(r-r(\tau))\right.\\ &+\left[\left(A_{n\bar{m}1}+A_{\bar{m}\bar{m}1}\right)\delta(r-r(\tau))\right]'\\ &\left.+\left[A_{\bar{m}\bar{m}2}\delta(r-r(\tau))\right]''\right\}, \end{align}\tag{13}\] where \(\mu\) is the mass of the particle and \(\gamma\) denotes its trajectory. The expressions of \(A_{nn0}\), \(A_{n\bar{m}0}\), \(A_{\bar{m}\bar{m}0}\), \(A_{n\bar{m}1}\), \(A_{\bar{m}\bar{m}1}\), \(A_{\bar{m}\bar{m}1}\) can be found in Appendix 7.

We can solve the inhomogeneous SN equation using the Green’s function method. Similarly, we need \(X^{\rm in}\) and \(X^{\rm up}\) that satisfy the purely ingoing boundary condition at the horizon and purely outgoing boundary condition at spatial infinity, respectively. Asymptotically, they are given by \[\label{Eq46Xin95Asymptotic} X^{\mathrm{in}}(r_{*})=\begin{cases} B^{\mathrm{trans}}_{\mathrm{SN}}e^{-i \kappa r_{*}}& r_{*}\to-\infty\\ B^{\mathrm{inc}}_{\mathrm{SN}}e^{-i\omega r_{*}}+B^{\mathrm{ref}}_{\mathrm{SN}}e^{i\omega r_{*}}& r_{*}\to\infty \end{cases},\tag{14}\] and \[\label{Eq46Xup95Asymptotic} X^{\mathrm{up}}(r_*)=\begin{cases} C^{\mathrm{ref}}_{\mathrm{SN}}e^{-i \kappa r_*}+C^{\mathrm{inc}}_{\mathrm{SN}}e^{i \kappa r_{*}}& r_{*}\to -\infty\\ C^{\mathrm{trans}}_{\mathrm{SN}}e^{i \omega r_{*}} & r_{*}\to\infty \end{cases},\tag{15}\] where \(B^{\rm trans,inc, ref}_{\rm SN}\) and \(C^{\rm trans,inc, ref}_{\rm SN}\) are the transmission, incidence and reflection coefficients for the \(X^{\rm in}\) and \(X^{\rm up}\) solutions, respectively. The inhomogeneous solution \(X^{\rm inhomo}(r_{*})\) is then given by \[\label{eq:full95inhomo95X95using95green95func} \begin{align} X^{\rm inhomo}_{\ell m\omega}(r_{*})=&\frac{X_{\ell m\omega}^{\mathrm{up}}(r_{*})}{W_X}\int_{-\infty}^{r_{*}} X_{\ell m\omega}^{\mathrm{in}}(\tilde{r}_*)\frac{\mathcal{S}_{\ell m\omega}(\tilde{r}_*) }{\eta}d\tilde{r}_* \\ &+\frac{X_{\ell m\omega}^{\mathrm{in}}(r_{*})}{W_X}\int_{r_{*}}^\infty X_{\ell m\omega}^{\mathrm{up}}(\tilde{r}_*)\frac{\mathcal{S}_{\ell m\omega}(\tilde{r}_*)}{\eta}d\tilde{r}_*, \end{align}\tag{16}\] where \(W_X\) is the scaled Wronskian³ defined by \[W_X = \dfrac{1}{\eta} \left[ X_{\ell m\omega}^{\mathrm{in}} \frac{dX_{\ell m\omega}^{\mathrm{up}}}{dr_{*}} - X_{\ell m\omega}^{\mathrm{up}} \frac{dX_{\ell m\omega}^{\mathrm{in}}}{dr_{*}} \right] =\frac{2i\omega}{c_0}B^{\mathrm{inc}}_{\mathrm{SN}}C^{\mathrm{trans}}_{\mathrm{SN}} .\]

In particular, when \(r_{*}\to \infty\), the inhomogeneous SN solution becomes \[\begin{gather} \label{Eq46X94Infty} X^{\rm inhomo}_{\ell m\omega}(r_{*}\to \infty) = \\ \underbrace{\frac{c_0}{2i\omega B_{\mathrm{SN}}^{\mathrm{inc}}}\int_{-\infty}^\infty\frac{X_{\ell m\omega}^{\mathrm{in}}(r_{*})\mathcal{S}_{\ell m\omega}(r_{*})}{\eta}dr_{*}}_{X_{\ell m\omega}^{\infty}} e^{i\omega r_{*}}. \end{gather}\tag{17}\] Using Eq. 10 , one can relate the asymptotic amplitude at infinity \(X^{\infty}_{\ell m \omega}\) with \(Z^{\infty}_{\ell m \omega}\), which means \[R^{\rm inhomo}_{\ell m\omega}(r \to \infty) = -\frac{4\omega^2}{c_0}X_{\ell m\omega}^{\infty} r^3 e^{i \omega r_{*}}.\] The GW polarizations \(h_{+}\) and \(h_{\times}\) can then be expressed as \[\begin{gather} \label{Eq46h} h_{+} - ih_{\times} = \\ -\frac{2}{r}\sum_{\ell m}\int_{-\infty}^\infty \frac{Z_{\ell m\omega}^{\infty}}{\omega^2}{_{-2}}S_{\ell m \omega}(\theta)e^{-i\omega(t-r_{*}) + im\varphi}d\omega, \end{gather}\tag{18}\] where \[\label{Eq46Z} Z_{\ell m\omega}^{\infty}=-\frac{4\omega^2}{c_0}X_{\ell m\omega}^{\infty}.\tag{19}\]

2.2 New scheme for solving the inhomogeneous Sasaki-Nakamura equation using integration by parts↩︎

Conventionally, solving for gravitational waveforms \(h_{+,\times}\) using the SN formalism requires first integrating Eq. 12 for \(\mathcal{W}\) [and hence \(\mathcal{S}\) through Eq. 11 ], often numerically, and then integrating the convolution integral of some homogeneous solution \(X\) with the source term \(\mathcal{S}\) in Eq. 17 for \(X^{\infty}_{\ell m \omega}\). Comparing with Eq. 6 , where the source term \(\mathcal{T}\) in the convolution integral often has an analytical expression, the SN formalism seems to be at a disadvantage. However, this does not have to be the case. Here, we show that by using IBP twice with the help of an auxiliary function, one can convert the convolution integral in the SN formalism to use the Teukolsky source term \(\mathcal{T}\).

We first define two auxiliary functions \(Y_{\ell m\omega}^{\rm in/up}(r)\), respectively, where \[\label{Eq46ODEforY} Y_{\ell m\omega}^{\rm in/up\;\prime\prime}(r)\equiv \frac{X_{\ell m\omega}^{\rm in/up}(r)}{r^2\sqrt{r^2+a^2}}\exp\left(-i\int^r\frac{K}{\Delta}dr\right).\tag{20}\] Furthermore, we replace \(X^{\rm in/up}_{\ell m \omega}\) with \(Y^{\rm in/up}_{\ell m \omega}\) in Eq. 16 , we can obtain \[\label{eq:full95inhomo95X95with95Y} \begin{align} X^{\rm inhomo}_{\ell m\omega}(r_{*})=&\frac{X_{\ell m\omega}^{\mathrm{up}}(r_{*})}{W_X}\int_{r_{+}}^{r(r_{*})} Y^{\rm in\;\prime\prime}_{\ell m \omega} \mathcal{W}(r) dr \\ &+\frac{X_{\ell m\omega}^{\mathrm{in}}(r_{*})}{W_X} \int_{r(r_{*})}^{\infty} Y^{\rm up\;\prime\prime}_{\ell m \omega} \mathcal{W}(r) dr. \end{align} \tag{21}\] Now Eq. 21 is written in a suggestive form. We can apply IBP twice to swap the differentiation (with respect to \(r\)) from \(Y\) to \(\mathcal{W}\), at the expense of picking up extra boundary terms where \[\begin{gather} \label{eq:Xinhomo95IBP} X^{\rm inhomo}_{\ell m\omega}(r_{*}) = \frac{X_{\ell m\omega}^{\mathrm{up}}(r_{*})}{W_X} \int_{r_{+}}^{r(r_{*})} Y^{\rm in}_{\ell m \omega} \frac{d^2\mathcal{W}(r)}{dr^2} dr \\ + \frac{X_{\ell m\omega}^{\mathrm{up}}(r_{*})}{W_X} \left[ Y^{\rm in \;\prime}_{\ell m\omega}(r)\mathcal{W}(r)-Y^{\rm in}_{\ell m\omega}(r)\mathcal{W}'(r) \right]_{r_{+}}^{r(r_{*})} \\ + \frac{X_{\ell m\omega}^{\mathrm{in}}(r_{*})}{W_X} \left[ Y^{\rm up \;\prime}_{\ell m\omega}(r)\mathcal{W}(r)-Y^{\rm up}_{\ell m\omega}(r)\mathcal{W}'(r) \right]^{\infty}_{r(r_{*})} \\ +\frac{X_{\ell m\omega}^{\mathrm{in}}(r_{*})}{W_X} \int_{r(r_{*})}^{\infty} Y^{\rm up}_{\ell m \omega} \frac{d^2\mathcal{W}(r)}{dr^2} dr. \end{gather}\tag{22}\] Specifically, we are interested in the case when \(r_{*}\to \infty\), i.e., \[\label{Eq46Xinf32after32IBP} \begin{align} X^{\infty}_{\ell m\omega} & = \frac{c_0}{2i\omega B^{\mathrm{inc}}_{\mathrm{SN}}}\left[Y^{\rm in \;\prime}_{\ell m\omega}(r)\mathcal{W}(r)-Y^{\rm in}_{\ell m\omega}(r)\mathcal{W}^{\prime}(r)\right]_{r_{+}}^\infty \\ & \; +\frac{c_0}{2i\omega B^{\mathrm{inc}}_{\mathrm{SN}}}\int_{r_{+}}^\infty Y^{\rm in}_{\ell m\omega}(r)\frac{d^2\mathcal{W}(r)}{dr^2}dr. \end{align}\tag{23}\] This is the key result of the paper—if we can discard the boundary terms (which later in the text we show that this is justified in some cases), then we can calculate \(X^{\infty}_{\ell m\omega}\) without having to solve for \(\mathcal{W}\). Additionally, the new auxiliary function \(Y^{\rm in}_{\ell m\omega}\) introduced here does not depend on the source term and can be constructed easily from the homogeneous solution \(X^{\rm in}_{\ell m\omega}\), which will be the subject of Sec. 2.2.1.

By inserting Eq. 12 into Eq. 23 , the convolution integral in our new scheme can be written as

where we define for convenience⁴ \[\label{eq:mathcalY} \mathcal{Y}(r)=r^2Y^{\rm in}(r)\exp\left(i \int^r\frac{K}{\Delta}d\tilde{r}\right).\tag{25}\] Notice that we have exchanged the order of the \(d\tau\) integral and the \(dr\) integral in the last equality of Eq. 24 . This allows us to eliminate the Dirac delta function and its derivative and evaluate the integral along the particle trajectory.

In fact, by simplifying the expression enclosed in the square brackets in the last equality of Eq. 24 , we can obtain a very elegant expression for \(I\) as

where \[\label{Eq46N95and95M} \begin{align} \mathcal{N} & =u^t-a\sin^2\theta u^\varphi+\frac{\Sigma}{\Delta}u^r,\\ \bar{\mathcal{M}} & =i a\sin\theta u^t-i\left(r^2+a^2\right)\sin\theta u^\varphi+ \Sigma u^\theta, \end{align}\tag{27}\] with \(u\) denoting the four velocity of the particle and \(\Sigma = r^2 + a^2 \cos^2 \theta\). Note that Eq. 26 holds also for nongeodesic motions. The \(W\) terms (not to be confused with \(\mathcal{W}\)) represent the coupling between the auxiliary function \(Y(r)\) and certain structures of the source, while \(\mathcal{N}\) and \(\bar{\mathcal{M}}\) contain the information about the motion of the particle along the \(n_{\mu}\) and \(\bar{m}_{\mu}\) direction in the NP tetrad, respectively. The expressions for the \(W\) terms can be found in Appendix 7.

2.2.1 \(Y(r)\) function↩︎

A core ingredient of our new SN-IBP approach is the auxiliary functions \(Y^{\rm in, up}\), which are constructed as the solutions to the second order ODE in Eq. 20 subjecting to different initial conditions, respectively. Here, we give a prescription on how to solve for these functions.

For the \(Y^{\rm in}\) function, we impose the initial conditions that \(Y^{\rm in}(r \to \infty) = {Y^{\rm in}}'(r \to \infty) = 0\).⁵ Unfortunately, Eq. 20 needs to be integrated numerically. To speed up the computation, we expand \({Y^{\rm in}}''\) near infinity as \[\begin{gather} \label{Eq46Yin3939AsymptoticExpansion} {Y^{\rm in}}''(r\to\infty) = \\ \frac{B_{\rm SN}^{\rm ref}}{r^3}\sum_{j=0}^\infty\frac{Y^\infty_{+,j}}{r^j}+\frac{B_{\rm SN}^{\rm inc}e^{4i\omega \ln 2-2i\omega r}}{r^{3+4i\omega}}\sum_{j=0}^\infty\frac{Y^\infty_{-,j}}{r^j}, \end{gather}\tag{28}\] where the coefficients \(Y^{\infty}_{\pm, j}\) are given in Appendix 8.1. This allows us to start the numerical integration for \(Y^{\rm in}\) at a smaller outer boundary since we can analytically integrate Eq. 28 to evaluate the proper initial values to use at the outer boundary.

In addition, we found that it is easier to integrate the ODE in \(r_{*}\) instead, which is now given by \[\begin{gather} \label{Eq46Y95rs95ODE} \frac{d^2Y}{dr_{*}^2}= \frac{2(r^2-a^2)}{(r^2+a^2)^2}\frac{dY}{dr_{*}}\\ +\frac{\Delta^2X(r_{*})}{r^2(r^2+a^2)^{5/2}}\exp\left(-i\int^r\frac{K}{\Delta}d\tilde{r}\right), \end{gather}\tag{29}\] where we have omitted the \(\left\{\rm{in, up}\right\}\) superscript for simplicity.⁶ As an example, Fig. 1 shows the solution for \(Y^{\rm in}(r_{*})\) and \({Y^{\rm in}}'(r_{*})\) with \(\ell=m=2\), \(a/M=0.9\), and \(M\omega=1\), \(0.5\), and \(0.1\). We see that both \(Y(r_{*}\to\infty)\) and \(Y'(r_{*}\to\infty)\) converge to zero, while \(Y(r_{*}\to-\infty)\) and \(Y'(r_{*}\to-\infty)\) are constants. Importantly, \(Y\) and \(Y'\) are nonoscillatory at both ends, unlike the Teukolsky variable \(R\) or the SN variable \(X\).

Similarly, for the \(Y^{\rm up}\) function, we impose the initial conditions that \(Y^{\rm up}(r = r_+) = {Y^{\rm up}}'(r = r_+) = 0\).⁵ We expand \({Y^{\rm up}}''\) near the horizon as \[\begin{gather} \label{Eq46Yup3939AsymptoticExpansion} {Y^{\rm up}}''(r\to r_{+})=C_{\rm SN}^{\rm inc}\sum_{j=0}^\infty Y_{+,j}^{\rm H}\left(r-r_{+}\right)^j \\ +C_{\rm SN}^{\rm ref}\left(r-r_{+}\right)^{i\frac{(ar_{+}m+2a^2\omega-4r_{+}\omega)}{r_{+}\sqrt{1-a^2}}}\sum_{j=0}^\infty Y_{+,j}^{\rm H}\left(r-r_{+}\right)^j, \end{gather}\tag{30}\] where the coefficients \(Y_{\pm,j}^{\rm H}\) are given in Appendix 8.2. Again, this allows us to start the numerical integration for \(Y^{\rm up}\) at a finite inner boundary (in \(r_{*}\)) since we can analytically integrate Eq. 30 to obtain the proper initial values to use at the inner boundary.

2.2.2 \(\mathcal{W}(r)\) function↩︎

Another ingredient that is needed for the SN formalism is the \(\mathcal{W}\) function, which in turns gives the actual source term for the inhomogeneous SN equation. While we refer readers to Appendix A of Ref. [23] where the general solution of \(\mathcal{W}(r)\) for a generic geodesic was presented, here we rederive these formulas using notations consistent with this paper for the sake of clarity.

The key to deriving the expression for \(\mathcal{W}(r)\) lies in decoupling the terms related to \(\mathcal{N}\) and \(\bar{\mathcal{M}}\), respectively, and decomposing the expression using IBP. This improves the convergence of the integrand. Meanwhile, the expressions outside the integral reflect the asymptotic behavior of \(\mathcal{W}(r)\) at infinity. Integrating the resultant expression inward from infinity then yields the final expression of \(\mathcal{W}(r)\).

In addition, we also need to use two crucial identities. One is associated with \(\mathcal{N}\) and the other with \(\bar{\mathcal{M}}\). In Ref. [23], these derivations and identities are restricted to geodesic motions. Here, we show that these identities and expressions remain valid in all cases including nongeodesics.

When solving for \(\mathcal{W}(r)\) in Eq. 12 [not to be confused with \(W\) defined in Eq. 26 ], we generally partition it into three terms, namely, \[\mathcal{W}(r)=\mathcal{W}_{nn}(r)+\mathcal{W}_{n\bar{m}}(r)+\mathcal{W}_{\bar{m}\bar{m}}(r).\] In particular, we use \(\mathcal{W}_{nn}\) as an example to present part of the derivation. This is because \(\mathcal{W}_{nn}\) will be used in our subsequent analysis of GWs excited by particles falling radially along the spin axis (see Sec. 3.2). The results for \(\mathcal{W}_{n\bar{m}}\) and \(\mathcal{W}_{\bar{m}\bar{m}}\) will be given without detailed derivation.

It is not difficult to show from Eq. 12 and Eq. 13 that \(\mathcal{W}_{nn}\) satisfies the ODE

where \(\mathscr{L}^\dagger_s\equiv\partial_\theta-m/\sin\theta+a\omega\sin\theta+s\cot\theta\) is a differential operator on the angular sector and \[\label{Eq46chi} \chi(r)\equiv\omega t(r)-m\varphi(r)+\int^r\frac{K}{\Delta}d\tilde{r}=\omega v(r)-m\tilde{\varphi}(r),\tag{32}\] with \(v=t+r_{*}\) and \(\tilde{\varphi}=\varphi+\int^r\frac{a}{\Delta}d\tilde{r}\) defined as Kerr ingoing coordinates. Note that Eq. 31 is a more general version of Eq. (A25) in Ref. [23], where there was no \(\sum_{j}^{r=r_j}\) summation in the expression. The summation here is defined such that the particle is located at \(r\) when \(\tau = \tau_1, \tau_2, \cdots, \tau_j\). For unidirectional trajectories, e.g., radial infalls and quasicircular plunges, we have \(j=1\) and the summation can be omitted. While for bound orbits, \(j=\infty\), and for scattering orbits, \(j=2\). Most of the previous works have considered the case when \(j=1\) only, i.e., a particle moves unidirectionally.

For cases where \(j>1\), we should divide the orbit into multi-unidirectional pieces. Divisions are at the turning points where \(u^r=0\). At these turning points, the denominator becomes zero, making the expression singular. However, this does not affect the subsequent integrations, as these singular points of the integrand can be transformed into a smooth form through changing the integration variable. For details, see Ref. [24] for scattering orbits and Ref. [25] for spherical-inclined bound orbits.

Here, we first suppose that the trajectory is unidirectional and omit the summation. By taking the \(r\) derivative of \(\chi(r)\), one can show that \[\chi'(r)=\omega\frac{\mathcal{N}}{u^r}+\left(a\omega\sin^2\theta-m\right)\tilde{\varphi}'.\] Therefore, we can obtain an identity related to \(\mathcal{N}\) that reads \[\label{Eq46N95identity} \begin{align} f(r)\frac{\mathcal{N}}{u^r}e^{i\chi(r)}=&\frac{1}{i\omega}\left\{\left[f(r)e^{i\chi(r)}\right]'\right.\\ &\left.-\left[f'(r)+i\xi(r) f(r)\right]e^{i\chi(r)}\right\}, \end{align}\tag{33}\] where \[\xi(r)=\left(a\omega\sin^2\theta-m\right)\tilde{\varphi}'(r)\sim\mathcal{O}\left(r^{-3/2}\right),\] and \(f(r)\) is an arbitrary smooth function of \(r\). By integrating Eq. 31 and using Eq. 33 twice, we obtain the expression of \(\mathcal{W}_{nn}\) as a three-term form. Similarly, one can obtain the expressions of \(\mathcal{W}_{n\bar{m}}\) and \(\mathcal{W}_{\bar{m}\bar{m}}\) with the help of the identity in Eq. 56 . Schematically, they are

where the expressions of \(f_{0,1,2}\), \(g_{0,1,2}\), and \(h_{0,1,2}\) can be found in Appendix 9.

Moreover, for bound orbits, Eq. 12 actually reads \[\begin{align} \mathcal{W}''=&-\frac{r^2}{\Delta^2}\mathcal{T}\exp\left(i\int^r\frac{K}{\Delta}d\tilde{r}\right)\\ &\times\Theta(r-r_{\rm min})\Theta(r_{\rm max}-r), \end{align}\] where \(\Theta(x)\) is the Heaviside step function, \(r_{\rm min}\) and \(r_{\rm max}\) are the inner and outer edges of the orbit. Thus, we need to multiply all of the \(f\), \(g\), and \(h\) functions in Eqs. 34 by \(\Theta(r-r_{\rm min})\Theta(r_{\rm max}-r)\). This implies that \[\label{Eq46W95bound95infinity} \mathcal{W}(r)=\mathcal{W}'(r)=0\qquad r>r_{\rm max}.\tag{38}\] In Sec. 2.3.1, we will see that this result allows us to discard boundary terms when using the SN-IBP method.

2.2.3 Boundary terms \(Y'\mathcal{W}\) and \(Y\mathcal{W}'\)↩︎

Recall from Eq. 22 that for our SN-IBP approach, we need to evaluate \(Y'(r)\mathcal{W}(r)\) and \(Y(r)\mathcal{W}'(r)\) at both the horizon and infinity, respectively. Here, we will discuss when it is justified to discard these boundary terms.

The general solutions of \(Y(r)\) and \(\mathcal{W}(r)\) can be written as \[\begin{align} Y(r) & =Y^{\rm part}(r)+y_1r+y_0,\\ \mathcal{W}(r) & =\mathcal{W}^{\rm part}(r)+w_1r+w_0, \end{align}\] where \(Y^{\rm part}\) and \(\mathcal{W}^{\rm part}\) are the particular solutions to Eq. 20 and Eq. 12 that we gave earlier in the paper, respectively, and \(y_{0,1}\), \(w_{0,1}\) are some constants. We can choose these constants to our advantages. Moreover, we will refer to the particular solution of \(Y(r)\) obtained in Sec. 2.2.1 [\(Y^{\rm part}(r\to\infty)={Y^{\rm part}}'(r\to\infty)=0\)] with \(y_0=y_1=0\) as the canonical solution of Eq. 20 , and similarly we will refer to the particular solution of \(\mathcal{W}(r)\) obtained in Sec. 2.2.2 with \(w_0=w_1=0\) as the canonical solution of Eq. 12 .

In general, the canonical solution \(\mathcal{W}^{\rm canonical}(r)\) has the asymptotic behaviors \[\mathcal{W}^{\rm canonical}(r)\sim\begin{cases} \mathcal{O}(1), & r\to r_{+}\\ \mathcal{O}(r^{1/2}), & r\to\infty \end{cases},\] and \[{\mathcal{W}^{\rm canonical}}'(r)\sim\begin{cases} \mathcal{O}(1),&r\to r_{+}\\ \mathcal{O}(1),&r\to\infty \end{cases}.\] If one chooses the canonical solution as the particular solution and set \[\label{Eq46w0w1} \begin{align} &w_0=-\mathcal{W}^{\rm part}(r_{+})+r_{+}{\mathcal{W}^{\rm part}}'(r_{+}),\\ &w_1=-{\mathcal{W}^{\rm part}}'(r_{+}), \end{align}\tag{39}\] then the boundary terms of \(\mathcal{W}(r)\) at the horizon vanish [26]. However, this assumes that the limits \(\mathcal{W}^{\rm part}(r\to r_{+})\) and \({\mathcal{W}^{\rm part}}'(r\to r_{+})\) exist. From the analysis in Sec. 2.2.2, this requirement translates into the condition that the limit \(\chi(r\to r_{+})\) exists.

When a particle is close enough to the horizon of a BH, all external forces will be negligible compared to the influence of spacetime curvature itself. This means that the particle moves along a geodesic when approaching the horizon. Therefore, we can use the geodesic equations and obtain \[\frac{d\chi}{dr_{*}}\sim\mathcal{O}(\Delta), \qquad r_{*}\to-\infty.\] As a result, we have \[\chi(r_{*}\to-\infty)={\rm const}.\] With this choice of \(w_{0, 1}\), we have the following asymptotic behaviors for \(\mathcal{W}\) as \[\mathcal{W}(r)\sim\begin{cases} \mathcal{O}(\Delta^2),&r\to r_{+}\\ \mathcal{O}(r),&r\to\infty \end{cases},\] and \[\mathcal{W}'(r)\sim\begin{cases} \mathcal{O}(\Delta),&r\to r_{+}\\ \mathcal{O}(1),&r\to\infty \end{cases},\] at the expense that \(\mathcal{W}(r\to\infty)\) becomes less convergent.⁷

Similarly, the canonical solution \(Y^{\rm canonical}(r)\) has the asymptotic behaviors (which can also be seen in Fig. 1) as \[Y^{\rm canonical}(r) \sim\begin{cases} \mathcal{O}(1), & r\to r_{+}\\ \mathcal{O}(1/r), & r\to\infty \end{cases},\] and \[{Y^{\rm canonical}}'(r) \sim\begin{cases} \mathcal{O}(1),&r\to r_{+}\\ \mathcal{O}(1/r^2),&r\to\infty \end{cases}.\] By setting \[\label{Eq46y0y1} \begin{align} y_0 & =-Y^{\rm part}(r_{+})+r_{+}{Y^{\rm part}}'(r_{+}),\\ y_1 & =-{Y^{\rm part}}'(r_{+}), \end{align}\tag{40}\] we can obtain a solution with the asymptotic behaviors where \[Y(r) \sim \begin{cases} \mathcal{O}(\Delta^2),&r\to r_{+}\\ \mathcal{O}(r),&r\to\infty \end{cases},\] and \[Y'(r) \sim \begin{cases} \mathcal{O}(\Delta),&r\to r_{+}\\ \mathcal{O}(1),&r\to\infty \end{cases}.\]

We summarize the three different approaches for computing \(\psi_4(r \to \infty)\), and in turn gravitational waveforms and fluxes at infinity, that are described in this paper using a flowchart in Fig. 2. The new approach introduced in this subsection is shown in the rightmost column. In the next subsection, we give the recipes to calculate gravitational waveforms and fluxes at infinity using the SN formalism.

2.3 Recipes for calculating waveforms and fluxes at infinity↩︎

2.3.1 Bound orbits↩︎

We have seen from Eq. 38 that the canonical solution of \(\mathcal{W}(r)\) for bound orbits, i.e., \(w_0=w_1=0\), vanishes at infinity. Therefore, we only need to

1. Solve for \(Y(r)\) following the scheme introduced in Sec. 2.2.1.

2. Extract the boundary values of \(Y(r)\) and \(Y'(r)\) at the horizon and calculate \(y_0\) and \(y_1\) following Eqs. 40 .

With these choices, all four of the boundary terms in Eq. 23 vanish. We can then calculate the inhomogeneous solution using Eq. 26 , and therefore the gravitational waveform and fluxes at infinity.

Here, we briefly introduce a procedure for calculating the amplitude for each harmonic of an EMRI waveform on a generic bound geodesic. The derivation is analogous to the one in Ref. [20], but under the SN formalism. Then, in Sec. 3.1, we show some examples of EMRIs waveforms on generic (eccentric-inclined) orbits using our SN-IBP scheme.

A generic bound geodesic orbit in the BL coordinates can be decoupled into harmonics of \(r\) and \(\theta\). This is because the Kerr metric components have no dependence on \(t\) and \(\varphi\). The general solutions to the timelike bound geodesic equation can be expressed as \[\label{Eq46KerrGeoOrbit} \begin{align} &t(\lambda)=\Gamma\lambda+\Delta t[r(\lambda),\theta(\lambda)],\\ &r(\lambda)=\sum_{n=-\infty}^\infty r_ne^{-in\Upsilon_r\lambda},\\ &\theta(\lambda)=\sum_{k=-\infty}^\infty\theta_ke^{-ik\Upsilon_\theta\lambda},\\ &\varphi(\lambda)=\Upsilon_\varphi\lambda+\Delta \varphi[r(\lambda),\theta(\lambda)], \end{align}\tag{41}\] where \(\Gamma\), \(\Upsilon_r\), \(\Upsilon_\theta\), \(\Upsilon_\varphi\) are frequencies parametrized by the Mino time \(\lambda\) which is defined by \(d\tau=\Sigma d\lambda\). To help with our calculations, we also introduce an open source julia package KerrGeodesics.jl for solving timelike Kerr geodesics, see Appendix 10 for details.

Therefore, we can write the Green’s function integral Eq. 26 as \[I=-\mu\int_{\gamma} J_{\ell m\omega}\left[r(\lambda),\theta(\lambda)\right]e^{i(\omega\Gamma-m\Upsilon_\varphi)\lambda}d\lambda.\] The integrand kernel is defined by \[\begin{align} J_{\ell m\omega} & = \frac{d\tau}{d\lambda}\left(W_{nn}\mathcal{N}^2+W_{n\bar{m}}\mathcal{N}\mathcal{M}+W_{\bar{m}\bar{m}}\mathcal{M}^2\right)\\ & = \sum_{k=-\infty}^\infty\sum_{n=-\infty}^\infty J_{\ell mkn}(\omega)e^{-i(k\Upsilon_\theta+n\Upsilon_r)\lambda}, \end{align}\] where \[\label{Eq46J95double95integral} J_{\ell mkn}=\int_0^{2\pi}\int_0^{2\pi}e^{i(k\phi_\theta+n\phi_r)}J_{\ell m\omega}(\phi_r,\phi_\theta)\frac{d\phi_\theta d\phi_r}{(2\pi)^2},\tag{42}\] with \(\phi_r=\Upsilon_r\lambda\), \(\phi_\theta=\Upsilon_\theta\lambda\) defined as the decoupled phases. Finally, we can rewrite the integral, with \(\gamma=(-\infty,\infty)\), as \[\begin{align} I =&\int_{-\infty}^\infty e^{i(\omega\Gamma-m\Upsilon_\varphi-k\Upsilon_\theta-n\Upsilon_r)\lambda}\sum_{k=-\infty}^\infty\sum_{n=-\infty}^\infty J_{\ell mkn}(\omega)d\lambda\\ =&\sum_{k=-\infty}^\infty\sum_{n=-\infty}^\infty2\pi\delta(\omega\Gamma-m\Upsilon_\varphi-k\Upsilon_\theta-n\Upsilon_r) J_{\ell mkn}(\omega). \end{align}\] Then, we insert it into Eq. 18 and Eq. 19 and obtain the gravitational waveform at infinity as \[\begin{align} h&=h_+-ih_\times\\ &=\frac{8}{r}\sum_{\ell m}\int_{-\infty}^\infty\frac{I}{2i\omega B^{\rm inc}_{\rm SN}}{_{-2}}S^{a\omega}_{\ell m}(\theta)e^{-i\omega(t-r_{*})+im\varphi}d\omega\\ &=\sum_{\ell mnk}h_{\ell mnk}, \end{align}\] where \[\omega_{mnk}=m\frac{\Upsilon_{\varphi}}{\Gamma}+n\frac{\Upsilon_r}{\Gamma}+k\frac{\Upsilon_\theta}{\Gamma}\] and \[\label{Eq46hlmnk95bound} h_{\ell mnk}=-\frac{2\mu}{r}\frac{Z_{\ell mnk}^\infty}{\omega_{mnk}^2}{_{-2}}S^{a\omega_{mnk}}_{\ell m}(\theta)e^{-i\omega_{mnk}(t-r_{*})+im\varphi},\tag{43}\] where \[\label{Eq46Zlmnk95bound} Z_{\ell mnk}^{\infty}=-\frac{4i\pi\omega_{mnk}}{B_{\rm SN}^{\rm inc}\Gamma}J_{\ell mnk}.\tag{44}\] The averaged energy flux, angular momentum flux, and Carter constant flux at infinity are given by \[\begin{align} \left\langle \dot{\mathcal{E}} \right\rangle^\infty=&\sum_{\ell mnk}\frac{\left|Z_{\ell mnk}^{\infty}\right|^2}{4\pi\omega_{mnk}^2},\tag{45}\\ \left\langle \dot{\mathcal{L}_{z}} \right\rangle^{\infty} = &\sum_{\ell mnk}\frac{m\left|Z_{\ell mnk}^{\infty}\right|^2}{4\pi\omega_{mnk}^3},\tag{46}\\ \left\langle \dot{\mathcal{Q}} \right\rangle^{\infty} = &\sum_{\ell mnk}\frac{\left(\mathcal{L}_{mnk}+k\Upsilon_\theta\right)\left|Z_{\ell mnk}^{\infty}\right|^2}{2\pi\omega_{mnk}^3},\tag{47} \end{align}\] where \[\begin{align} \mathcal{L}_{mnk} & = m\langle\cot^2\theta\rangle\mathcal{L}_z-a^2\omega_{mnk}\langle\cos^2\theta\rangle\mathcal{E},\\ \langle\cot^2\theta\rangle & =\frac{1}{\pi}\int_0^{\pi}\left[\cot\theta(\phi_\theta)\right]^2d\phi_\theta,\\ \langle\cos^2\theta\rangle & =\frac{1}{\pi}\int_0^{\pi}\left[\cos\theta(\phi_\theta)\right]^2d\phi_\theta. \end{align}\]

2.3.2 Unbound orbits↩︎

Unlike bound orbits, \(\mathcal{W}(r \to \infty)\) and \(\mathcal{W}'(r \to \infty)\) may not vanish for unbound orbits. An example of this would be the radial infall of a particle from infinity. In this case, the natural thing to do with our SN-IBP approach would be choosing \(y_0 = y_1 = 0\) such that the boundary terms at infinity in Eq. 23 vanish, but one still needs to solve Eq. 12 to evaluate \(\mathcal{W}(r = r_+)\) and \(\mathcal{W}(r = r_+)\). While for the original formulation (without using IBP), we also need to solve for \(\mathcal{W}(r)\) and integrate Eq. 17 . Therefore, one needs to solve for \(\mathcal{W}(r)\) either way, and our IBP approach does not have any advantage over the original formulation. A question naturally arises as to which method should be used for unbound orbits? We can answer this question by analyzing the convergence of the integrands in both formulations.

As an example, we derive how one calculates the waveform induced by a particle falling radially into a Kerr BH along its spin axis without using IBP in the original SN formulation. The 4-velocity is given by \[\begin{align} u^t&=\mathcal{E}\frac{r^2+a^2}{\Delta},\\ u^r&=-\frac{\sqrt{\mathcal{E}^2(r^2+a^2)^2-\Delta(r^2+a^2\mathcal{E}^2)}}{r^2+a^2},\\ u^\theta&=u^\varphi=0, \end{align}\] where \(\mathcal{E}\) is the orbital energy per mass. One instantly find that \(\bar{\mathcal{M}}=0\) from the definition in Eq. 27 and therefore only \(\mathcal{W}_{nn}\) is nonvanishing.

Specifically, we see that when \(\mathcal{E}=1\), we have \(u^r(r\to\infty)\sim\mathcal{O}(r^{-1/2})\) and therefore \(\mathcal{W}(r\to\infty)\sim f_0\sim\mathcal{O}(1/r^{1/2})\). When \(\mathcal{E}>1\), we have \(u^r(r\to\infty)\sim\mathcal{O}(1)\) as \(r\to\infty\) and therefore \(\mathcal{W}(r\to\infty)\sim f_0\sim\mathcal{O}(1)\).⁸ By setting \(y_0=y_1=0\), we have \(Y(r\to\infty)\sim\mathcal{O}(1/r)\). From the definition in Eq. 20 and Eq. 12 , we already know that \(Y''(r\to\infty)\sim\mathcal{O}(1/r^3)\), \(\mathcal{W}''(r\to\infty)\sim\mathcal{O}(r^{1/2})\) for \(\mathcal{E}=1\), and \(\mathcal{W}''(r\to\infty)\sim\mathcal{O}(1)\) for \(\mathcal{E}>1\). As a result, we obtain the convergence of the integrands in the SN-IBP and the original SN method, respectively, as \[\label{Eq46IBP47non-IBP95convergence} \begin{align} &\text{\gls{SN}-\gls{IBP}:}\quad Y(r)\mathcal{W}''(r)\sim\begin{cases} \mathcal{O}(1/r^{1/2}),&\mathcal{E}=1\\ \mathcal{O}(1/r),&\mathcal{E}>1 \end{cases},\\ &\text{original \gls{SN}:}\quad Y''(r)\mathcal{W}(r)\sim\begin{cases} \mathcal{O}(1/r^{7/2}),&\mathcal{E}=1\\ \mathcal{O}(1/r^3),&\mathcal{E}>1 \end{cases}. \end{align}\tag{48}\] From Eq. 48 , we conclude that the integrand of the non-IBP method converge way faster than that of the IBP method and suggest using the original SN formulation for unbound orbits. We will show the convergence speed clearly in Sec. 3.2.

For an unbound orbit, there is no discrete frequency spectrum as in the case for a bound orbit. The frequency domain waveform \(\tilde{h}(\omega)\) and the now-continuous energy spectrum \(d\mathcal{E}/d\omega\), in the case of radial infall, can be expressed as \[\label{Eq46amp95spectra95unbound} \begin{align} \tilde{h}_{\ell}(\omega)&=-\frac{2\mu}{r}\frac{Z^\infty_{\ell0\omega}}{\omega^2}{_{-2}}S^{a\omega}_{\ell 0}(\theta)=\frac{8\mu}{r}\frac{X_{\ell 0\omega}^\infty}{c_0}{_{-2}}S^{a\omega}_{\ell 0}(\theta),\\ \left(\frac{d\mathcal{E}}{d\omega}\right)_\ell^\infty &= \frac{\mu^2}{2\omega^2}\left(\left|Z_{\ell 0\omega}^\infty\right|^2+\left|Z_{\ell 0-\omega}^\infty\right|^2\right)\notag\\ & =8\omega^2\mu^2\left(\left|\frac{X_{\ell 0\omega}^\infty}{c_0}\right|^2+\left|\frac{X_{\ell 0-\omega}^\infty}{c_0}\right|^2\right). \end{align}\tag{49}\] The corresponding time-domain waveform is given by \[\label{Eq46Waveform95time95domain95unbound} h_+-ih_\times=\sum_{\ell}\int_{-\infty}^\infty\tilde{h}_\ell(\omega)e^{-i\omega u}d\omega,\tag{50}\] where \(u=t-r_{*}\) is the retarded time.

3 Results↩︎

In this section, we present some example waveforms and energy flux calculations for both bound and unbound orbits. Specifically, for bound orbits, we use the SN-IBP approach introduced in this paper to compute the EMRI waveform snapshot (à la Ref. [20]) for a generic timelike geodesic. For unbound orbits, we consider particles falling radially from infinity along the spin axis and covering two cases—the rest limit (\(\mathcal{E}=1\)) and the ultrarelativistic limit (\(\mathcal{E}\to\infty\)).

3.1 Generic bound stable orbits↩︎

Numerous studies have already calculated the gravitational radiation from particles on bound Kerr geodesic orbits using the Teukolsky formalism, including eccentric-equatorial orbits [27], inclined-spherical orbits [28], and generic orbits [20]. Prior to this work, there were also calculations using the SN formalism on circular-equatorial orbits [29], eccentric-equatorial orbits [30], and inclined-spherical orbits [25]. No calculation for generic orbits has been done with the SN formalism. We present our results and compare them with the literature and codes using the Teukolsky formalism, namely the Teukolsky package from BHPToolkit [31] and pybhpt [32], [33].

Here, we show the results of Eq. 43 –45 . We set \(a=0.9M\), \(p=6M\), \(e=0.7\), \(x=\cos\pi/4\) as our fiducial parameters for a generic geodesic orbit.⁹ For higher values of \(n\) and \(k\), Eq. 42 becomes a highly oscillatory double integral, which is hard to integrate numerically. To achieve a better precision and speed, we employ Levin’s method, which converts a quadrature problem into an ODE problem. The algorithm is introduced in Appendix 11.

To verify our codes, we calculate the energy flux using the SN-IBP method in this work (implemented in GeneralizedSasakiNakamura.jl¹⁰) and pybhpt for the \(\ell=m=2\) and \(\ell=m=4\) modes with \(k=0\) and \(n=0\) to \(n=70\) in Fig. 3.

In addition, we tabulate the total energy flux for each \(\ell\) mode from the two codes, which is defined as \[\label{Eq46flux95l95mode} \left\langle\dot{\mathcal{E}}\right\rangle^\infty_\ell=\sum_{mnk}\left\langle\dot{\mathcal{E}}\right\rangle^\infty_{\ell mnk}.\tag{51}\] The truncation rules¹¹ for the summation in Eq. 51 are specified as follows:

1. For each \(\ell\) mode, we manually set the truncation limits as \(n_{\rm max}^\ell=80+20\ell\) and \(k_{\rm max}^\ell=8+2\ell\).

2. For fixed \(\ell\), \(m\), and \(n=0\), if three consecutive values of \(\langle\dot{\mathcal{E}}\rangle^\infty_{\ell mnk}\) are smaller than \(10^{-6}\times\langle\dot{\mathcal{E}}\rangle^\infty_{\ell}\) (i.e., the current value of the summation of that \(\ell\) mode), then we truncate the \(k\) summation.

3. For fixed \(\ell\), \(m\), and \(k\), if three consecutive values of \(\langle\dot{\mathcal{E}}\rangle^\infty_{\ell mnk}\) are smaller than current \(10^{-6}\times\langle\dot{\mathcal{E}}\rangle^\infty_{\ell}\) (i.e., current value of the summation of this \(\ell\) mode), then we truncate the \(n\) summation.

Following these rules, we calculate the energy fluxes for the \(\ell=2\), \(3\), \(4\), \(5\), and \(6\) modes. These values are tabulated Table 1, together with the total number of modes summed in those calculations.¹² The two sets of numbers agree to the twelve digit, and disagreement only appears after the thirteenth digit (indicated by the brackets in Table 1). Moreover, in Fig. 4, we show the waveform snapshot with the fiducial parameters, using the amplitude data from Table 1. In total, \(5874\) modes were used for the generation of the waveform.

3.2 Radial infalls↩︎

As discussed in Sec. 2.3.2, there are two cases—\(\mathcal{E}=1\), where the particle has no initial velocity at infinity (also referred to as the rest limit) and \(\mathcal{E}>1\). In addition, \(\mathcal{E} \gg 1\) or the ultrarelativistic limit corresponds to a particle moving nearly at the speed of light and hitting a Kerr BH along its spin axis.

The asymptotic behavior of \(u^r\) at infinity differs in these two cases, leading to distinct asymptotic behaviors of \(f_{0,1,2}\) in Eq. 35 . This further results in differences in the asymptotic behavior of \(\mathcal{W}(r)\) and the integrand within the Green’s function integral involved in the calculations. Consequently, the results obtained for these two cases are also significantly different. The behaviors of \(f_{0,1,2}\), \(\mathcal{W}\) and the Green’s function integrand are shown in Fig. 5 for \(\mathcal{E}=1\) and Fig. 6 for \(\mathcal{E}=3\). One can see that \(f_0\sim\mathcal{O}(1/r^{1/2})\), \(f_1\sim\mathcal{O}(1/r^{3/2})\), \(f_2\sim\mathcal{O}(1/r^{3})\) for \(\mathcal{E}=1\) and \(f_0\sim\mathcal{O}(1)\), \(f_1\sim\mathcal{O}(1/r^{2})\), \(f_2\sim\mathcal{O}(1/r^{3})\) for \(\mathcal{E}>1\).

Since the convergence of \(\mathcal{W}\) is controlled by \(f_0\) following Eq. 35 , therefore the overall convergence of the integrand for the original SN formulation [cf. Eq. 17 ] is \(\sim\mathcal{O}(1/r^{7/2})\) for \(\mathcal{E}=1\) and \(\sim\mathcal{O}(1/r^{3})\) for \(\mathcal{E}>1\). However, the integrand in the SN-IBP approach in Eq. 26 behaves as \(\sim\mathcal{O}(1/r^{1/2})\) for \(\mathcal{E}=1\) and \(\sim\mathcal{O}(1/r)\) for \(\mathcal{E}>1\). All of the asymptotic behaviors above agree with our theoretical analyses in Sec. 2.3.2.

Figures 7 and 8 show the amplitude \(\left|X_{\ell0\omega}^\infty/c_0\right|\) and energy spectrum \(\left(d\mathcal{E}/d\omega\right)_\ell^\infty\) in the rest limit and the ultrarelativistic limit (using \(\mathcal{E} = 100\) as an approximation), respectively. We can see that Fig. 7 agrees well with Fig. 1 in Ref. [34], which shows the amplitude and the energy spectrum for a radial infall into a Schwarzchild BH in the rest limit. We also find the same power-law behavior of the amplitudes as \[\left|\frac{X_{\ell0\omega}^\infty}{c_0}\right|\sim \omega^{(\ell-3)/3}\] at the ZFL, i.e. \(\omega\to0\), consistent with the result reported in Ref. [34]. In the ultrarelativistic limit, one can see from Fig. 8 that the power-law behavior of all the \(\ell\) modes are \[\left|\frac{X_{\ell0\omega}^\infty}{c_0}\right|\sim 1/\omega\] at the ZFL. This makes the energy spectra nonvanishing at the ZFL by definition in Eq. 49 .

Therefore, we can extract the value for the energy spectrum values in the ZFL per \(\ell\) mode numerically from our calculations. Theoretically, the total energy spectrum (summed over all \(\ell\) modes) in the ZFL was derived in Ref. [35], which is given by \[\left(\frac{d\mathcal{E}}{d\omega}\right)^{\rm ZFL}=\sum_{\ell=2}^\infty\left(\frac{d\mathcal{E}}{d\omega}\right)_\ell^{\rm ZFL}=\frac{4}{3\pi}\mathcal{E}^2\mu^2,\] while the per-\(\ell\) mode value was also given in Ref. [36] as \[\label{Eq46ZFL95l} \left(\frac{d\mathcal{E}}{d\omega}\right)^{\rm ZFL}_\ell=\frac{4\mathcal{E}^2\mu^2}{\pi}\frac{(2\ell+1)(\ell-2)!}{(\ell+2)!}.\tag{52}\]

Here we show the numerical values extracted from Fig. 8 (b) and the theoretical predictions using Eq. 52 . The results are tabulated in Table 2. Our numerical results match the theoretical results within an error of \(0.1\%\) and are consistent with those shown in Ref. [37].

Finally, we show in Fig. 9 the time-domain waveform in the rest limit by performing an inverse FT in Eq. 50 .

4 Discussions↩︎

4.1 Waveform modeling for extreme mass ratio inspirals↩︎

One obvious application of our SN-IBP approach would be computing EMRI waveforms, which we have already demonstrated in Sec. 3.1 (and Fig. 4). However, such waveform generation requires many—around thousands of—modes to be calculated and summed up, which can take upwards of seconds per waveform and thus too slow for the purpose of LISA data analysis.

Fortunately, the FastEMRIWaveforms framework [38]–[41] solves this problem by generating EMRI waveforms using precomputed waveform amplitude and flux data and thus decouples the waveform generation for data analysis from the relatively expensive waveform calculation. As mentioned in Ref. [41], the framework can be easily extended handle eccentric and inclined orbits around a Kerr BH once the corresponding amplitude and flux data are available, which we can easily generate with the SN-IBP approach.

Since our formalism and code implementation are independent from the Teukolsky formalism, one can also use our code to cross-validate the adiabatic—or 0PA—amplitude and flux data in the literature. In terms of performance, our implementation is comparable with the state-of-the-art pybhpt. We benchmark the performance of three different codes on calculating the waveform amplitude at infinity using the same set of fiducial parameters, i.e., \(a=0.9M\), \(p=6M\), \(e=0.7\), \(x=\cos\pi/4\), that are used throughout the paper, namely, ours, the Teukolsky package from BHPToolkit, and pybhpt. The single-core CPU times are tabulated in Table 3.¹³ Note that no attempt was made to optimize our current implementation, and there is still room for improvement. For example, the current bottleneck of the calculation is actually in solving the homogeneous solutions \(X^{\rm in, up}\). Optimization of our implementation is planned but it is outside the scope of this paper.

4.2 Current limitations and future extensions↩︎

While the IBP approach we presented here drastically simplifies the computation of the waveform amplitude and fluxes at infinity for bound orbits when using the SN formalism, there are still some limitations to our formulation. For instance, the IBP approach has no advantage in computing those quantities near the BH horizon over the original formulation. This is because, when near the horizon (or \(r \to r_{+}\)), the inverse transformation from the SN variable to the Teukolsky variable [cf. Eq. 9 ] is, in fact, dominated by the contribution coming from \(\mathcal{S}\) when using the canonical solution \(\mathcal{W}^{\rm canonical}\) (and by extension \(\mathcal{S}^{\rm canonical}\)). Therefore, we still need to compute \(\mathcal{W}(r = r_{+})\) when computing fluxes down the horizon.

Note that the SN formalism itself is perfectly valid in this case. In fact, one can choose \(w_{0, 1}\) [cf. Eq. 39 in Sec. 2.2.3] such that it is the contribution coming from the \(\Lambda^{-1}\) operator acting on the inhomogeneous SN solution that dominates the transformation [26]. However, we have already imposed the boundary conditions that \(w_{0,1} = 0\) to make the boundary terms in Eq. 23 at infinity vanish. A workaround to this issue is to solve for the inhomogeneous solution with a spin weight of \(s = +2\) instead, which allows us to still use the IBP approach to simplify calculations. We leave this for a separate publication [42].

Another extension to our work here is to consider also generic bound plunge orbits into a Kerr BH. We believe that the IBP approach is still advantageous over the original SN formulation. These kinds of problems also serve as an analytical model for studying and understanding more about the physics and mechanism of the quasinormal mode excitation in binary black hole mergers using the SN formalism [23], [43], [44]. We again leave this for a future publication. An interesting avenue to employ the SN formalism is computing gravitational radiation from scattered orbits around a Kerr BH. This is particularly exciting since the same calculation can also be done with scattering amplitude techniques with post-Minkowskian expansions [45]. The calculation from the BH perturbation theory has only be done in the nonspinning limit (e.g., Ref. [46]). However, due to their unbound nature, we expect that our SN-IBP approach will not be advantageous over the original SN formulation for those scattered orbits (cf. Sec. 2.3.2).

5 Conclusions↩︎

In this work, we introduce a new scheme for solving the inhomogeneous SN equation using integration by parts. When computing gravitational waveforms and fluxes at infinity coming from Kerr BHs perturbed by particles in bound orbits, this simple trick eliminates the need for performing yet another radial integration to obtain the source term that needs to be convolved with a Green’s function as in the original SN formulation. Our approach enables the efficient computation of gravitational waveforms within the SN formalism now in all cases, from bound to unbound orbits, without having to transform between the Teukolsky and SN formalisms in intermediate steps.

Specifically, we define a new auxiliary variable \(Y\) in place of the SN variable \(X\) that we convolve with the source term \(\mathcal{T}\) that one would use in the Teukolsky formalism. This new variable \(Y\) is independent of the source term and therefore only needs to be computed once per frequency. Furthermore, it is nonoscillatory and regular at the BH horizon and spatial infinity, thus allowing for easy numerical calculations. As a byproduct of this work, we also derive a source term for the SN formalism that is valid for arbitrary motion and not just for geodesic motions.

We demonstrate that our approach and code implementation yield waveform amplitude and flux data that are consistent with the literature, while already achieving comparable speed without any optimization attempt. Getting these amplitude and flux data accurately and efficiently is crucial as they enable the rapid generation of waveforms for future LISA data analysis, especially for EMRI waveforms with generic (eccentric and inclined) bound orbits.

Data availability↩︎

The data that support the findings of this article are openly available [47].

6 Deriving the relation between the Teukolsky and Sasaki-Nakamura source terms↩︎

In this Appendix, we rederive the relation between the Teukolsky source term \(\mathcal{T}\) and the SN source term \(\mathcal{S}\) following Refs. [17], [18].

We first consider the variable \(\mathcal{X}\), which is related to the SN variable \(X\) by \(X(r) = \sqrt{\left(r^2 + a^2 \right)/\Delta^{2}} \mathcal{X}\) (cf. Ref. [13]). It satisfies a Regge-Wheeler-like equation given by \[\label{eq:inhomogeneous95RWlikeEqn} \Delta^{2} \left( \frac{1}{\Delta} \mathcal{X}' \right)' - \Delta F_{1} \mathcal{X}' - U_{1} \mathcal{X} = \mathscr{S},\tag{53}\] where \(\mathscr{S}\) is the source term for the \(\mathcal{X}\) variable. Note that this equation reduces to the usual Regge-Wheeler equation when \(a = 0\) since in this case \(\eta(r) = c_0\) is just a constant. We use Eq. 53 to write \(\mathcal{X}''\) in terms of \(\mathcal{X}\) and \(\mathcal{X}'\), which is \[\mathcal{X}'' = \dfrac{\mathscr{S}}{\Delta} + \dfrac{U_1}{\Delta} \mathcal{X} + \left[ F_1 + \dfrac{\Delta'}{\Delta} \right] \mathcal{X}',\] where setting \(\mathscr{S} = 0\) recovers the sourceless case.

Following Refs. [17], [18], we modify the inverse transformation from SN variables \(X\) to Teukolsky variables \(R\) as \[R = \dfrac{1}{\eta} \left[ \left( \alpha + \dfrac{\beta'}{\Delta} \right)\mathcal{X} - \dfrac{\beta}{\Delta} \mathcal{X}' \right] + \dfrac{\mathscr{S}}{\eta},\] and setting \(\mathscr{S} = 0\) recovers the homogeneous case.

We then evaluate \(R'\) in terms of \(\mathcal{X}\) and \(\mathscr{S}\) (and their derivatives) and substitute them back to the inhomogeneous Teukolsky equation in Eq. 1 . As a result, we obtain an ODE for \(\mathscr{S}\), which is given by \[\label{eq:inhomo95Teukolsky95like95eqn} \begin{align} & \Delta^{2} \left[ \dfrac{1}{\Delta} \left( \dfrac{\mathscr{S}}{\eta} \right)' \right]' \\ & + \Delta^{2} \left[ -\dfrac{\beta}{\Delta^{3}} \left(\dfrac{\mathscr{S}}{\eta}\right) \right]' + \left( \alpha - V_{\rm T}\right) \dfrac{\mathscr{S}}{\eta} = -\mathcal{T}. \end{align}\tag{54}\]

Our goal is to solve for \(\mathscr{S}\) in terms of \(\mathcal{T}\). Note that we can rewrite Eq. 54 into a much more compact form as \[\label{eq:compact95form95for95inhomo95Teukolsky95like95eqn} \mathscr{J}^{\dagger} \left[ \mathscr{J}^{\dagger} \left( \dfrac{r^2}{\Delta} \dfrac{\mathscr{S}}{\eta} \right) \right] = - \dfrac{r^2}{\Delta^2} \mathcal{T},\tag{55}\] where \(\mathscr{J}^\dagger \equiv \partial_r+iK/\Delta\) is a differential operator.¹⁴ If we introduce an auxiliary variable \(\mathcal{W}\) such that \[\mathcal{W}(r) = f(r) \exp \left( \int^{r} i\dfrac{K}{\Delta} d\tilde{r} \right),\] for any differentiable function \(f(r)\), then \(\mathcal{W}'\) can be written as \[\label{eq:identity95with95Wprime95and95Jdagger} \mathcal{W}'(r) = \exp \left( \int^{r} i\dfrac{K}{\Delta} d\tilde{r} \right) \mathscr{J}^{\dagger} \left[ f(r) \right].\tag{56}\] If we define \[\mathcal{W}(r) = \dfrac{r^2}{\Delta} \dfrac{\mathscr{S}}{\eta} \exp \left( \int^{r} i\dfrac{K}{\Delta} d\tilde{r} \right),\] then by using the identity in Eq. 56 twice, we have \[\left(\mathcal{W}' \right)' = \exp \left( \int^{r} i\dfrac{K}{\Delta} d\tilde{r} \right) \mathscr{J}^{\dagger} \left[ \mathscr{J}^{\dagger}\left( \dfrac{r^2}{\Delta} \dfrac{\mathcal{S}}{\eta} \right) \right].\] Using Eq. 55 , we have \[} \mathcal{W}'' = -\dfrac{r^2}{\Delta^2} \mathcal{T} \exp \left( \int^{r} i\dfrac{K}{\Delta} d\tilde{r} \right),\] which is the ODE that one needs to solve to obtain \(\mathscr{S}\) from \(\mathcal{T}\).

Given the source term \(\mathscr{S}\) for the variable \(\mathcal{X}\), we can convert that to the source term needed the SN equation \(\mathcal{S}\) simply with \[\label{Eq46SN95source95term} \mathcal{S} = \dfrac{1}{(r^2 + a^2)^{3/2}} \mathscr{S},\tag{57}\] as \(\mathcal{X}\) solutions are related to the corresponding \(X\) solutions by \(X = \sqrt{(r^2 + a^2)/\Delta^{2}} \mathcal{X}\). Putting everything together and the subscript back, we have \[} \mathcal{S}_{\ell m\omega}=\frac{\eta\Delta\mathcal{W}}{(r^2+a^2)^{3/2}r^2}\exp\left(-i\int^r\frac{K}{\Delta}d\tilde{r}\right).\]

7 \(A\) terms and \(W\) terms↩︎

The source term components \(A\) in the Teukolsky formalism are given by \[\begin{align} &A_{nn0}=\frac{\mathscr{A}}{2}\rho\bar{\rho}^2\mathcal{N}^2\mathscr{L}_1^\dagger\left[\rho^{-4}\mathscr{L}_2^\dagger\left(\rho^3S\right)\right],\\ &A_{n\bar{m}0}=\mathscr{A}\bar{\rho}^2\mathcal{N}\bar{\mathcal{M}}\left[\left(\mathscr{L}_2^{\dagger}S\right)\left(\frac{i K}{\Delta}-\rho-\bar{\rho}\right)\right.\\ &\qquad\qquad\qquad\qquad\qquad \left.-a\sin\theta S\frac{K}{\Delta}\left(\rho-\bar{\rho}\right)\right],\notag\\ &A_{\bar{m}\bar{m}0}=\frac{\mathscr{A}}{2}\bar{\rho}^2\bar{\mathcal{M}}^2S\left[-i\left(\frac{K}{\Delta}\right)_{,r}-\frac{K^2}{\Delta^2}-2i\rho\frac{K}{\Delta}\right],\\ &A_{n\bar{m}1}=\mathscr{A}\bar{\rho}^2\mathcal{N}\bar{\mathcal{M}}\left[\mathscr{L}_2^{\dagger}S+i a\sin\theta\left(\rho-\bar{\rho}\right)S\right],\\ &A_{\bar{m}\bar{m}1}=\mathscr{A}\bar{\rho}^2\bar{\mathcal{M}}^2S\left(i\frac{K}{\Delta}-\rho\right),\\ &A_{\bar{m}\bar{m}2}=\frac{\mathscr{A}}{2}\bar{\rho}^2\bar{\mathcal{M}}^2S, \end{align}\] where \(S\) are the SWSH with all of its subscripts suppressed to avoid confusion.

The source term components in the SN formalism are given by \[\begin{align} W_{nn}=&\mathscr{A}\frac{\rho\bar{\rho}^2}{2}\mathscr{L}_1^\dagger\left[\rho^{-4}\mathscr{L}_2^\dagger\left(\rho^3 S\right)\right]r^2Y\mathrm{phase},\\ W_{n\bar{m}}=&-\mathscr{A}r\bar{\rho}^2\left\{\left(\mathscr{L}_2^\dagger S\right)\left(\rho+\bar{\rho}\right)rY\right.\\ &\left.+\left[\mathscr{L}_2^\dagger S+i a\sin\theta\left(\rho-\bar{\rho}\right)S\right]\left(2Y+rY'\right)\right\}\mathrm{phase},\notag\\ W_{\bar{m}\bar{m}}=&\mathscr{A}S\bar{\rho}^2\biggl[\frac{X}{2\sqrt{r^2+a^2}}+\left(Y+2rY'\right)\mathrm{phase}\\ &+\rho r\left(2Y+rY'\right)\mathrm{phase}\biggr],\notag\\ \text{phase}=&\exp\left(i \int^r\frac{K}{\Delta}d\tilde{r}\right) \label{eq:KoverDeltaIntegral}\\ =&\exp\left(i \omega r_{*}-\frac{i am}{2\sqrt{1-a^2}}\ln\frac{r-r_{+}}{r-r_{-}}\right)\notag. \end{align}\tag{58}\]

7.1 Normalization conventions↩︎

The value of the constant \(\mathscr{A}\) above depends on the normalization conventions on the FT and SWSHs adopted. Specifically for the FT, there are canonically two conventions for the normalization, namely the unitary FT convention where \[\begin{align} F(\omega)&=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^\infty f(t)e^{i\omega t}dt,\\ f(t)&=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^\infty F(\omega)e^{-i\omega t}d\omega, \end{align}\] and the non-unitary FT convention where \[\begin{align} F(\omega) & =\frac{1}{2\pi}\int_{-\infty}^\infty f(t)e^{i\omega t}dt,\\ f(t) & =\int_{-\infty}^\infty F(\omega)e^{-i\omega t}d\omega. \end{align}\] As for SWSHs, there are also two normalization conventions where either \[\label{eq:SWSH95normalization95scheme1} \int_0^\pi \left|{_s}S_{\ell m}^{a\omega}(\theta)\right|^2\sin\theta d\theta=1,\tag{59}\] which we will refer to as the SWSH normalization scheme 1, and \[\label{eq:SWSH95normalization95scheme2} \int_0^\pi \left|{_s}S_{\ell m}^{a\omega}(\theta)\right|^2\sin\theta d\theta=\frac{1}{2\pi},\tag{60}\] which we will refer to as the SWSH normalization scheme 2, respectively.

Table 4 shows the value of \(\mathscr{A}\) with different choices of normalization conventions. Although the exact choice does not have any impact on physics, care should be taken when comparing results from different papers since the value of \(\mathscr{A}\) may vary literature to literature. In this work, we use the nonunitary FT convention and the SWSH normalization scheme 2, and therefore \(\mathscr{A}=-1\).

8 Asymptotic expansions for \(Y^{\rm in, up}\)↩︎

In this Appendix, we derive the asymptotic expansions for \(Y^{\rm in, up}\), respectively, for speeding up the numerical integration of Eq. 20 , which is repeated here for reference as \[} Y_{\ell m\omega}^{\rm in/up\;\prime\prime}(r)\equiv \frac{X_{\ell m\omega}^{\rm in/up}(r)}{r^2\sqrt{r^2+a^2}}\exp\left(-i\int^r\frac{K}{\Delta}dr\right).\]

8.1 \(Y^{\rm in}\) solution↩︎

Recall that in Ref. [13], we have shown that asymptotically as \(r \to \infty\), \[X^{\rm in}(r\to\infty)=B_{\rm SN}^{\rm ref}e^{i\omega r_{*}}\sum_{w=0}^\infty\frac{\mathcal{C}^\infty_{+,w}}{r^w}+B_{\rm SN}^{\rm inc}e^{-i\omega r_*}\sum_{w=0}^\infty\frac{\mathcal{C}^\infty_{-,w}}{r^w}.\] The expressions of \(\mathcal{C}^{\infty}_{\pm,1,2,3}\) can be found in the Appendix G of Ref. [13] (note that \(\mathcal{C}^{\infty}_{\pm, 0} = 1\)). To write down an asymptotic expansion of \(Y^{\rm in}(r \to \infty)\), we also need the series expansions of the two other terms, which are given by \[\begin{align} \exp\left(-i\int^r\frac{K}{\Delta}d\tilde{r}\right) & =e^{-i\omega r_{*}}\sum_{j=0}^\infty\frac{a_j}{r^j},\\ \frac{1}{r^2\sqrt{r^2+a^2}} & =\frac{1}{r^3}\sum_{j=0}^\infty\frac{b_j}{r^j}, \end{align}\] where \[\begin{align} a_j & = \frac{1}{j!}B_j( P_1,\dots,P_j),\\ P_j & = \frac{iam\left(r_{+}^j-r_{-}^j\right)\Gamma(j)}{r_{-}-r_{+}},\\ b_j & = \frac{1+(-1)^j}{2} a^j \binom{-1/2}{j/2}, \end{align}\] and \(B_j\) denotes the \(j\)th complete exponential Bell polynomial. Therefore, the piece that is proportional to \(B_{\rm SN}^{\rm ref}\), which we denote as \({Y^{\rm in}_+}''(r\to\infty)\), can be expressed as \[\label{Eq46Y393995inf9543} \begin{align} &{Y^{\rm in}_+}''(r\to\infty)\\ =&\frac{X^\infty_+(r)}{r^2\sqrt{r^2+a^2}}\exp\left(-i\int^r\frac{K}{\Delta}d\tilde{r}\right)\\ =&\frac{B_{\rm SN}^{\rm ref}}{r^3}\left(\sum_{j=0}^\infty\frac{a_j}{r^i}\right) \left(\sum_{v=0}^\infty\frac{b_v}{r^j}\right) \left(\sum_{w=0}^\infty\frac{\mathcal{C}^\infty_{+,w}}{r^k}\right)\\ =&B_{\rm SN}^{\rm ref}\sum_{j=0}^\infty\frac{Y^{\infty,+}_{j}}{r^{j+3}}, \end{align}\tag{61}\] where \[Y^{\infty}_{+,j}=\sum_{v=0}^j\sum_{w=0}^{j-v}a_{v} b_{w} \mathcal{C}^{\infty}_{+,j-v-w}.\] Notice that it is not oscillatory because \(e^{-i \int^r\frac{K}{\Delta}d\tilde{r}} \sim e^{-i\omega r_{*}}\) [cf. Eq. 58 ] cancels out the phase term \(e^{i \omega r_{*}}\) coming from \(X^{\rm in}\).

The other piece that is proportional to \(B_{\rm SN}^{\rm inc}\), which we denote as \({Y^{\rm in}_-}''(r\to\infty)\), is more complicated because the phase terms do not cancel out each other. Here, we need to expand also \(r_{*}(r \to \infty)\), which is given by \[r_*=r+2\ln \frac{r}{2}-\sum_{v=1}^\infty\frac{2}{vr^v}\left(\sum_{j=0}^vr_+^jr_-^{v-j}\right).\] Therefore, we have \[e^{-2i\omega r_{*}}=e^{4i\omega\ln 2}\frac{e^{-2i\omega r}}{r^{4i\omega}}\sum_{j=0}^\infty\frac{d_j}{r^j},\] where \[\begin{align} &d_j=\frac{1}{j!}B_j(Q_1,\dots,Q_j),\\ &Q_j=4i\omega\Gamma(j)\left(\sum_{v=0}^j r_+^vr_-^{j-v}\right). \end{align}\] Finally, we have \[\label{Eq46Y393995inf95-} \begin{align} &{Y_-^{\rm in}}''(r\to\infty)\\ =&B_{\rm SN}^{\rm inc}\frac{e^{4i\omega\ln 2-2i\omega r}}{r^{3+4i\omega}}\left(\sum_{j=0}^\infty\frac{a_j}{r^j}\right) \left(\sum_{w=0}^\infty\frac{b_w}{r^w}\right)\\ & \left(\sum_{v=0}^\infty\frac{\mathcal{C}^\infty_{-,v}}{r^v}\right) \left(\sum_{u=0}^\infty\frac{d_u}{r^u}\right)\\ =&B_{\rm SN}^{\rm inc}\frac{e^{4i\omega\ln 2-2i\omega r}}{r^{4i\omega}}\sum_{j=0}^\infty\frac{Y_{-,j}^{\infty}}{r^{j+3}}, \end{align}\tag{62}\] where \[Y_{-,j}^{\infty}=\sum_{v=0}^j\sum_{w=0}^{j-v}\sum_{u=0}^{j-v-w}a_v b_w \mathcal{C}^\infty_{-,u} d_{j-v-w-u}.\] Combining Eq. 61 and Eq. 62 , we obtain Eq. 28 . It is not difficult to show that \(Y_{+,0}^{\infty}=Y_{-,0}^{\infty}=1\). Here, we also give the explicit expressions of the next three coefficients, which are given by

The initial values of \(Y^{\rm in}(r)\) and \({Y^{\rm in}}'(r)\) at large \(r = r_{\rm out}\) can then be obtained by integrating Eq. 28 . The \({Y_+^{\rm in}}\) piece is straightforward and is given by \[\begin{align} &{Y_+^{\rm in}}'(r_{\rm out})=-B_{\rm SN}^{\rm ref}\sum_{j=0}^\infty\frac{Y_{+,j}^{\infty}}{j+2}\frac{1}{r_{\rm out}^{j+2}},\\ &Y_+^{\rm in}(r_{\rm out})=B_{\rm SN}^{\rm ref}\sum_{j=0}^\infty\frac{Y_{+,j}^{\infty}}{(j+1)(j+2)}\frac{1}{r_{\rm out}^{j+1}}. \end{align}\]

While the \(Y_-^{\rm in}\) piece is more complicated because the phase term is nonvanishing. We define

where \(_1F_2(a;b,c;x)\) is the hypergeometric function. With these, we can write \[\begin{align} {Y_-^{\rm in}}'(r_{\rm out})=&-e^{4i\omega \ln 2}B_{\rm SN}^{\rm inc}\sum_{j=0}^{\infty}Y_{-,j}^{\infty}y_{j+3}(r_{\rm out}),\\ Y_-^{\rm in}(r_{\rm out})=&e^{4i\omega \ln 2}B_{\rm SN}^{\rm inc}\sum_{j=0}^{\infty}Y_{-,j}^{\infty}\left[y_{j+2}(r_{\rm out})\right.\notag\\ &\left.-r_{\rm out}\cdot y_{j+3}(r_{\rm out})\right]. \end{align}\]

Unfortunately, using the hypergeometric function implemented in HypergeometricFunctions.jl [48] is too time-consuming for an acceptable precision when \(\omega r_{\rm out}\) is a relatively large value. Therefore, we switch to an asymptotic expansion given by \[\begin{gather} \frac{\Gamma(a_1)}{\Gamma(b_1)\Gamma(b_2)}{}_1F_2(a_1;b_1,b_2;-z)\\ = {}_1H_2(z)+{}_1E_2(ze^{-\pi i})+{}_1E_2(ze^{\pi i}), \end{gather}\] where \[\begin{align} {}_1H_2(z)=&\sum_{j=0}^\infty\frac{(-1)^j}{j!}\frac{\Gamma(a_1+j)}{\Gamma(b_1-a_1-j)\Gamma(b_2-a_1-j)}z^{-a_1-j},\\ {}_1E_2(z)=&\frac{e^{2\sqrt{z}}}{\sqrt{\pi}}\sum_{j=0}^\infty c_k\frac{z^{(\nu-j)/2}}{2^{j+1}}\\ =&\begin{cases} \frac{(-z)^{\nu/2}e^{-2i\sqrt{-z}}}{\sqrt{\pi}}\sum_{j=0}^\infty\frac{c_k(-z)^{-j/2}}{2^{j+1}}& z\to ze^{-\pi i}\\ \frac{(-z)^{\nu/2}e^{2i\sqrt{-z}}}{\sqrt{\pi}}\sum_{j=0}^\infty\frac{c_k(-z)^{-j/2}}{2^{j+1}}& z\to ze^{\pi i} \end{cases},\notag \end{align}\] with \[\begin{align} \nu=&a_1-b_1-b_2+\frac{1}{2},\\ c_0=&1,\\ c_j=&-\frac{1}{4j}\sum_{w=0}^{j-1}c_w e_{j,w},\\ e_{j,w}=&\frac{(1-\nu-2b_1+w)_{2+j-w}(a_1-b_1)}{(b_2-b_1)(1-b_1)}\notag\\ &+\frac{(1-\nu-2b_2+w)_{2+j-w}(a_1-b_2)}{(b_1-b_2)(1-b_2)}\\ &+\frac{(w-1-\nu)_{2+j-w}(a_1-1)}{(1-b_1)(1-b_2)}\notag. \end{align}\] In our case, we have \(b_2=a_1+1\). Therefore, \[\frac{\Gamma(a_1)}{\Gamma(b_1)\Gamma(b_2)}=\frac{1}{a_1\Gamma(b_1)}.\]

Now we can construct the initial conditions for the ODE in Eq. 29 for \(Y^{\rm in}\) and solve it inward to the horizon to get the values of \(Y^{\rm in}\) and \({Y^{\rm in}}'\) across the entire domain of definition with \[\begin{align} \left.Y\right|_{r_{*}=r_{*}^{\rm out}}& =Y^{\rm in}_+(r_{\rm out})+Y^{\rm in}_-(r_{\rm out}),\\ \left.\frac{dY}{dr_{*}}\right|_{r_{*}=r_{*}^{\rm out}} & =\frac{\Delta}{r_{\rm out}^2+a^2}\left[{Y^{\rm in}_+}'(r_{\rm out})+{Y^{\rm in}_-}'(r_{\rm out})\right]. \end{align}\]

Figure 10 shows the relative error between the asymptotic expansion of \({Y^{\rm in}}''(r\to\infty)\) defined in Eq. 28 and its definition in Eq. 20 by expanding up to \(\sim\mathcal{O}(1/r^6)\) order. We can see that for larger values of \(\omega\), the asymptotic expansion converges rapidly to the definition. However, for smaller values of \(\omega\), the convergence decreases, and we need to have a larger \(r_{*}^{\rm out}\), or equivalently, increase the expansion order. From our calculations, we find that generally setting \(r_{*}^{\rm out}={\rm max}(1000, 10\pi/|\omega|)\) is sufficient to reach the \(10^{-12}\) relative tolerance if we truncate the expansion at \(\sim\mathcal{O}(1/r^6)\) order.

8.2 \(Y^{\rm up}\) solution↩︎

Recall that in Ref. [13], we derived the asymptotic expansion of \(X^{\rm up}\) for \(r \to r_+\) as \[\begin{gather} X^{\rm up}(r\to r_{+}) = C_{\rm SN}^{\rm inc}e^{i\kappa r_{*}}\sum_{w=0}^\infty\mathcal{C}^{\rm H}_{+,w}\left(r-r_{+}\right)^w\\ +C_{\rm SN}^{\rm ref}e^{-i\kappa r_{*}}\sum_{w=0}^\infty\mathcal{C}^{\rm H}_{-,w}\left(r-r_{+}\right)^w. \end{gather}\] Following the same procedure in Appendix 8.1, we obtain \[\tag{63} \begin{align} {Y^{\rm up}_+}''(r\to r_{+}) & = C_{\rm SN}^{\rm inc}\sum_{j=0}^\infty Y^{\rm H}_{+,j}\left(r-r_{+}\right)^j,\tag{64}\\ {Y^{\rm up}_-}''(r\to r_{+}) & = C_{\rm SN}^{\rm ref}\sum_{j=0}^\infty Y^{\rm H}_{-,j}\left(r-r_{+}\right)^{j+iq},\tag{65} \end{align}\] where \[q=\frac{(ar_{+}m+2a^2\omega-4r_{+}\omega)}{r_{+}\sqrt{1-a^2}}.\] Unfortunately, the expressions of \(Y^{\rm H}_{\pm,0,1,2}\) and \(\mathcal{C}^{\rm H}_{\pm, 0,1,2}\) are too long to show directly here and are available in a Mathematica notebook [47]. By combining Eq. 64 and Eq. 65 , we obtain Eq. 30 .

We then integrate Eqs. 63 to get the initial values as \[\begin{align} &{Y_+^{\rm up}}'(r_{\rm in})=C_{\rm SN}^{\rm inc}\sum_{j=0}^\infty Y_{+,j}^{\rm H}\frac{\left(r_{\rm in}-r_{+}\right)^{j+1}}{j+1},\\ &Y_+^{\rm up}(r_{\rm in})=C_{\rm SN}^{\rm inc}\sum_{j=0}^\infty Y_{+,j}^{\rm H}\frac{\left(r_{\rm in}-r_{+}\right)^{j+2}}{(j+1)(j+2)},\\ &{Y_-^{\rm up}}'(r_{\rm in})=C_{\rm SN}^{\rm ref}\sum_{j=0}^\infty Y_{-,j}^{\rm H}\frac{\left(r_{\rm in}-r_{+}\right)^{j+1+iq}}{j+1+iq},\\ &Y_-^{\rm up}(r_{\rm in})=C_{\rm SN}^{\rm ref}\sum_{j=0}^\infty Y_{-,j}^{\rm H}\frac{\left(r_{\rm in}-r_{+}\right)^{j+2+iq}}{(j+1+iq)(j+2+iq)}. \end{align}\]

Now we can also construct the initial conditions for the ODE in Eq. 29 for \(Y^{\rm up}\) and solve it outward to infinity to get the initial values of \(Y^{\rm up}(r)\) and \({Y^{\rm up}}'(r)\) across the entire domain of definition using \[\begin{align} \left. Y^{\rm up} \right|_{r_{*}=r_{*}^{\rm in}} & =Y^{\rm up}_+(r_{\rm in})+Y^{\rm up}_-(r_{\rm in}),\\ \left. \frac{dY^{\rm up}}{dr_{*}}\right|_{r_{*}=r_{*}^{\rm in}} & =\frac{\Delta}{r_{\rm in}^2+a^2}\left[{Y^{\rm up}_+}'(r_{\rm in})+{Y^{\rm up}_-}'(r_{\rm in})\right]. \end{align}\]

9 Expressions in \(\mathcal{W}\) integrals↩︎

with \[\begin{align} &w_{nn}^{(0)}(r)=\frac{1}{2}r^2\rho\bar{\rho}^2u^r \mathscr{L}_1^\dagger\left[\rho^{-4}\mathscr{L}_2^\dagger\left(\rho^3S\right)\right],\\ &w_{nn}^{(1)}(r)=w_{nn}^{(0)}(r)\left(\frac{\mathcal{N}}{u^r}\right)'\frac{u^r}{\mathcal{N}}+{w_{nn}^{(0)}}'(r)+i\xi(r)w_{nn}^{(0)}(r),\\ &w_{n\bar{m}}^{(0)}(r)=\frac{r^2\bar{\rho}^3}{\rho}\bar{\mathcal{M}}\left[\mathscr{L}_2^\dagger S+ia\left(\rho-\bar{\rho}\right)\sin\theta S\right],\\ &w_{n\bar{m}}^{(1)}(r)=\frac{r^2\bar{\rho}\bar{\mathcal{M}}}{2}\mathscr{L}_2^\dagger\left[\rho^3S\left(\bar{\rho}^2\rho^{-4}\right)'\right],\\ &w_{n\bar{m}}^{(2)}(r)=\bar{\rho}\bar{\mathcal{M}}\left\{\frac{r^2\bar{\rho}^2}{\rho}\left[\mathscr{L}_2^\dagger S+i a\left(\rho-\bar{\rho}\right)\sin\theta S\right]\right\}'. \end{align}\]

10 Solving for geodesic motions for Kerr black holes↩︎

The motion of a particle in a Kerr background is determined by the four constants of motion, namely, the mass \(\mu\), energy \(E\), angular momentum along the spin axis \(L_z\), and Carter constant \(Q\). In test mass limit, i.e., \(\mu\ll1\), the motion can be described by the following equations of motion in Kerr spacetime: \[\label{Eq464velocities} \begin{align} &\Sigma\frac{dt}{d\tau}=-a\left(a\mathcal{E}\sin^2\theta-\mathcal{L}_z\right)+\frac{r^2+a^2}{\Delta}P,\\ &\Sigma\frac{dr}{d\tau}=\pm\sqrt{R},\\ &\Sigma\frac{d\theta}{d\tau}=\pm\sqrt{\Theta},\\ &\Sigma\frac{d\varphi}{d\tau}=-\left(a\mathcal{E}-\frac{\mathcal{L}_z}{\sin^2\theta}\right)+\frac{a}{\Delta}P, \end{align}\tag{67}\] where \[\label{Eq464veloFunctions} \begin{align} &P=\mathcal{E}(r^2+a^2)-a\mathcal{L}_z,\\ &R=P^2-\Delta\left[r^2+\left(\mathcal{L}_z-a\mathcal{E}\right)^2+\mathcal{Q}\right],\\ &\Theta=\mathcal{Q}-\cos^2\theta\left[a^2\left(1-\mathcal{E}^2\right)+\frac{\mathcal{L}_z^2}{\sin^2\theta}\right]. \end{align}\tag{68}\] Here the constants are rescaled by \(\mathcal{E} \equiv E/\mu\), \(\mathcal{L}_z \equiv L_z/\left(M\mu\right)\), and \(\mathcal{Q} \equiv Q/\left(M\mu\right)^2\).

We follow the procedure in Ref. [49] to solve the equations. Here we outline the algorithm:

1. For a given set of orbital parameters, namely, the spin parameter \(a\), semi-latus rectum \(p\), eccentricity \(e\), and inclination parameter \(x \equiv \cos\theta_{\rm inc}\), we follow Ref. [50] to map them to the constants of motion \((\mathcal{E},\;\mathcal{L}_z,\;\mathcal{Q})\).

2. Then we calculate the main frequencies of the motions by integrating the geodesic equations. By using the Mino time \(\lambda\) where \(d\lambda \equiv d\tau/\Sigma\), we can decouple the \(r\)- and \(\theta\)-direction motions [51] and obtain the Mino frequencies \(\Upsilon_r\) and \(\Upsilon_\theta\) using elliptic integrals. We calculate the \(t\)- and \(\varphi\)-direction Mino frequencies \(\Gamma\) and \(\Upsilon_\varphi\), respectively, based on the \(r\) and \(\theta\) motions.

3. Finally, we solve the geodesic equations by integrating them over one period, analytically expressing them as elliptic integrals using the Mino time \(\lambda\), and extending them to full domain from \(-\infty\) to \(\infty\).

Following the above three steps, we obtain the solution of a generic timelike bound geodesic motion which can be written as Eqs. 41 .

We implemented KerrGeodesics.jl¹⁵, utilizing the high performance of julia to numerically compute the elliptic integrals (using the package Elliptic.jl). One can obtain all the ingredients for calculating the solution in microseconds. Figure 11 shows the trajectory with \(a=0.9M\), \(p=6M\), \(e=0.7\), \(x=\cos\pi/4\), which corresponds to the waveform in Fig. 4.

Currently, unbound geodesics are not available in the package. We plan to include plunge orbits (see Ref. [52]) and scattering orbits (see Ref. [53]) in the future.

11 Levin’s method↩︎

The evaluation of highly oscillatory integrals, such as Eq. 42 when \(n\) and \(k\) are large, is difficult and often encountered in a wide range of problems. To tackle this issue, we employ Levin’s method [54], which converts the quadrature problem into an equivalent system of ODEs that gives the antiderivative function of the integrand kernel. In the following, we briefly illustrate Levin’s method for one-dimensional integrals.

For a one-dimensional integral of the form \[\label{eq:Levin951D95integral} \mathbb{I}=\int_a^b f(r)e^{ig(r)}dr,\tag{69}\] where the phase function \(g(r)\) varies rapidly while the kernel function \(f(r)\) varies slowly, we want to find the solution of \(p(r)\) that satisfies the following ODE \[\label{eq:ODE95Levin} p'(r)+ig'(r)p(r)=f(r).\tag{70}\] With \(p(r)\), \(\mathbb{I}\) can be evaluated using simply \[\mathbb{I}=p(b)e^{ig(b)}-p(a)e^{ig(a)}.\] Following Ref. [55], we solve Eq. 70 for \(p(r)\) using a Chebyshev spectral method. The ODE problem is further transformed into a problem of solving a system of linear equations given by \[\left[\overleftrightarrow{D}+i\overleftrightarrow{g'}\right] \; \vec{p}=\vec{f},\] where \(\overleftrightarrow{D}\) is the differentiation matrix, \(\overleftrightarrow{g'}\) and \(\vec{f}\) are a diagonal matrix and a vector evaluated at the collocation points, respectively. We refer readers to a detailed exposition of the algorithm for one-dimensional integrals and the two-dimensional generalization in Refs. [55] and [56], respectively.

To facilitate the calculations in this work, we implemented an optimized version of the adaptive Levin’s algorithm following Refs. [55], [56] in julia, which is publicly available as AdaptiveLevin.jl¹⁶.

References↩︎