Fast evolution and waveform generator for extreme-mass-ratio inspirals in equatorial-circular orbits

An extreme-mass-ratio inspiral (EMRI) arises following the capture of a compact star with stellar mass (white dwarf, neutron star or black hole) by a supermassive black hole. The orbital radius of such an EMRI is less than $O({10}^{1})$ the Schwarzschild radius of the supermassive black hole. Because of gravitational radiation, the orbit of the small body shrinks toward the central black hole in a long timescale. EMRIs are very important for revealing the properties of supermassive black holes since people may observe the gravitational signals with many wave cycles from the region near the horizon of the black hole. EMRIs are potential sources for the evolved Laser Interferometer Space Antenna (eLISA), a space gravitational wave observatory supported by the European Space Agency [1]. A pathfinder will be launched in 2015 to pave the way for the eLISA.

Due to the match-filter technology employed in gravitational wave detection, people must have a huge amount of theoretical waveform templates with enough accuracy in a very large parameter space. Up to now, there are usually three methods to compute the theoretical waveforms: the first one is post-Newtonian (PN) approximation, the second one is numerical relativity, and the last one is the black hole perturbation theory. As analyzed in a great deal of literature (for example see [2]), for EMRIs, the PN expansion would lose accuracy greatly in such a highly relativistic region, and even the factorized-resummed PN waveforms in an effective-one-body (EOB—an analytical approach which aims at providing an accurate description of the motion and radiation of coalescing binary black holes with arbitrary mass ratio, see review [3] for details) frame also do not have good performance for spinning black holes (see for example [4] and references therein). Numerical relativity still cannot handle the binary black hole systems with extreme mass ratio. A recent paper documents the simulation of a case of mass ratio 1:100 (without spin) with only two orbits' evolution [5].

Therefore, for such small mass-ratio binaries, the black hole perturbation theory is a good tool to study EMRIs (intermediate-mass-ratio inspirals; mass ratio $\sim {10}^{-7}-{10}^{-4}$) [6–14, 16, 17] and even IMRIs (mass ratio $\sim {10}^{-3}-{10}^{-2}$) [15, 18–22]. It means that the small body can be treated as a perturbation of the background field of the central supermassive black hole. The black hole perturbation theory was built by Regge, Wheeler and Zerilli in Schwarzschild spacetime [23, 24], and by Teukolsky in Kerr background [25, 26].

Usually, the eLISA observes wave-signals over a span of several years. This means that one needs to model the waveforms of EMRIs over a length of several years. Though the cost of calculating the Teukolsky equation is much less than numerical relativity, it is still unaffordable for evolving $O({10}^{5}-{10}^{6})$ orbits in various parameters. Therefore, researchers are still challenged about how to greatly reduce the CPU time of the simulation of EMRI waveforms, while at the same time maintaining enough accuracy. For the quasi-circular orbits, Yunes et al used Teukolsky-based energy fluxes to fit higher-order PN fluxes [27, 28]. Using the self-force data of the Schwarzschild black hole, Lackeos and Burko fitted polynomials with PN expressions to calculate the orbital evolution of IMRIs, but they still numerically solved the Teukolsky equation to get the waveforms [29]. Fujita gave out the expressions of gravitational radiation up to the 14th PN order by computing the Teukolsky equation analytically [33]. However, in this paper, we completely abandon using the analytical PN expansions. We use the Teukolsky-based numerical data to directly fit out a set of polynomials for energy fluxes and waveforms. This method is very simple and efficient, and the fitted polynomials can give very highly-accurate energy fluxes and waveforms. In principle, these polynomials are also a kind of PN expansion, but all the coefficients of such 'PN expansions' are obtained numerically from the fitting of Teukolsky-based data.

In the next section, we introduce our EOB-Teukolsky (ET) codes. The details of our method for fitting polynomials are presented in section 3. In section 4, results and comparisons are shown. Finally, conclusions and remarks are given in section 5. Throughout the paper, we use units $G=c=1$ and the metric signature $(-,+,+,+)$. Distance and time are measured by the central black hole mass M.

Our ET codes include two main parts: one is the EOB dynamics driver, and the other is the Teukolsky equation solver. The EOB part gives the orbital parameters to the Teukolsky equation, and then the later one produces waveforms and energy fluxes. Next the Teukolsky-based energy fluxes source the EOB dynamics to drive the orbital evolution of the small body around the supermassive black hole [21, 34]. A detailed introduction of the EOB dynamics will be presented in section 4.

The Teukolsky equation in ET codes can be solved in the frequency domain for the inspiralling phase and in the time domain for plunge and merge states. In the present work, because we focus on only the inspiralling process, we just need the frequency-domain Teukolsky calculation. We employ a semi-analytical method to solve the frequency-domain equation which was developed in [36–39] to replace the previous numerical integration method [40].

After decomposing the Weyl curvature (complex) scalar ${\psi }_{4}$ in a Fourier series [25],

