{'fraction_non_alphanumeric': 0.037144796818931236, 'fraction_numerical': 0.016244787120550868, 'mean_word_length': 4.452934955050238, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fast Evolution and Waveform Generator for Extreme-Mass-Ratio Inspirals in Equatorial-Circular Orbits

In this paper we discuss the development of a fast and accurate waveform model for the quasi-circular orbital evolution of extreme-mass-ratio-inspirals (EMRIs). This model simply employs the data of a few numerical Teukoulsky-based energy fluxes and waveforms to fit out a set of polynomials for the entire fluxes and waveforms. These obtained polynomials are accurate enough in the entire evolution domain, and much more accurate than the resummation post-Newtonian (PN) energy fluxes and waveforms, especially when the spin of a black hole becomes large. The dynamical equation we adopted for orbital revolution is the effective-one-body (EOB) formalism. Because of the simplified expressions, the efficiency of calculating the orbital evolution with our polynomials is also better than the traditional method which uses the resummed PN analytical fluxes. Our model should be useful in calculation of waveform templates of EMRIs for the gravitational wave detectors such as the evolved Laser Interferometer Space Antenna (eLISA).


I. INTRODUCTION
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 EMRI is about or less than O(10 1 ) 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 time scale. 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 eLISA (evolved Laser Interferometer Space Antenna), a space gravitational wave observatory supported by the European Space Agency now [1]. A pathfinder has been launched in 2015 to pave the way for 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 the post-Newtonian (PN) approximation, the second one is numerical relativity, and the last one is black hole perturbation theory. As analyzed in a great deal of literatures (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 a good performances for spinning black holes (see for example Ref. [4] and references insides); Numerical relativity still can not handle the binary black hole systems with extreme mass ratio. An recent paper documents the simulation of a case of mass ratio 1:100 (without spin) with only two orbits' evolution [5].
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 calculation of 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 challenged about how to greatly reduce the CPU time of the simulation of EMRI waveform, 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 Schwarzschild black hole, Lackeos and Burko fitted polynomials mixed PN expressions to do 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 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 expansions, but all the coefficients of such "PN expansions" are obtained numerically from the fitting of the Teukolsky-based data.
In the next section, we introduce our EOB-Teukolsky codes (ET codes) shortly. The details of our fitting polynomial method are presented in the 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 .

II. EOB-TEUKOLSKY FREQUENCY-DOMAIN CODES
Our ET codes include two main parts: one is the EOB dynamics driver, 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 Teukolskybased 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 the section 4.
The Teukolsky equation in ET codes can be solved in frequency-domain for inspiralling phase and in time-domain for plunge and merge states. In the present work, we focus on only 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][37][38][39] to replace the previous numerical integration method [40].
After decomposing the Weyl curvature (complex) scalar ψ 4 in a Fourier series [25], where ρ = −1/(r − ia cos θ), the Teukolsky equation is divided into two parts. One is the radial master equation where T lmω (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 where K = (r 2 + a 2 )ω − ma, λ = E lm + a 2 ω 2 − 2amw − 2. The other is the angular equation where −2 S aω lm (θ) is the spin-weighted angular function. If T lmω (r) = 0, Eq. (2) becomes a homogeneous equation, then can be solved quickly and accurately by a semi-analytical numerical method developed by Fujita and Tagoshi [38,39].
The homogeneous solutions of Teukolsky radial equation are expressed in terms of two kinds of series of special functions. The first one consists of series of hypergeometric functions and is convergent at the horizon where p in is expanded in a series of hypergeometric functions as and x = ω(r + − r)/ κ, = 2M ω,κ = √ 1 − a 2 ,τ = ( − ma)/κ. The hypergeometric function F (α, β; γ; x) can be found in mathematic handbooks.
The second one consists of series of Coulomb wave functions which is convergent at infinity. The homogeneous solution of Teukolsky equation is where f ν (z) is expressed in a series of Coulamb wave functions as and z = ω(r −r − ), (a) n = Γ(a+n)/Γ(a), F N (η, z) is a Coulomb wave function. The outgoing homogeneous solution can be expressed in Coulomb wave functions as, where the coefficient Note that Eqs.(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 the R ∞ lmω (r) and R H lmω (r) in hand, take into the source term T lmω (for an equatorialcircular case, see our previous work [40]), we have where All the quantities (A nn0 , Am n0 , etc.) shown in the above equations were given explicitly in [9], and the harmonic frequency is ω m = mΩ φ . The calculation of coefficients A nn0 , Am n0 , · · · involves the solution of the angular equation (4). There are several routes to calculate the spin weighted spheroidal function −2 S aω lm . In this paper, we adopt the method described in [9].
Then the amplitudes Z H,∞ lm fully determine the fluxes of gravitational radiations to infinity and horizon,Ė where the overdot stands for d/dt, E and L z mean energy and angular momentum respectively. The notations ∞ and H on E and L z mean the fluxes to infinity and the horizon respectively. The gravitational waveform can be expressed as: III. FITTING POLYNOMIAL FROM TEUKOLSKY-BASED ENERGY FLUXES

AND WAVEFORMS
As we mentioned, the CPU expense for simulating the EMRI evolution by numerical solving the Teukolsky equation is quite huge. Adopting PN fluxes to do evolution is fast but will lost the accuracy when the small body approaching it's innermost stable circular orbit (ISCO) around the central black hole. In this section, we use the Teukolsky-based flux 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, the calculate the Teukolsky-based fluxes and waveforms at these points, and 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 set of polynomials for the orbital evolution and waveform extraction. Two of them are for fluxes down to the infinity and horizon: where n is the order of polynomials. For waveforms, from Eq. (17), the term can be fitted by polynomials is Considering this term is a complex function, we need to use 2 × (n + 1) fitting polynomials for each (l, m) mode: where R i lm and I i lm are the polynomial coefficients. For reducing interpolation errors, we use the Chebshev nodes to produce the interpolating points in our model: where r a , r b are the boundaries of the calculation area and n is the total number of nodes, and i goes from 0 to n − 1.
In Fig. 1 and Fig. 2, we compare the fitted polynomials of different orders with the factorized-resummation PN fluxes. We can find that the 6th polynomials do not give a good performance. And the factorized-resummation PN fluxes which are used in the EOB model perform worst. Both the 10th and 12th order polynomial can give very good fit to the numerical Teukolsky-based fluxes. It is assumed that LISA will observe gravitational waves of EMRIs at the typical frequency ∼ 10 2 Hz and the total wave cycle is about N ∼ 10 5 for 1 yr. Thus, the relative error of energy luminosity required to establish the accuracy for the cycle ∆N ≤ 1 must be ≤ 10 −5 in circular orbit cases [41]. Therefore, the 6th polynomials cannot satisfy the requirement of accuracy. The 8th one is at the edge of this requirement.
The 10th polynomials can meet the requirement well. Considering the simplification and saving CPU time, we decide to use the 10th polynomial to fit energy fluxes and waveform in this paper. We list the 10th polynomial coefficients of energy fluxes in Tab. I and Tab. II  which are used in Fig. 1 and Fig. 2. The polynomial coefficients for calculating (2, 2)-mode waveform displayed in Fig. 6 are listed in Tab. III and IV.

IV. ORBITAL EVOLUTION AND WAVEFORM
During the inspiralling process, the orbit of small body is semi-circular, and the frequencydomain Teukolsky based waveform is highly accurate in this process. As discussed in the last section, we use the 10th polynomials to replace the original numerical Teukolsky fluxes and waveforms. For producing the 10th polynomials, firstly we need flux and waveform data of the 11 points during the evolution. These data are calculated by 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 numerically only takes 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 fastly. 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,47] where F φ =Ė/φ, andĖ is the energy-flux of gravitational radiation. For a non-spinning test particle, the Hamiltonian is where µ = m 1 m 2 /M , m 1 , m 2 are the masses of two bodies respectively, M = m 1 + m 2 and g µν is the inverse Kerr metric.
In our previous ET codes, the energy fluxes are obtained from Eq. (15) (19) or (22). We list our numerical algorithm here for clarity: (1) determine the calculating area of the EMRIs by distance (r 0 , r end ) or GW frequency (2) choose n 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 > n).
(3) calculate the dynamics of the small body by solving the EOB dynamical Eqs. (23)- (26) with the source-term obtained in (2) from the initial point (r 0 , φ 0 , p r 0 , p φ 0 ); (4) compute the waveforms by the polynomials obtain in (2) at every evolution temporal point until to the end.
(5) if necessary, we need an iteration method for the correction of non-circular orbit: take the evolved data Ω φ , v r at the N points instead of the original Ω circ φ , v circ r = 0 into the step (2) and repeat all the remainning steps.
We use the orbital frequency Ω circ φ from 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 Ω circ φ . This is why we need an iteration procedure in the step (5). For massratio µ/M 10 −3 , the iteration is not necessary, but for the one of 10 −2 , we recommend to do the iteration. Please see Fig. 4 for the comparison of the orbital frequency and radial velocity.
In Fig. 5, as an example, we show the orbital evolution of EMRIs with µ/M = 10 −3 for a = 0. Because the plunge process is dominated by the conservation dynamics, the plunge orbits passed the ISCO are archive here just by ignoring the radiation-reaction.
The corresponding waveforms are demonstrated in Fig. 6. All these evolutions are archive in 1 second by one CPU of a desktop. In the final year of the inspiral, an EMRI waveform has 10 5 circles [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 × (10 5 ) 8 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 years 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.
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 the matched-filter technology (an optimal method when searching for known signals in noisy data), we use the spectral noise density of LISA [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, where the product (a|b) is given as The overhead tildes stand for the Fourier transform and the star stands for a complex conjugation. The quantity Sn(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%, 97.01%, 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 confirm our previous claim: for the mass-ratio 10 −3 , we need not an iteration to correct the quasi-circular approximation. The results of comparison are plotted in Fig. ??. However, as we have claimed, for mass-ratio around 1:100, we suggest that one should take the step (5) to obtain the more accurate waveforms because of the large radial velocity.

V. CONCLUSIONS AND DISCUSSIONS
In this paper, we use the Teukolsky-based fluxes and waveforms at a few points on the evolution route of the EMRI to fit out a set of polynomials for fluxes and waveforms. A circular orbital condition is adopted to solve the orbital parameters for 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 Fujitas 11th PN results for EMRIs [33]. As we can see from Fig. 1, our 10th or 12th-order polynomials are much better than the resummation PN fluxes especially for the spin black hole cases. Comparing the results shown in Fig. 4 of [33], we can 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 wave-forms 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 (I-III) 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, Taiji [48] and Tianqin [49].

ACKNOWLEDGMENTS
This work was supported by the NSFC (No. U1431120, No.11273045).