{'fraction_non_alphanumeric': 0.02833610397461085, 'fraction_numerical': 0.01214044632512216, 'mean_word_length': 4.355138926355544, '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': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Modeling a circular equatorial test-particle in a Kerr spacetime

Extreme Mass Ratio Inspirals (EMRIs) are one of the main gravitational wave (GW) sources for a future space detector, such as eLISA/NGO, and third generation ground-based detectors, like the Einstein Telescope. These systems present an interest both in astrophysics and fundamental physics. In order to make a high precision determination of their physical parameters, we need very accurate theoretical waveform models or templates. In the case of a circular equatorial orbit, the key stumbling block to the creation of these templates is the flux function of the GW. This function can be modeled either via very expensive numerical simulations, which then make the templates unusable for GW astronomy, or via some analytic approximation method such as a post-Newtonian approximation. This approximation is known to be asymptotically divergent and is only known up to 5.5PN order for the Schwarzschild case and to 4PN order for the Kerr case. A way to improve the convergence of the flux is to use re-summation methods. In this work we extend previous results using the Pad\'e and Chebyshev approximations, first by taking into account the absorption of the GWs by the central black hole which was neglected in previous studies, and secondly by using the information from the Schwarzschild and absorption terms to create a Kerr flux up to 5.5PN order. We found that these two additions both improve the convergence. We also demonstrate that the best re-summation method for improving the flux model is based on a flux function which we call the"inverted Chebyshev approximation".

later extended to Kerr orbits [17]. In both cases, they found that in general, the Padé flux has a better and more monotonous convergence than the PN flux. Even so, the Padé solution was not complete due to the existence of poles in the approximation at certain PN orders, and also due to the fact that at some PN orders the new approximation was worse than the original PN series. In previous articles it was demonstrated that one could create both sub and super-diagnonal Padé approximations, depending on the number of terms in either the numerator or the denominator. The problem of poles could be cured by moving to a sub-diagonal approximation should a pole occur in the superdiagonal approximation and vice versa. However, due to the inherent form of the approximation, there is no cure for the cases where the approximation simply performs badly. Furthermore, in the Kerr case, a variety of spins were investigated, but the flux was kept at the 4PN order, which we now know to be divergent. The second solution which was proposed was to replace the divergent power series by a Chebyshev series, and was tested in the non-spinning test-mass and comparable mass cases [18,19]. As a Chebyshev series is closely related to the minimax polynomial, they were demonstrated to obtain a better convergence than using the PN or Padé fluxes. This improvement was due to two aspects : the first is that while power series have a radius of convergence that may not cover the entire domain of interest, Chebyshev series have an ellipse of convergence that spans the entire domain of interest. Secondly, as we shall see later, the PN factorized flux contains logarithmic terms which begin at the 3PN order. This part of the series is untreatable with Padé approximation as it yields coefficients that are either null or infinite. The Chebyshev approximation has no problems with the fact the the first non-zero term of the series is at 3PN order and we can thus define a Chebyshev series here as well.
In this work we extend these previous results, first, by taking into account the absorption of the GW by the central black hole. Secondly, in the spinning case, as the 4PN flux is divergent, we include the information that we have in the Schwarzschild case up 5.5PN in the Kerr flux function. Then, we search for the best factorized flux function, from the function proposed in Ref. [16] where the velocity at the LSO is used as a normalising velocity, and of some of its variants, which include using the velocity at the photon ring as a normalising velocity or inverting the factorized flux function. Finally, we study the convergence of this more complete PN flux, and of the Padé and Chebyshev flux models, for different orders of approximation and different spins.
We will see that the flux needs to be treated differently for different types of orbits : the non-spinning, retrograde and prograde cases. This separation stems from the differing relativistic nature involved as the test-particle approaches the central black hole and will influence our choice of normalizing velocity in the flux function. We will use the invariant velocity at the photon ring for the retrograde case, and the invariant velocity at the last stable orbit for the non-spinning and prograde cases. We will also see that in most situations, it is better to use an inverted factorized flux, introduced in Ref [17], rather than the original factorized flux, proposed by Ref [16]. And finally we will see that for the Chebyshev re-summation, we will not necessarily use the highest known PN order, especially if that order less convergent than at lower orders.
This work is organized as follow, in Section II we will describe the PN flux. In Section III, we will describe the Padé and Chebyshev approximations, and we will show how to construct these respective fluxes. In Section IV, we will present our results for the different spins. Finally, in Section V we will provide our conclusion.

II. POST-NEWTONIAN GW FLUX.
Search algorithms for GW astronomy are based on the technique of matched filtering [20], which cross-correlates the received signal with theoretical templates of known parameter values. In order to determinate as precisely as possible the physical parameters of the observed source, we need very accurate templates. As matched filtering is very sensitive to the GW phase, it is therefore imperative that we have a very accurate phase model. The GW phase can be written in the form where φ 0 is a constant reference GW phase (which can be taken to be when the signal is first seen in the detector), E(v) is the binding energy and E ′ (v) = dE(v)/dv. F (v) is the flux of the GW, and v is the local linear velocity given by v ≡ M/r, where r is the radial coordinate in Boyer-Lindquist coordinates and M is the mass of the central black hole. In the case of equatorial and circular orbits, we have an exact analytical expression for the binding energy [21] E(v, q) = η 1 − 2v 2 + qv 3 where η = µM/(µ + M ) 2 ∼ µ/M is the symmetric mass ratio, µ is the mass of the stellar mass object, and q is the dimensionless spin parameter defined as q ≡ a/M .
From this expression we can determine three particular orbits. First, we can determinate the position of the last stable orbit (LSO) which corresponds to the condition E ′ (v) = 0 [21] where The signs + and − correspond to prograde and the retrograde orbits respectively. Then we can obtain the last unstable orbit, also known as the photon ring, where the derivative of the binding energy, E ′ (v), has a pole [21,22], Finally, we can determine the position of the black hole horizon, where the energy is zero The position of these three orbits varies with the sign and magnitude of q. In the Schwarzschild case the horizon, the photon ring and the last stable orbit are respectively at 2M , 3M and 6M . When q is positive and increasing, these orbits move to smaller radius until the coincide at radius r = M in the limiting case of a maximally spinning BH, i.e. q = 1. When q is negative and decreasing, the LSO moves away from the black hole, while the photon ring and the horizon again approach the central BH.
While we have a exact expression for the binding energy in the case of circular equatorial orbit, it is not the case for the GW flux function F (v). Instead we use a PN approximation where we have a series approximation in powers of v/c. This approximation has been calculated up to 5.5PN order in the case of a Schwarzschild BH and up to 4PN in the case of a Kerr BH. The general form of the flux in both cases is given by the expression [5][6][7][8][9][10][11][12] where n = 8 and 11 in the Kerr and Schwarzschild cases respectively, and F N (x) is the Newtonian flux given by Here, x is the invariant PN velocity parameter observed at infinity, which is related to the orbital angular frequency Ω by x = (M Ω) 1/3 (Note that this is a different use of notation to the comparable mass literature where we define The relation between the invariant velocity x and local linear velocity v is given by [23] x which reduces to x = v in the Schwarzschild limit. The coefficients a k (q) and b k (q) appearing in the flux function are given in Appendix A. They are the sum of three contributions : a spinless term, a spin dependent term, and an absorption term [23]. The absorption term comes form the absorption of the GW by the central BH. Until now, in the Kerr case, the two series in Eqn (8) have been stopped at n = 8 as we do not know the spindependent term for higher orders. However, this 4PN flux function is divergent. In this work we propose an extension of this series up to n = 11 where we use the spin-independent and the absorption terms for n > 8. While this flux function is still incomplete, due to the lack of knowledge of spin dependent coefficients at greater than 4PN order, we shall demonstrate that this flux function has a better convergence when compared with a numerical flux that the truncated 4PN Kerr flux function.

III. APPROXIMATING THE GW FLUX FUNCTION.
The field of numerical analysis presents a number of ways of improving the convergence of a series. Padé approximation is a powerful method due to the fact that it is the best way of representing a power series by means of a rational polynomial, i.e.
with n + 1 = N + M + 1 and the power series of the right hand side being identical to the original power series. Theory suggests that if we think of the Padé approximation as a matrix, the best convergence is achieved by staying as close to the diagonal as possible, i.e. keep the number of terms in the numerator and denominator equal. At odd orders of approximation, it is possible to allow either the number of terms in the numerator (super-diagonal) or the denominator (sub-diagonal) increase by one term. Thus, if we choose M = N + ǫ, where ǫ = 0 or 1, we have the diagonal and sub-diagonal approximations, respectively. The advantage of the sub-diagnoal Padé approximant is that we can write it in the form of a continued fraction In this form, when we go from the order n to the order n+1, only the (n+1) th coefficient needs to be calculated [24,25], all lower order coefficients remain constant. In contrast, with a super-diagonal approximation, each higher order of approximation requires a re-calculation of all coefficients. Padé approximation is a very useful method for accelerating the convergence of a series as not only can it well represent zeros and poles, but it can provide new information regarding the underlying function. The main problem with Padé approximation is when it gets the approximation wrong. As with any power series, the Padé approximation is ultimately based on the Weierstrass theorem which only guarantees convergence with a high number of terms, it is impossible to know at any particular order if the Padé approximation will perform well. On the other hand, the best way of increasing the convergence of a power series by using a polynomial representation is to use the minimax polynomial, which unfortunately is not that easy to find. The minimax polynomial approximates the underlying function such that we are presented with an oscillating equal error curve. This approximation has the ability to allow the error to grow in parts of the domain where we already have high accuracy in order to improve the accuracy in parts of the domain where the error begins to grow. Closely related to the minimax polynomial is a polynomial based on Chebyshev polynomials of the first kind which can be defined using the iterative relation with The Chebyshev polynomials are orthogonal polynomials and are a special case of the Gegenbauer polynomials (see Ref [18,26] for a more in depth explanation). They are defined in the interval [−1, 1], where at the n th order they have n zeros and n + 1 extrema such that |T n (x)| = 1. As most problems are not limited to this interval and are defined in the interval [a, b], so we define the shifted Chebyshev polynomials of the first kind, T * (x) = T (s), by The shifted Chebyshev polynomials can be calculated using the following iterative expression with The Chebyshev approximation again tries to minimize the maximum error, albeit less successfully than the minimax polynomial, but sufficiently so to give a more convergent series at lower orders of approximation than the truncated power series. While a Padé approximation can provide new information on an underlying function, a Chebyshev approximation essentially re-packages the information that we have from the power series into a more useful form. Another beauty of the Chebyshev series is that each coefficient at higher order is smaller than the previous coefficient.
As the Chebyshev polynomials have a maximum magnitude of unity, this means that the truncation error associated with a Chebyshev series is essentially the value of the truncated coefficient. Ultimately, the performance of the Chebyshev series is dependent on the generating power series. If a series at a particular order is particularly divergent, then the corresponding Chebyshev series with perform better, but will still be divergent. We will see that in this case, starting at a lower order of approximation in the power series can result in a more convergent Chebyshev approximation.
The first method that was introduced in order to try and improve the convergence of the flux function was based on a Padé approximation [16]. Due to a number of issues with the PN approximation to the flux, a direct application of Padé approximation is not possible. The first problem with the PN flux are the logarithmic terms which appears at 3PN, which effect the convergence of the approximation. In order to reduce this effect, it was shown that one could normalize the velocity parameter appearing in the logarithm and also factor out the logarithmic terms from the original series. This led to the creation of a new flux function where and In Ref [16], the invariant velocity at the last stable orbit is used as the normalizing velocity (i.e. x N = x LSO ). This value was used as the intention was to model the inspiral of the test particle only as far as the LSO. As we approach this point the effect of the logarithmic terms are greatly diminished. However, as there are no constraints in choosing the value of the normalizing velocity, we will also investigate other possibilities. A second issue is that the linear term of the PN flux is null, which creates a problem when we use the Padé approximation. A null term in a series can result in null and infinite Padé coefficients when we write the expansion in continued fraction form. Furthermore, we know that the flux function has a pole at the photon ring, so in order that we use the full potential of the Padé approximation it was required that the function exhibit this pole [24,27]. Therefore, by multiplying the flux by 1 − x xpr we obtain the factorized flux where and where x pr is the invariant velocity at the photon ring. There exits another factorized flux introduced in Ref [17], that we will call the inverted factorized flux, which is obtained by analytically inverting the previous equation, which gives, where n k=0 and the d k coefficients are calculated using As the sum over the l k 's in Eqn (21) and Eqn (24) only starts at k = 6, we cannot apply the Padé approximation as it would generate null, undefined and infinite coefficients in the Padé approximation. So the Padé approximation is only applied to the non-logarithmic sum in the above equations, and we hence obtain the Padé flux [16] and the inverted Padé flux [17] F The Padé approximant was studied in the non-spinning [16] and spinning [17] cases, and it was shown that in both case the Padé flux better modeled a numerical flux than the PN flux. However it was also seen in both cases that the Padé approximation can suffer from singularities at different PN orders. In principle we can recover from this by moving to a super-diagonal approximant. Furthermore, in some cases the Padé approximation was a worse fit than the original PN series. This was due to the fact that the factorized flux is highly divergent at certain PN orders and performed much worse than the PN series. While the Padé approximation was able to improve the convergence of the factorized flux, it still remained less convergent than the PN flux. This also turned out to be the case in the Kerr study. The convergence of the flux was studied for different spin values of the central BH, but only at the highest Kerr order of 4PN. However, we now know that this order is more divergent than lower orders and thus the study may not be truly reflective of just how well we can model the Kerr flux function.

B. Chebyshev flux approximation.
To model the Chebyshev flux, we defined the interval between the initial and LSO velocities [x ini , x lso ], where we now further define two variables From Eqn (16), we have the new iterative relation with For our study we took x ini = 0, however we could just as easily take the velocity when the system is first seen in the detector.
Contrary to the Padé approximation, the Chebyshev re-summation is not effected by null coefficients and can be applied to both sums appearing in the factorized flux function. We therefore define the Chebyshev flux function as [18] with the Chebyshev series As the inverted Padé flux was more convergent when compared with the numerical flux, we thus similarly define an inverted Chebyshev flux function The Chebyshev re-summation was studied in the Schwarzschild case [18], and in this case it has already shown a more stable behavior which assured the minimization of the error between the approximated and numerical fluxes when the order of approximation increased. In general the Chebyshev approximation provided better accuracy at lower orders than either the PN or Padé fluxes. The only approximation at which it did not perform the best was at the 5.5PN order where it was better than the PN flux but worse than the Padé flux, which for unknown reasons is particularly good. We believe that as the Padé approximation had a highly oscillatory nature of performance, it just happened to match the numerical flux closely at this PN order. We do not believe that this was a sign of convergence of the Padé approximation as the 6PN approximation (if we had it) could again have been worse than the PN approximation. While it is not truly reflective as a metric due to the fact that a bad approximation could luckily coincide with the numerical flux, the error at the LSO was also investigated for each approximation. While the PN and Padé fluxes oscillated wildly in their performance, the Chebyshev approximation essentially reached the same constant error value at all orders of approximation which was at least a factor of two better, and sometimes a factor of 13 better, than either the PN or Padé flux models.

IV. RESULTS.
We know that depending on the value of the spin, the position of the LSO is different. For prograde orbits, the test-particle approaches closer to the central black and the system is more relativistic, whereas for retrograde orbits the small body falls into the black hole earlier, and the system is less relativistic. Based on this physical consideration, we decided to organize our study into three parts : the non-spinning case, the retrograde case and the prograde case. In each case, we found the best factorized flux for each re-summation technique, and using these best-choice fluxes we compared the two re-summed fluxes with the PN approximation at different orders and for different spins. A. The Schwarzschild Case.
Our first study concerns the non-spinning case, and before working on any resummation methods we first need to find the best factorized flux which will form the basis of these methods. We saw in Section III A that two aspects of the factorized flux can been tuned : the choice of the normalization in the logarithmic term (we tried x lso and x pr ) and the choice to use an inverted or an non-inverted factorized flux. We tried these four cases at 5.5PN and we applied on each of them the Padé and the Chebyshev re-summations, the results are shown in Fig 1. The first thing that we see is that for both Padé and Chebyshev fluxes, and for the two normalizations, the inverted flux models have a smaller error than the non-inverted flux. Depending on which normalization factor is used, the improvement ranges between a factor of two and an order of magnitude. This validates the proposition made in Ref [17] to use the inverted flux. We also see that for both the Padé and Chebyshev fluxes, using x lso for the normalization, which was the choice made in [16], produces a smaller maximum error. Consequently, the best flux, in this sense, is the inverted flux using x lso for the normalization. While the results presented here are for the 5.5PN case, we have confirmed a similar behaviour at lower orders of approximation.
Using this factorized flux we compare the PN (black solid curves), Padé (red dot-dashed curves), Chebyshev (green double-dot-dashed curves) and inverted Chebyshev (blue dashed curves) fluxes with a numerical flux. This comparison ranges from 2PN to 5.5PN and is shown in Fig 2. The first thing that we see is that, except at the 2PN order where the PN flux has a fortunate coincidence with the numerical flux and the orders where the Padé flux has a singularity, the three re-summations are always more convergent than the PN flux. If we compare the performance of each flux model as we approach the last stable orbit (and again keeping in mind that this is not on its own the best metric of convergence) we can see that the PN approximation returns a relative error ǫ = |1 − F A /F N | (where F A is the approximated flux, and F N is the numerical flux) of 0.03 ≤ ǫ ≤ 0.4 at various orders of approximation. The Padé approximation, when not suffering from a singularity, has relative errors of between 3 × 10 −2 ≤ ǫ ≤ 8 × 10 −2 , except at the 5.5PN order where it has an error of 1.7 × 10 −3 . The Chebyshev approximation, as expected, tries to minimize the maximum error and produce an oscillating error curve. As was shown in Ref [18], the Chebyshev approximation converges to an almost constant error value of 1.8×10 −2 at the LSO. The inverted Chebyshev approximation performs fantastically in this case achieving a constant error value of approximately ǫ = 10 −3 at all orders of approximation. Therefore, we can see from the figure that the inverted Chebyshev flux at 2PN order is as accurate as the 5.5PN Padé approximation (that is difficult to explain) and at least an order of magnitude better than all other models at each order of approximation.
At wider separations, the Chebyshev approximation allows the error to float a little, but in this domain, we do not require the same level of accuracy as we do closer to the central BH.

B. Retrograde Orbits.
As in the Schwarzschild case, we start by searching for the best factorized flux for each re-summation method. We demonstrate in the bottom two cells of Fig 3, the comparison between the inverted/non-inverted fluxes, using different normalization constants, at 5.5PN for a spin q = −0.5. The result from this case is representative of the results for all other retrograde spins. As in the Schwarzschild case, we see that once again for both the Padé and Chebyshev flux models, the inverted flux is the more convergent model. One could argue in the Padé case that the standard flux looks to be the best choice as it has a smaller error as we approach the LSO, but it is not the case due to the presence of a singularity in the flux. This pole is also present at all other spins and other differing orders of approximation, which is why we reject using this flux. We also can see that using x pr is a marginally better choice than using x lso as a normalizing velocity. However, the results are so close that one could continue to use the LSO value as a normalizing factor. Now that we have a model for the optimal factorized flux, we now compare the PN flux (black solid curves), the inverted Padé (red dot-dashed curves) and the inverted Chebyshev (blue dashed curves), from 2PN to 5.5PN for the four retrograde spin values of q = −0.25, −0.50, −0.75, −0.95, which are shown in Fig 4-7. The first thing that we notice is that the addition of the Schwarzschild and absorption terms to the PN flux to 5.5PN order provides a more convergent PN flux model than the standard PN Kerr model at 4PN. This result has the knock-on effect that the ensuing 5.5PN factorized flux is better performing than the corresponding factorized flux at 4PN order, thus ensuring that the re-summed fluxes starting at 5.5PN order are more convergent than previously studied.
If we first investigate the performance of the PN flux, we see that this approximation has a slightly improved convergence at 5.5PN order as we move from a spin of -0.25 to a spin of -0.95. This is due to the fact that the higher retrograde spins being less relativistic and the test-particle plunging further out from the central BH. In general, at q = −0.25 the error at the LSO is between 3 × 10 −2 ≤ ǫ ≤ 0.3, improving to between 10 −2 ≤ ǫ ≤ 0.2 at a spin of q = −0.95. While the performance of the PN flux illustrates the oscillatory convergence that we would expect, we should point out that as we move towards a maximally counter-rotating spin, the 4.5PN order PN approximation performs better and better with an errors of 3 × 10 −2 , 10 −2 , 2 × 10 −4 and 3 × 10 −3 at increasing orders of spin. This mimics the unexplainable behaviour we observed in the Schwarzschild case for the Padé approximation at 5.5PN order, and in this case the PN flux just happens to have a fortunate coincidence with the numerical flux as we approach the LSO.
The inverted Padé flux does not perform particularly well in the case of retrograde motion. At all spin values, when not suffering from singularities, there are levels of approximation where the Padé approximation is inferior to the original PN flux function. This is due to the fact that the factorized flux is divergent at these orders and performs much worse than the original PN flux function. In the cases where it does provide an improvement, the overall gain is a factor approximately 1.5 over the PN approximation with errors of between 2 × 10 −2 ≤ ǫ ≤ 0.2 and 10 −2 ≤ ǫ ≤ 0.1 at spins of -0.25 and -0.95 respectively. The inverted Chebyshev approximation performs consistently well at all spin values, with relative errors of ǫ = 7 × 10 −3 at q = −0.25 and ǫ = 10 −2 at a spin of q = −0.95. Regardless of the spin value, the inverted Chebyshev flux is approximately an order of magnitude better at low orders of approximation and at least a factor of two better at higher orders. More importantly is the fact that this approximation is almost always better than both the PN and inverted Padé fluxes, and essentially converges to the same order of magnitude error at the LSO for all retrograde spin values. In the cases where it is less convergent than the other approximations, we again highlight the fact at times the PN approximation has a fortunate coincidence with the numerical flux as we approach the LSO, and the Padé approximation at 5PN order seems to be very good for high retrograde spin values. However, neither of these facts can be considered as global trends, so it is clear that for retrograde orbits, the inverted Chebyshev flux provides a more convergent flux model than either of the other two methods. Furthermore, the inverted Chebyshev flux function at 2.5 or 3PN order appears to have a sufficiently small error at all values of separation to be sufficient as a basis for a template.

C. Prograde Orbits.
We finish our study by looking at the prograde case. As previously, we first search for the best factorized flux. In the top two panels of Fig 3 we plot the comparison between the Padé and Chebyshev fluxes at 5.5PN order for q = 0.5. We see that in both cases the normalization using x lso is marginally better than a normalization using x pr . In the Padé case the best choice of flux is the standard Padé flux due to the presence of the poles in both of the inverted flux models. On further investigation, these poles appear for all values of prograde spin, so we choose the non-inverted factorized flux as the basis for the Padé approximation. On the other hand, for the Chebyshev flux, there is not much difference between the four flux models until we approach the LSO. In this case that standard Chebyshev flux, with either normalization, begins to diverge. We will therefore use the inverted Chebyshev flux with x lso as a normalizing velocity as our preferred model. As for the retrograde case we compare the PN (black solid curves), Padé (red dot-dashed curves) and the inverted Chebyshev fluxes (blue dashed curves), from 2PN to 5.5PN and for q = 0.25, 0.50, 0.75, 0.95, in Fig 8-11. First of all, we can see that, contrary to the retrograde case, when the spin increases the convergence of the flux decreases. This is explained by the fact that when the spin increases, the test-particle is now orbiting much closer to the central black hole and the system becomes much more relativistic. In this case the PN flux breaks down faster. Before discussing the results in full, we should mention here that during our investigation, we noticed that the inverted Chebyshev flux was not performing as well as we had expected. When delving deeper into this issue, we discovered that the factorized flux function at 5.5PN order behaved much worse than the original PN approximation. The Chebyshev approximation, while being able to improve on this somewhat, was unable to recover from the high divergence of the factorized flux. We then studied the behaviour of the factorized flux at different PN orders and discovered that the factorized flux at 5PN order was better performing than the 5.5PN PN flux. As the performance of the Chebyshev approximation is highly dependent on the starting information, we decided to construct the inverted Chebyshev flux using the factorized flux at 5PN order. In the following figures, a new curve represented this truncated inverted Chebyshev flux is represented (green double-dot-dashed curves). The case of prograde orbits is sufficiently difficult to describe that we divide it into two regimes. We noticed that even up to a spin of q = 0.5 the orbits are only slightly more relativistic than the Schwarzschild case, while from spins of q = 0.75 onwards, the orbits become truly relativistic to the point that the PN flux begins to break down sooner and sooner. However, one global observation is that, even in the prograde case, the PN flux at 4.5 to 5.5PN order is more convergent than the standard Kerr 4PN flux. This again justifies the inclusion of the Schwarzschild and absorption terms at higher PN orders.
Firstly, we will discuss the results at spins of 0.25 and 0.5. In both these cases, the PN flux performs as expected with fractional errors of between 8 × 10 −3 ≤ ǫ ≤ 0.5 and 2 × 10 −2 ≤ ǫ ≤ 0.9 respectively. In the lower spin case the error of ǫ = 8 × 10 −3 is due again to a lucky coincidence with the numerical flux as we approach the LSO. In general the mean error is on the order of ǫ = 0.1. In the q = 0.5 case, things begin to become a little clearer and we start to oscillate wildly between highly convergent and highly divergent. While the Padé flux presents a large improvement on the PN flux, in terms of stability, we still see some of the wild fluctuations associated with the PN flux, with relative errors of between 4 × 10 −3 ≤ ǫ ≤ 5 × 10 −2 and 3 × 10 −4 ≤ ǫ ≤ 6 × 10 −2 . Here as in the PN case, the extremely convergent values are due coincidences with the numerical fluxes as we approach the LSO. In general the errors are of the order of ǫ = 5 × 10 −2 and ǫ = 3 × 10 −2 respectively for both spin values. At this level of spin, the Padé flux suffers from singularities at the 3.5PN order which rule out using these orders for template creation. The inverted Chebyshev flux starting from the divergent factorized flux at 5.5PN order is almost always better than the PN flux at a spin of q = 0.25 achieving an almost constant error of ǫ = 2 × 10 −2 at the LSO. For this low spin, it in general also outperforms the Padé approximation. However for a spin of q = 0.5, while still outperforming the PN flux function, the 5.5PN inverted Chebyshev flux is constantly worse than the Padé approximation, achieving a maximum spin error of only ǫ = 0.1. This is also reflected at higher spin values and demonstrates the increasing divergence of the 5.5PN factorized flux. However, if we now concentrate on the 5PN inverted Chebyshev flux, we can see that it approaches an almost constant error value of ǫ = 10 −2 for a spin of q = 0.25 and ǫ = 5 × 10 −3 for a spin of q = 0.5. For both spin values, it is consistently better than all other approximation models. While in one or two cases it is beaten by a fortunate PN or Padé flux model, it clearly provides a more stable model, even at higher values of r/M . We next move onto the more relativistic spins of q = 0.75 and q = 0.95. For these two spin values, there is not much point in discussing the relative errors for the various approximations as none of them can cope with the high relativistic nature of the orbits. In this case, a more indicative metric may be the separation at which the relative error reaches ǫ = 0.1. With this in mind, the PN flux is the least performing approximation achieving the threshold error at separations of between 5 ≤ r/M ≤ 10 for both spins. In each case the PN flux at 5.5PN performs exceptionally well, getting as close as r = 3.5M before having an error of 10 −1 . However, this is clearly not reflective of convergence, as a 6PN flux could perform much worse. The Padé flux again improves over the PN flux, except at the 5.5PN order, reaching the threshold error at separations of between 3 ≤ r/M ≤ 5 in both cases, with the odd exception performance in the highest spin case at 5PN where a separation of r = 2.5M is reached. However, the Padé flux approximation now begins to suffer from singularities not only at 3.5PN order as we saw for lower spins, but also now at the 3PN order. The 5.5PN inverted Chebyshev flux has a relative error of 10 −1 at separations of around r = 5M in both cases. While it achieves a consistent separation before reaching the threshold, it can only be marginally considered to be better than the PN flux, and is constantly outperformed by the Padé flux. However, the 5PN inverted Chebyshev flux reaches the threshold value at separations of approximately r = 4M in both spin cases. While sometimes outperformed at a particular approximation by either the PN or Padé flux models, we believe that the 5PN inverted Chebyshev flux represents the most stable and consistently performing model for prograde orbits.

V. CONCLUSION.
In this work we extended the studies of [16][17][18] by taking into account the absorption of the GW by the central black hole, and by adding the information from the non-spinning case, we created a 5.5PN Kerr flux. We saw that, because the 4PN flux is divergent, adding these higher order terms improves the convergence of the PN flux in all cases.
By first looking at the non-spinning case, we searched for the most convergent factorized flux at 5.5PN for each re-summation method. We tried a number of different combinations : the standard factorized flux given in Ref. [16] or the inverted factorized flux proposed in Ref. [17], and we also tried x lso , the invariant speed at the last stable orbit and x pr , the invariant speed at the photon ring, as a normalization factor in the logarithmic term of the flux.  It appeared that for both Padé and Chebyshev flux models, the best choice was to use the inverted flux normalized using x lso , which confirms the choices made in [16] and [17], respectively. Then we compared the two resummation methods and the PN flux, at different orders, with a numerical flux. We observed that, except when it presents a pole, the Padé approximation is in general better than the PN flux. However, while the standard Chebyshev flux performed extremely well, the better approximation was clearly the inverted Chebyshev flux function.
Next we studied the case of retrograde orbits, and once again we started by searching for the best factorized flux.  Contrary to the previous case, for negative spins the marginally better normalization choice is using x pr , and as in the Schwarzschild case, the inverted flux is the best choice for the two techniques. When comparing the different fluxes we saw that, except for a few particular cases, the re-summation methods are better than the PN flux. At times however, the PN flux has a fortunate crossing with the numerical flux near the LSO and therefore appears to work well. However, the fluctuations is convergence from one order to another means that it is difficult to rely on the PN approximation. In general, the Padé approximation improves on the PN flux, but again there are orders of approximation where it performs worse than the PN approximation. The inverted Chebyshev approximation outperforms both other flux models, improving the accuracy of the flux model by a factor of two. Our final study involved prograde orbits and once again we found that the best normalization is using x lso , but while the inverted flux was the best choice for the Chebyshev approximation, it appeared that the standard flux is a better choice for the Padé flux. Here we found that the PN approximation performs worse and worse as we increase the spin magnitude, except at the 5.5PN order where it performs surprisingly well. The Padé flux is a great improvement over the PN flux, but again oscillates between good and bad performance. In this case, we found that the inverted Chebyshev flux at 5.5PN order did not perform as well as expected. Upon further investigation, we found that this was due to the factorized flux being divergent at this order. However, by using the 5PN factorized flux, we were once again able to construct a very convergent Chebyshev flux model. While more marginal, the inverted Chebyshev flux is the most stable model, converging to the same order of error at the LSO.
This work improves upon previous studies and introduces a new inverted Chebyshev flux model. This model is highly convergent and stable. In most cases, this flux model is more convergent at the 2.5PN level than either the PN or Padé fluxes at 5.5PN order. We believe that this model provides a good basis for the creation of theoretical templates for GW astronomy.
where ψ (0) (z) and ψ (1) (z) are respectively the digamma and the trigamma functions. For ease and speed of calculation, the digamma and trigamma functions were generated using a Chebyshev-Padé approximation such that the maximum error between the approximation and the analytic expressions is less than 10 −10 for all spin values.