$$\begin{eqnarray}&&{\psi }_{4}={\rho }^{4}{\displaystyle \int }_{-\infty }^{+\infty }{\rm{d}}\omega \displaystyle \sum _{{lm}}{R}_{{lm}\omega }{(r)}_{-2}{S}_{{lm}}^{a\omega }(\theta ){e}^{{\rm{i}}m\phi }{e}^{-{\rm{i}}\omega t},\end{eqnarray} \tag{ 1 }$$

where $\rho =-1/(r-{ia}\mathrm{cos}\theta )$, the Teukolsky equation is divided into two parts. One is the radial master equation

$$\begin{eqnarray}&&{{\rm{\Delta }}}^{2}\displaystyle \frac{{\rm{d}}}{{\rm{d}}r}\left(\displaystyle \frac{1}{{\rm{\Delta }}}\displaystyle \frac{{{\rm{d}}R}_{{lm}\omega }}{{\rm{d}}r}\right)-V(r){R}_{{lm}\omega }=-{{\mathcal{T}}}_{{lm}\omega }(r),\end{eqnarray} \tag{ 2 }$$

where ${{\mathcal{T}}}_{{lm}\omega }(r)$ is the source term which is connected with the energy-momentum tensor of the test particle around a black hole, and the potential is

$$\begin{eqnarray}&&V(r)=-\displaystyle \frac{{K}^{2}+4i(r-M)K}{{\rm{\Delta }}}+8{\rm{i}}\omega r+\lambda ,\end{eqnarray} \tag{ 3 }$$

where $K=({r}^{2}+{a}^{2})\omega -{ma},\quad \lambda ={E}_{{lm}}+{a}^{2}{\omega }^{2}-2{amw}-2$. The other is the angular equation

$$\begin{eqnarray}&&\displaystyle \frac{1}{\mathrm{sin}\theta }\displaystyle \frac{{\rm{d}}}{{\rm{d}}\theta }\left(\mathrm{sin}\theta \displaystyle \frac{{{\rm{d}}}_{-2}{S}_{{lm}}^{a\omega }}{{\rm{d}}\theta }\right)\\ &&+\left[{(a\omega )}^{2}{\mathrm{cos}}^{2}\theta +4a\omega \mathrm{cos}\theta -\left(\displaystyle \frac{{m}^{2}-4\;m\mathrm{cos}\theta +4}{{\mathrm{sin}}^{2}\theta }\right)+{E}_{{lm}}\right]{}_{-2}{S}_{{lm}}^{a\omega }=0,\end{eqnarray} \tag{ 4 }$$

where ${}_{-2}{S}_{{lm}}^{a\omega }(\theta )$ is the spin-weighted angular function.

If ${{\mathcal{T}}}_{{lm}\omega }(r)=0$, equation (2) becomes a homogeneous equation, and can be solved quickly and accurately by a semi-analytical numerical method developed by Fujita and Tagoshi [38, 39]. The homogeneous solutions of the radial Teukolsky equation are expressed in terms of two kinds of series of special functions. The first one consists of a series of hypergeometric functions and is convergent at the horizon

$$\begin{eqnarray}&&{R}_{{lm}\omega }^{{\rm{H}}}={{\rm{e}}}^{{\rm{i}}\epsilon \kappa x}{(-x)}^{-s-{\rm{i}}(\epsilon +\tau )/2}{(1-x)}^{{\rm{i}}(\epsilon -\tau )/2}{p}_{\mathrm{in}}(x),\end{eqnarray} \tag{ 5 }$$

where ${p}_{\mathrm{in}}$ is expanded in a series of hypergeometric functions as

$$\begin{eqnarray}&&{p}_{{\rm{in}}}(x)=\displaystyle \sum _{n=-\infty }^{\infty }{a}_{n}F(n+\nu +1-{\rm{i}}\tau ,-n-\nu -{\rm{i}}\tau ;1-s-{\rm{i}}\epsilon -{\rm{i}}\tau ;x),\end{eqnarray} \tag{ 6 }$$

and $x=\omega ({r}_{+}-r)/\epsilon \kappa$, $\epsilon =2\;M\omega$,$\kappa =\sqrt{1-{a}^{2}}$, $\tau =(\epsilon -{ma})/\kappa$. The hypergeometric function $F(\alpha ,\beta ;\gamma ;x)$ can be found in mathematics handbooks.

The second one consists of a series of Coulomb wave functions, and is convergent at infinity. The homogeneous solution of the Teukolsky equation is

$$\begin{eqnarray}&&{R}_{C}={z}^{-1-s}{(1-\epsilon \kappa /z)}^{-s-{\rm{i}}(\epsilon +\tau )/2}{f}_{\nu }(z),\end{eqnarray} \tag{ 7 }$$

where ${f}_{\nu }(z)$ is expressed in a series of Coulomb wave functions as

$$\begin{eqnarray}&&{f}_{\nu }(z)=\displaystyle \sum _{-\infty }^{\infty }{(-{\rm{i}})}^{n}\displaystyle \frac{{(\nu +1+s-{\rm{i}}\epsilon )}_{n}}{{(\nu +1-s+{\rm{i}}\epsilon )}_{n}}{a}_{n}{F}_{n+\nu }(-{\rm{i}}s-\epsilon ,z),\end{eqnarray} \tag{ 8 }$$

and $z=\omega (r-{r}_{-}),{(a)}_{n}={\rm{\Gamma }}(a+n)/{\rm{\Gamma }}(a)$, ${F}_{N}(\eta ,z)$ is a Coulomb wave function. The outgoing homogeneous solution can be expressed in Coulomb wave functions as,

$$\begin{eqnarray}&&{R}_{{lm}\omega }^{\infty }={A}_{-}^{\nu }{z}^{-1-2s}{e}^{{\rm{i}}(z+\epsilon \mathrm{ln}z)},\end{eqnarray} \tag{ 9 }$$

where the coefficient ${A}_{-}^{\nu }$ is

$$\begin{eqnarray}&&{A}_{-}^{\nu }={2}^{-1-s+{\rm{i}}\epsilon }{e}^{-\pi \epsilon /2}{e}^{-{\rm{i}}\pi (\nu +1-s)/2}\displaystyle \sum _{-\infty }^{+\infty }{(-1)}^{n}\displaystyle \frac{{(\nu +1+s-{\rm{i}}\epsilon )}_{n}}{{(\nu +1-s+{\rm{i}}\epsilon )}_{n}}{f}_{n}^{\nu }.\end{eqnarray} \tag{ 10 }$$

Note that equations (6)–(10) involve a parameter ν, the so-called renormalized angular momentum. The key part of this semi-analytical method is to search the renormalized angular momentum numerically (see [38, 39] for details).

With ${R}_{{lm}\omega }^{\infty }(r)$ and ${R}_{{lm}\omega }^{H}(r)$ in hand, and the source term ${{\mathcal{T}}}_{{lm}\omega }$ (for an equatorial-circular case, see our previous work [40]),

$$\begin{eqnarray}{{\mathcal{T}}}_{{lm}\omega }(r) & = & \displaystyle \int {\rm{d}}t{{\rm{\Delta }}}^{2}\{({A}_{{nn}0}+{A}_{n\bar{m}0}+{A}_{\bar{m}\bar{m}0})\delta (r-{r}_{0})\\ & & +{\partial }_{r}[({A}_{n\bar{m}1}+{A}_{\bar{m}\bar{m}1})\delta (r-{r}_{0})]+{\partial }_{r}^{2}[{A}_{\bar{m}\bar{m}2}\delta (r-{r}_{0})]\}.\end{eqnarray} \tag{ 11 }$$

we have

$$\begin{eqnarray}&&{Z}_{{lm}}^{H}=\displaystyle \frac{\pi }{{\rm{i}}{\omega }_{m}{B}_{{lm}\omega }^{\mathrm{in}}{T}_{r}}{{ \mathcal I }}^{H}({r}_{0}),\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{Z}_{{lm}}^{\infty }=\displaystyle \frac{\pi {B}^{\mathrm{hole}}}{{\rm{i}}{\omega }_{m}{B}_{{lm}\omega }^{\mathrm{in}}{D}_{{lm}\omega }^{\infty }{T}_{r}}{{ \mathcal I }}^{\infty }({r}_{0}),\end{eqnarray} \tag{ 13 }$$

where

$$\begin{eqnarray}{{ \mathcal I }}_{{lm}\omega }^{H,\infty } & = & \left[{R}_{{lm}\omega }^{H,\infty }({A}_{{nn}0}+{A}_{\bar{m}n0}+{A}_{\bar{m}\bar{m}0})-\displaystyle \frac{{{\rm{d}}R}_{{lm}\omega }^{H,\infty }}{{\rm{d}}t}({A}_{\bar{m}n1}+{A}_{\bar{m}\bar{m}1})\right.\\ & & +{\left.\displaystyle \frac{{{\rm{d}}}^{2}{R}_{{lm}\omega }^{H,\infty }}{{{\rm{d}}t}^{2}}{A}_{\bar{m}\bar{m}2}\right]}_{{r}_{0},\theta =\pi /2}.\end{eqnarray} \tag{ 14 }$$

All quantities (${A}_{{nn}0},{A}_{\bar{m}n0}$, etc) shown in the above equations were given explicitly in [9], and the harmonic frequency is ${\omega }_{m}=m{{\rm{\Omega }}}_{\phi }$. The calculation of coefficients ${A}_{{nn}0},{A}_{\bar{m}n0},\cdots$ involves the solution of the angular equation (4). There are several routes to calculate the spin-weighted spheroidal function ${}_{-2}{S}_{{lm}}^{a\omega }$. In this paper, we adopt the method described in [9].

Then the amplitudes ${Z}_{{lm}}^{H,\infty }$ fully determine the fluxes of gravitational radiations to infinity and the horizon,

$$\begin{eqnarray}&&{\dot{E}}^{\infty ,{\rm{H}}}=\displaystyle \sum _{{lm}}\displaystyle \frac{{| {Z}_{{lm}}^{{\rm{H}},\infty }| }^{2}}{4\pi {\omega }_{m}^{2}},\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{\dot{{L}_{z}}}^{\infty ,{\rm{H}}}=\displaystyle \sum _{{lm}}\displaystyle \frac{m{| {Z}_{{lm}}^{{\rm{H}},\infty }| }^{2}}{4\pi {\omega }_{m}^{3}},\end{eqnarray} \tag{ 16 }$$

where the overdot stands for ${\rm{d}}/{\rm{d}}t$, E and L_z mean energy and angular momentum respectively. The notations $\infty$ and ${\rm{H}}$ on E and L_z mean the fluxes to infinity and the horizon respectively. The gravitational waveform can be expressed as:

$$\begin{eqnarray}&&{h}_{+}-{{\rm{i}}h}_{\times }=\displaystyle \frac{2}{r}\displaystyle \sum _{{lm}}\displaystyle \frac{{Z}_{{lm}\omega }^{{\rm{H}}}}{{\omega }_{m}^{2}}{S}_{im}^{a{\omega }_{m}}(\theta ){e}^{-{\rm{i}}{\omega }_{m}t+{\rm{i}}m\phi }.\end{eqnarray} \tag{ 17 }$$

As we mentioned, the CPU expense for simulating the EMRI evolution by numerically solving the Teukolsky equation is quite huge. Adopting PN fluxes to do evolution is fast but will lose accuracy when the small body approaches its innermost stable circular orbit (ISCO) around the central black hole. In this section, we use the Teukolsky-based data of a few points to fit out a set of polynomials for replacing the PN fluxes or original Teukolsky-based fluxes. The idea is very simple: we select a few points between the initial radius and ISCO, and calculate the Teukolsky-based fluxes and waveforms at these points; we then use these data to fit out a set of polynomials. The fitted polynomials are used as fluxes and waveforms; they are functions of the radius r to the black hole.

We choose 7, 9, 11 and 13 points (Chebshev nodes) to fit out the 6th, 8th, 10th and 12th-order polynomials and compare the results with factorized PN ones (of course, one also can use more points to fit out these polynomials). For completeness, we need four sets of polynomials for the orbital evolution and waveform extraction. Two of them are for fluxes down to infinity and the horizon:

$$\begin{eqnarray}&&{\dot{E}}^{\infty }=\displaystyle \sum _{i=0}^{n}{a}_{i}/{r}^{i},\quad {\dot{E}}^{{\rm{H}}}=\displaystyle \sum _{i=0}^{n}{b}_{i}/{r}^{i},\end{eqnarray} \tag{ 18 }$$

where n is the order of polynomials. For waveforms, from equation (17), the term which can be fitted by polynomials is

$$\begin{eqnarray*}&&{H}_{{lm}}\equiv \displaystyle \frac{{Z}_{{lm}\omega }^{{\rm{H}}}}{{\omega }_{m}^{2}}{S}_{{lm}}^{a{\omega }_{m}}(\theta ).\end{eqnarray*}$$

Considering this term is a complex function, we need to use $2\times (n+1)$ fitting polynomials for each (l, m) mode:

$$\begin{eqnarray}&&\mathrm{Re}[{H}_{{lm}}]=\displaystyle \sum _{i=0}^{n}{R}_{{lm}}^{i}/{r}^{i},\quad \mathrm{Im}[{H}_{{lm}}]=\displaystyle \sum _{i=0}^{n}{I}_{{lm}}^{i}/{r}^{i},\end{eqnarray} \tag{ 19 }$$

where Rⁱ_lm and Iⁱ_lm are the polynomial coefficients.

For reducing interpolation errors, we use the Chebshev nodes to produce the interpolating points in our model:

$$\begin{eqnarray}&&{r}_{i}=\displaystyle \frac{{r}_{b}+{r}_{a}}{2}+\displaystyle \frac{{r}_{b}-{r}_{a}}{2}\mathrm{cos}\displaystyle \frac{2i+1}{2{n}^{\prime }+1}\Pi ,\end{eqnarray} \tag{ 20 }$$

where ${r}_{a},\;{r}_{b}$ are the boundaries of the calculation area and ${n}^{\prime }$ is the total number of nodes, and i goes from 0 to n'−1.

In figures 1 and 2, we compare the fitted polynomials of different orders with the factorized-resummation PN fluxes. We can find that the 6th-order polynomials do not have good performance, but the factorized-resummation PN fluxes which are used in the EOB model perform worst. Both the 10th and 12th-order polynomial can fit very well to the numerical Teukolsky-based fluxes. It is assumed that eLISA will observe gravitational waves of EMRIs at a typical frequency $\sim {10}^{2}\;{\rm{Hz}}$ and the total wave cycle is about $N\sim {10}^{5}$ for one yr. Thus, the relative error of energy luminosity required to establish the accuracy for the cycle ${\rm{\Delta }}N\leqslant 1$ must be $\leqslant {10}^{-5}$ in circular orbit cases [41]. Therefore, the 6th-order polynomials cannot satisfy the requirement of accuracy. The 8th-order one is at the edge of this requirement. Both the 10th and 12th-order polynomials can meet this accuracy requirement. We also observe that the 10th and 12th-order polynomials for a = 0.9 do not perform as well as the other cases, especially approaching the inner boundary (see figures 1 and 2). This is because for the cases of a approaching 1, the energy fluxes increase very steeply when the small body approaches the ISCO which is very close to the horizon. This problem should be solved if we decompose the domain and use piecewise polynomials. Nevertheless, the 10th-order polynomial fluxes are still much better than the PN ones, and can meet the requirement well.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Considering simplification and saving CPU time, in this paper we decide to use the 10th-order polynomials to fit energy fluxes and waveforms. We list the 10th-order polynomial coefficients of energy fluxes in tables 1 and 2 which are used in figures 1 and 2. The polynomial coefficients for calculating the (2, 2)-mode waveform displayed in figure 6 are listed in tables 3 and 4.

Table 1.  Polynomial parameters for infinity fluxes.

	a0	a1	a2	a3	a4	a5	a6	a7	a8	a9	a10
a = 0.9	1.51e-7	−1.22e-5	4.40e-4	−9.43e-3	1.36e-1	4.90e0	−5.27e0	−1.05e1	2.98e1	−4.34e1	3.99e1
a = 0.7	6.61e-7	−5.98e-5	2.41e-3	−5.69e-2	8.76e-1	−2.90e0	5.22e1	−2.75e2	8.81e2	−1.64e3	1.44e3
a = 0.0	1.46e-6	−1.72e-4	9.08e-3	−2.83e-1	5.76e0	−7.39e1	7.63e2	−5.02e3	2.21e4	−5.77e4	7.19e4
$a=-0.9$	1.26e-6	−1.85e-4	1.22e-2	−4.78e-1	1.22e1	−2.09e2	2.63e3	−2.21e4	1.25e5	−4.20e5	6.78e5

Table 2.  Polynomial parameters for horizon fluxes.

	b0	b1	b2	b3	b4	b5	b6	b7	b8	b9	b10
a = 0.9	−3.44e-9	2.93e-7	−1.09e-5	2.25e-4	−2.66e-3	1.34e-2	1.07e-1	−3.19e0	8.84e0	−8.33e0	2.37e0
a = 0.7	8.88e-8	−7.97e-6	3.17e-4	−7.39e-3	1.11e-1	−1.14e0	8.00e0	−3.89e1	1.17e2	−2.04e2	1.58e2
a = 0.0	4.84e-7	−5.70e-5	3.00e-3	−9.28e-2	1.88e0	−2.59e1	2.48e2	−1.62e3	6.97e3	−1.79e4	2.11e4
$a=-0.9$	5.30e-7	−7.76e-5	5.11e-3	−1.99e-1	5.08e0	−8.90e1	1.08e3	−9.10e3	5.05e4	−1.69e5	2.63e5

Table 3.  Polynomial coefficients for waveform (2, 2) mode: real part.

	${R}_{22}^{0}$	${R}_{22}^{1}$	${R}_{22}^{2}$	${R}_{22}^{3}$	${R}_{22}^{4}$	${R}_{22}^{5}$	${R}_{22}^{6}$	${R}_{22}^{7}$	${R}_{22}^{8}$	${R}_{22}^{9}$	${R}_{22}^{10}$
a = 0.9	4.00e-5	4.96e-1	−9.20e-1	−3.43e-1	−6.23e0	1.73e2	−1.02e3	3.20e3	−6.14e3	6.78e3	−3.27e3
a = 0.7	2.88e-5	4.96e-1	−9.20e-1	1.03e0	−1.98e1	2.75e2	−1.57e3	5.33e3	−1.18e4	1.54e4	−9.09e3
a = 0.0	1.78e-5	4.97e-1	−8.67e-1	5.44e0	−6.22e1	6.73e2	−4.23e3	1.89e4	−5.58e4	9.83e4	−7.12e4
$a=-0.9$	1.38e-5	4.97e-1	−7.79e-1	1.17e1	−1.17e2	1.34e3	−9.37e3	4.84e4	−1.46e5	2.09e5	1.13e5

Table 4.  Polynomial coefficients for waveform (2, 2) mode: imaginary part.

	${I}_{22}^{0}$	${I}_{22}^{1}$	${I}_{22}^{2}$	${I}_{22}^{3}$	${I}_{22}^{4}$	${I}_{22}^{5}$	${I}_{22}^{6}$	${I}_{22}^{7}$	${I}_{22}^{8}$	${I}_{22}^{9}$	${I}_{22}^{10}$
a = 0.9	−8.83e-5	1.20e-2	−1.72e0	2.85e0	4.48e1	−3.08e2	1.16e3	−2.89e3	4.82e3	−5.07e3	2.49e3
a = 0.7	−6.68e-5	1.03e-2	−1.66e0	1.52e0	5.76e1	−4.20e2	1.91e3	−6.10e3	1.36e4	−1.88e4	1.18e4
a = 0.0	−3.52e-5	7.18e-3	−1.52e-1	−2.30e0	1.02e2	−9.80e2	6.81e3	−3.42e4	1.19e5	−2.53e5	2.54e5
$a=-0.9$	−2.06e-5	5.32e-3	−1.41e0	−6.21e0	1.59e2	−2.02e3	1.86e4	−1.24e5	5.72e5	−1.62e6	2.20e6

The polynomials above are obtained from the numerical data of 11 points between the ISCO and $10\;M+{r}_{\mathrm{ISCO}}$. Actually one can choose any range which is interesting for the research to fit out the corresponding polynomials. One can also choose the number of points and the order of polynomials based on the accuracy requirement. However, the polynomials we obtained may be used directly to the neighbor area (except for the area inside the ISCO). In figure 3, we show the validity of our polynomials (fitted from the data between ${r}_{{\rm{ISCO}}}$ and ${r}_{{\rm{ISCO}}}+10\;M$) in the further area. We can see that the polynomials for ${\dot{E}}^{\infty }$ fitted from the numerical data can be used directly to the further area. At the same time, the polynomials for the energy flux down to the horizon do not perform very well in the further area. However, ${\dot{E}}^{{\rm{H}}}$ is much less than ${\dot{E}}^{\infty }$ (only 10⁻⁵ of the latter). In many cases, one can just simply omit ${\dot{E}}^{{\rm{H}}}$.

Zoom In Zoom Out Reset image size

In addition, one can also fit the energy fluxes and waveforms by the polynomials with the PN parameter $x\equiv {(v/c)}^{2}={({GM}{\rm{\Omega }}/{c}^{3})}^{2/3}$, then equations (18) and (19) are transferred to

$$\begin{eqnarray}&&{\dot{E}}^{\infty }=\displaystyle \sum _{i=0}^{n}{a}_{i}^{\prime }{x}^{i},\quad {\dot{E}}^{{\rm{H}}}=\displaystyle \sum _{i=0}^{n}{b}_{i}^{\prime }{x}^{i},\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&\mathrm{Re}[{H}_{{lm}}]=\displaystyle \sum _{i=0}^{n}R{{}^{\prime }}_{{lm}}^{i}{x}^{i},\quad \mathrm{Im}[{H}_{{lm}}]=\displaystyle \sum _{i=0}^{n}I{{}^{\prime }}_{{lm}}^{i}{x}^{i}.\end{eqnarray} \tag{ 22 }$$

The equations (18–22) are essentially PN expansions. However, all coefficients are obtained from the numerical fitting of the Teukolsky-based data to guarantee the accuracy, in contrast to the analytical expressions of the PN approximation such as the factorized-resummation ones.

During the inspiralling process, the orbit of a small body is semi-circular, and the frequency-domain Teukolsky-based waveform is highly accurate in this process. As discussed in the previous section, we use the 10th-order polynomials to replace the original numerical Teukolsky fluxes and waveforms. For producing the 10th-order polynomials, firstly we need flux and waveform data of 11 points during the evolution. These data are calculated from the Teukolsky equation. Once the evolution area is decided, the 11 interpolation points are generated by the Chebyshev nodes in this work. Calculating the Teukolsky-based fluxes and waveforms at 11 points only takes a few seconds by a desktop.

Using the data on these 11 points, we can give out the 10th-order polynomials immediately. With the flux polynomials at hand, we can use the EOB dynamics to evolve the orbits very quickly. The well-known EOB formalism was first introduced by Buonanno and Damour more than ten years ago to model comparable-mass black hole binaries [42, 43], and was also applied in small mass-ratio systems [15, 18–20, 44, 45]. The EOB dynamical evolution equations under radiation reaction for a quasi-circular orbit can be given as [46–48]

$$\begin{eqnarray}&&\dot{r}=\displaystyle \frac{\partial {H}_{\mathrm{EOB}}}{\partial {p}_{r}},\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&\dot{\phi }=\displaystyle \frac{\partial {H}_{\mathrm{EOB}}}{\partial {p}_{\phi }},\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&\dot{{p}_{r}}=-\displaystyle \frac{\partial {H}_{\mathrm{EOB}}}{\partial r}+{{ \mathcal F }}_{\phi }\displaystyle \frac{{p}_{r}}{{p}_{\phi }},\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&\dot{{p}_{\phi }}={{ \mathcal F }}_{\phi },\end{eqnarray} \tag{ 26 }$$

where ${{ \mathcal F }}_{\phi }=\dot{E}/\dot{\phi }$, and $\dot{E}$ is the energy flux of gravitational radiation. For a non-spinning test particle, the Hamiltonian is

$$\begin{eqnarray}&&{H}_{\mathrm{NS}}={\beta }^{i}{p}_{i}+\alpha \sqrt{{\mu }^{2}+{\gamma }^{{ij}}{p}_{i}{p}_{j}},\end{eqnarray} \tag{ 27 }$$

where $\mu ={m}_{1}{m}_{2}/M$, ${m}_{1},\;{m}_{2}$ are the masses of two bodies respectively, $M={m}_{1}+{m}_{2}$ and

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{1}{\sqrt{-{g}^{{tt}}}},\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{\beta }^{i}=\displaystyle \frac{{g}^{{ti}}}{{g}^{{tt}}},\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&{\gamma }^{{ij}}={g}^{{ij}}-\displaystyle \frac{{g}^{{ti}}{g}^{{tj}}}{{g}^{{tt}}},\end{eqnarray} \tag{ 30 }$$

${g}^{\mu \nu }$ is the inverse Kerr metric.

In our previous ET codes, the energy fluxes are obtained from equation (15) by calculating the Teukolsky equation numerically at every time step. Then, the fluxes are sourced to the EOB dynamical equations (23)–(26) to drive the particle inspirals into the central black hole. This method will cost a lot of CPU time for calculating the Teukolsky equation. Now we use the flux polynomials (18) or (21) as the source in the EOB dynamical equation. For getting the coefficients of (18) or (21), we also need to calculate the Teukolsky equation but only at a few points (11 points are chosen in this paper; calculation is finished in a few seconds). With the flux polynomials at hand, the EOB dynamical equations can be driven in a very fast way. At every time step, we calculate the waveforms by the polynomials (19) or (22). We list our numerical algorithm here for clarity:

• 1.  
  determine the calculating area of the EMRIs by distance (${r}_{0},{r}_{\mathrm{end}}$) or GW frequency (${f}_{0}^{\mathrm{GW}},{f}_{\mathrm{end}}^{\mathrm{GW}}$);
• 2.  
  choose ${n}^{\prime }$ points (Chebyshev nodes) in this area to calculate the Teukolsky-based energy fluxes and waveform, then to fit out n order flux and waveform polynomials based on the accuracy requirement (${n}^{\prime }\gt n$).
• 3.  
  Calculate the dynamics of the small body by solving the EOB dynamical equations (23)–(26) with the source term obtained in (2) from the initial point (r₀, ${\phi }_{0}$, ${{p}_{r}}_{0}$,  ${{p}_{\phi }}_{0}$);
• 4.  
  compute the waveforms by the polynomials obtained in (2) at every evolution temporal point until the end.
• 5.  
  If necessary, we need an iteration method for the correction of the non-circular orbit: take the evolved data ${{\rm{\Omega }}}_{\phi }$, v_r at the ${N}^{\prime }$ points instead of the original ${{\rm{\Omega }}}_{\phi }^{{\rm{circ}}},{v}_{r}^{{\rm{circ}}}=0$ into step (2) and repeat all the remaining steps.

We use the orbital frequency ${{\rm{\Omega }}}_{\phi }^{{\rm{circ}}}$ from the circular orbital condition to calculate the Teukolsky equation. However, during the evolution process, the small body has radial velocity and the orbit is not exactly circular. Therefore the practical orbital frequency calculated by (24) is different from ${{\rm{\Omega }}}_{\phi }^{{\rm{circ}}}$. This is why we need an iteration procedure in step (5). For mass ratio $\mu /M\lesssim {10}^{-3}$, the iteration is not necessary, but for the one of 10⁻², we recommend to do the iteration. Please see figure 4 for the comparison of the orbital frequency and radial velocity.

Zoom In Zoom Out Reset image size

In figure 5, as an example, we show the orbital evolution of EMRIs with $\mu /M={10}^{-3}$ for a = 0. Because the plunge process is dominated by the conservation dynamics, the plunge orbits past the ISCO are achieved here just by ignoring the radiation reaction.

Zoom In Zoom Out Reset image size

The corresponding waveforms are demonstrated in figure 6. All these evolutions are achieved in 1 second by one CPU of a desktop. In the final year of the inspiral, an EMRI waveform has 10⁵ cycles [30]. The computation time of such waveform is less than 300 seconds (single CPU) for an EMRI. However, the complexity of EMRI waveforms makes this procedure challenging. The inspiral waveform depends on 14 different parameters [31]. Based on the analysis of [30], if we assume that only about eight of these 14 parameters affect the phase evolution, it will need $300\times {({10}^{5})}^{8}\;{\rm{s}}$ to produce all waveform templates with a single CPU! Only when we assume there are two to three parameters, this computation time becomes acceptable by using thousands of CPUs with parallel technology. If we use the standard Teukolsky techniques, it will take about a few tens of days to compute one year's evolution by using our codes with one CPU! It is much longer than the new procedure developed in this paper. In this sense, our new technique makes great progress in saving CPU time, though it is still far from practical demands.

Zoom In Zoom Out Reset image size

To confirm the validation of our polynomial waveforms, we compare the fitting polynomials with the numerical Teukolsky waveform near the ISCO. We find that for the mass ratio 1:1000, the match of our model with the numerical waveforms looks quite good (see the left panel of figure 6). For confirmation, based on matched-filter technology (an optimal method when searching for known signals in noisy data), we use the spectral noise density of eLISA [32] to calculate the overlap between our polynomial waveforms and the direct Teukolsky ones. According to the basic set-up in matched filtering, the overlap ${O}_{a,b}$ between two time series of signals a and b is defined as,

$$\begin{eqnarray}&&{O}_{a,b}=\displaystyle \frac{(a| b)}{\sqrt{(a| a)(b| b)}},\end{eqnarray} \tag{ 31 }$$

where the product $(a| b)$ is given as

$$\begin{eqnarray}&&(a| b)=2{\displaystyle \int }_{0}^{\infty }\displaystyle \frac{\tilde{a}(f){\tilde{b}}^{*}(f)+{\tilde{a}}^{*}(f)\tilde{b}(f)}{{S}_{{\rm{n}}}(f)}{\rm{d}}f.\end{eqnarray} \tag{ 32 }$$

The overhead tildes stand for the Fourier transform and the star stands for a complex conjugation. The quantity S_n(f) is the spectral noise density curve taken from [32]. From these equations, the overlaps of the two kinds of waveforms (our waveform polynomials and the Teukolsky one) are $99.95\%,\;99.99\%,\;99.76\%$ and 100% for $a=0.9,\;0.7,\;0$ and −0.9 respectively. This means that our polynomial waveforms are faithful during the whole evolution process. This also confirms our previous claim: for the mass ratio $\lesssim {10}^{-3}$, we do not need an iteration to correct the quasi-circular approximation. However, as we have claimed, for a mass ratio around 1:100, we suggest that one should take step (5) to obtain more accurate waveforms because of the breaking down of adiabatic approximation.

In this paper, we use the Teukolsky-based fluxes and waveforms at a few points on the evolution route of EMRIs to fit out a set of polynomials for fluxes and waveforms. A circular orbital condition is adopted to solve the orbital parameters for the numerical calculation of the Teukolsky equation in a frequency domain.

Essentially, these flux and waveform polynomials are also a kind of PN expansion but all coefficients are obtained from the fitting of numerical data. Usually the PN coefficients are calculated from analytical expressions, such as the resummation PN waveform and Fujita's 11th PN results for EMRIs [33]. As we can see from figure 1, our 10th or 12th-order polynomials are much better than the resummation PN fluxes, especially for the spin black hole cases. Comparing them with the results shown in figure 4 of [33], we find that our results are also more accurate than the 11th PN analytical fluxes which have very long expressions.

Yunes et al used the numerical Teukolsky-based fluxes to calibrate higher-order PN coefficients and add them to resummed PN analytical expressions [27, 28]. Though their results are quite good, the accuracy of our flux polynomials is a little better than theirs. Furthermore, the fitting of the coefficients of polynomials costs less CPU time than their calibration method. Solving the Teukolsky equation to give a few sets of fluxes and waveforms by the semi-analytical method introduced in section 2 just needs a few seconds, and the fitting process almost does not add extra CPU time. However, Yunes has mentioned that they need O(10) minutes to complete calibration [28].

Just for demonstration, in the present paper we use a total of 11 points to fit out a set of 10th-order polynomials for fluxes and waveforms. However, using more points can fit out more accurate polynomials. For examples, 10th or 12th-order polynomials obtained from 20 points will give out better fluxes and waveforms. This will only add a little CPU time (twice of the case of 11 points).

With these polynomials at hand, EMRIs can be evolved in a very fast way and the waveforms can be extracted at the same time. The parameters listed in tables (1)–(3) can be used directly for the GW data analysis. In the present paper, we only fit out the one-dimensional polynomials of r. In principle, one can try to fit out two-dimensional polynomials of both r and a. Unfortunately we find that the coefficients of some polynomials vary suddenly while a just changes a little. We will leave this to future work. For many different a of Kerr black holes, principally we can build a database of flux and waveform polynomials by scanning a large number of different Kerr parameters (each a only takes a few seconds to produce the polynomials). Our fast and accurate polynomial model could be useful in the waveform-template calculations for future GW detectors such as eLISA.

This work was supported by the National Natural Science Foundation of China (No. U1431120, No.11273045).