Spinning test body orbiting around a Kerr black hole: Eccentric equatorial orbits and their asymptotic gravitational-wave fluxes
SSpinning test body orbiting around a Kerr black hole: eccentric equatorial orbits andtheir asymptotic gravitational-wave fluxes
Viktor Skoup´y , and Georgios Lukes-Gerakopoulos Astronomical Institute of the Czech Academy of Sciences,Boˇcn´ı II 1401/1a, CZ-141 00 Prague, Czech Republic and Institute of Theoretical Physics, Faculty of Mathematics and Physics,Charles University in Prague, 18000 Prague, Czech Republic
We use the frequency and time domain Teukolsky formalism to calculate gravitational wave fluxesfrom a spinning body on a bound eccentric equatorial orbit around a Kerr black hole. The spinningbody is represented as a point particle following the pole-dipole approximation of the Mathisson-Papapetrou-Dixon equations. Reformulating these equations we are not only able to find the tra-jectory of a spinning particle in terms of its constants of motion, but also to provide a method tocalculate the azimuthal and the radial frequency of this trajectory. Using these orbital quantities, weintroduce the machinery to calculate through the frequency domain Teukolsky formalism the energyand the angular momentum fluxes at infinity, and at the horizon, along with the gravitational strainat infinity. We crosscheck the results obtained from the frequency domain approach with the resultsobtained from a time domain Teukolsky equation solver called Teukode.
I. INTRODUCTION
An Extreme Mass Ratio Inspiral (EMRI) is one ofthe most promising events expected to be detected withfuture space-based gravitational wave (GW) detectorslike Laser Interferometer Space Antenna (LISA) [1]. AnEMRI occurs when a stellar mass compact object such asa black hole (BH) or a neutron star (secondary object)is trapped in the vicinity of a supermassive black hole(SMBH) (primary object). Due to gravitational radia-tion reaction the secondary is slowly spiralling into theprimary while emitting GWs. From these GWs it is pos-sible to extract information about the EMRI system suchas the masses of the objects, their spins etc. On a morefundamental physics level, EMRIs detection are expectedto allow us to probe the strong gravity regime around aSMBH [2].Currently in order to extract information from a GWsignal, when it is detected by the terrestrial observato-ries, it has to be uncovered from a dominating noise back-ground. To achieve this, matched filtering is employed,i.e. waveform templates for a wide range of parametersare matched with the detected time series. It is expectedthat we will have to use matched filtering for GW signalreceived by LISA as well, but not to uncover the sig-nal from the noise; in LISA’s case we will use them todisentangle overlapping GW signals from simultaneouslydetected sources. Because of this, accurate models of theGW waveform templates are planned to be produced fora wide range of parameters.To model GWs from an EMRI, first the trajectory ofthe secondary object must be reproduced. The standardway to do this is to apply the two timescale approxi-mation [3]. In an EMRI the mass ratio q ≡ µ/M liesbetween 10 − and 10 − , where µ is the secondary massand M is the primary mass. The energy changes at rate˙ E/E = O ( q ) which is very small. The timescale of theinspiral is, thus, large of the order O ( q − ), i.e q − times larger than the orbital timescale. This allows us to breakour analysis in two timescales, the fast orbital and theslow adiabatic dissipation in the constants of motion. Inthe fast one, the trajectory of the secondary over one or-bital period is close to a trajectory calculated without adissipation. The secondary is actually drifting betweenorbits characterized by a set of constants of motion. Inthis setup, the azimuthal coordinate of the inspiral canbe expanded as φ = q − φ (0) ( qt ) + φ (1) ( qt ) + O ( q ). Thefirst term of the expansion is of adiabatic order and in-cludes the contribution from the time-averaged dissipa-tive part of the first-order self-force. The second term,which is of the order of radians is called post-adiabatic and contains contributions from the conservative part ofthe first-order self-force, oscillating part of the dissipativepart of the first order self-force as well as time-averageddissipative part of the second-order self-force. The spinof the secondary contributes to the post-adiabatic termas is of the order of O ( q ) [4, 5]. In particular, for the spinmagnitude S of a secondary compact object, like a KerrBH or a neutron star, holds that S (cid:46) µ , hence the di-mensionless spin parameter defined as σ ≡ S/ ( µM ) ≤ q is of the same order as the mass ratio [6]. The phase φ is approximately proportional to the phase of the GW.Hence, to accurately model the GW fluxes, all the afore-mentioned terms must be taken into account.In this work, we deal with the contribution of the sec-ondary spin to the post-adiabatic term, in the case ofbounded equatorial orbits around a Kerr BH. The de-scription of a spinning test body moving on a curvedbackground was for the first time studied in [7–9]. Inparticular, Mathisson [10] managed to write the stress-energy tensor of an extended test body as a sum of mul-tipolar moments. When the body is sufficiently smalland compact, then it is sufficient to take into accountonly the mass (monopole) and the spin (dipole) leadingto what is known as the pole-dipole approximation , whichessentially reduces the body to a spinning test particle.Later on Papapetrou [11, 12] was able to employ the con- a r X i v : . [ g r- q c ] F e b servation law of the stress energy tensor ∇ µ T µν = 0 toderive the equations of motion for a spinning particle. Fi-nally, these equations were rewritten by Tulczyjew [13],Dixon [14–17] and Wald [18] bringing them to their mod-ern form. MPD equations have been studied in severalworks, see e.g. [6, 19–22]. Particularly, these equationssimplify when the particle is confined into the equatorialplane of the Kerr spacetime [23]. In this case, the motioncan be determined by the following constants of motion:the energy E , the component of the total angular mo-mentum parallel to the axis of the central BH J z , themass of the secondary µ and the magnitude of its spin S .In the present work, we rederive the equations of mo-tion for a spinning particle in the equatorial plane ina reduced form. This allows us to find analytical for-mulas for the constants of motion dependence on theeccentricity and the semi-latus rectum and to providea method to numerically calculate the fundamental fre-quencies. These results are then used to calculate theGW fluxes. To achieve this, we employ the Teukolskyformalism and solve the GWs perturbatively. Namely,we solve the Teukolsky equation (TE) both in the fre-quency and in the time domain with a spinning-particleas a source. In the frequency domain, the formulas pro-viding the energy and the angular momentum fluxes toinfinity and to the horizon from a spinning particle fol-lowing equatorial trajectories are novel. While, for thecalculations in the time domain, we introduce a new ap-proach to simulate the spinning source making the com-putations more efficient. Due to the GW flux balance lawin an EMRI, these fluxes equal to the rate of change ofthe constants of motion of an inspiraling spinning particle[5, 24]. Hence, once these fluxes are obtained, then theadiabatic term with the spinning-particle contribution tothe post-adiabatic term can be reconstructed.This paper is organized as follows. Sec. II briefs thedynamics of a spinning particle moving in a curved space-time. After covering the basics, the equations of mo-tion of a spinning particle are rederived in a reducedform appropriate for eccentric equatorial orbits in a KerrBH background. Subsequently, the constants of motionand the frequencies are calculated. Sec. III reviews theTeukolsky formalism calculating the GW fluxes both inthe frequency and the time domain. Finally, the fre-quency domain results are compared with the time do-main results. To make the main text more readable, wehave concentrated in a list all the dimensionless quanti-ties we use in Appendix A, Appendix B provides all theexplicit formulas for the frequency domain fluxes, whilein Appendix C our frequency domain results for a non-spinning object are compared with the ones of [25]. Fi-nally, Appendix D provides tables from the frequency do-main calculations aiming to serve as reference for futureworks.Throughout this paper, we use geometrized unitswhere the speed of light and the gravitational constantare c = G = 1. The Riemann tensor is defined as R µνκλ = Γ µνλ,κ − Γ µνλ,κ + Γ µρκ Γ ρνλ − Γ µρλ Γ ρνκ where the comma denotes partial derivative U µ,ν = ∂ ν U µ . A co-variant derivative is denoted by a semicolon U µ ; ν = ∇ ν U µ and D U µ / d τ = U µ ; ν d x ν / d τ . The signature of the metricis ( − , + , + , +). Symmetrization of indices is denoted byround brackets Φ ( µν ) = (Φ µν + Φ νµ ) /
2. For some quan-tities we prefer to use their dimensionless counterparts.They are denoted by a hat, e.g. energy ˆ E = E/µ , radialcoordinate ˆ r = r/M etc (see Appendix A). II. A POLE-DIPOLE PARTICLE MOVING ONTHE EQUATORIAL PLANE OF A KERR BLACKHOLE
The motion of a spinning test object in a curved back-ground is governed by the Mathisson-Papapetrou-Dixon(MPD) equations [9, 11, 14] which readD P µ d τ = − R µνρσ v ν S ρσ , (1)D S µν d τ = P µ v ν − P ν v µ , where P µ is the four-momentum of the particle, R µνρσ isthe Riemann tensor of the background spacetime, v µ =d x µ / d τ is the four-velocity, S µν is the spin tensor of theparticle and τ is the proper time.The stress-energy tensor T µν for a spinning particlewith its trajectory parametrized by the coordinate time t reads [26] T µν = 1 √− g (cid:18) P ( µ v ν ) v t δ − ∇ α (cid:18) S α ( µ v ν ) v t δ (cid:19)(cid:19) , (2)where for Boyer-Lindquist (BL) coordinates δ = δ ( r − r p ( t )) δ ( θ − θ p ( t )) δ ( φ − φ p ( t )) is the delta function locatedat the particle position ( r p ( t ) , θ p ( t ) , φ p ( t )) parametrizedby coordinate time. Note that by using the conserva-tion law T µν ; ν = 0, it is possible to retrieve the MPDequations.Actually, the MPD system of equations is underdeter-mined. The physical implication of the latter fact is thatthe center of the mass of the spinning object is not de-fined. To close the system of equations and to definethe centre of the mass, a spin supplementary condition(SSC) in the form S µν V µ has to be specified, where V µ is atime-like vector field. In this work, we use the Tulczyjew-Dixon (TD) SSC [13, 15] S µν P µ = 0 . (3)Under the TD SSC, the rest mass of the particle withrespect to the four-momentun µ = − P µ P µ (4)and the magnitude of the spin S = 12 S µν S µν (5)are conserved quantities (see e.g. [19]). The conservationof the above quantities is independent of the spacetimebackground. The symmetries of the spacetime introducefor each Killing vector ξ µ a specific quantity C = ξ µ P µ − ξ µ ; ν S µν , (6)which is conserved upon the evolution of the MPD equa-tions.Instead of the spin tensor, it is sometimes more conve-nient to use the spin four-vector S µ = − (cid:15) µνρσ u ν S ρσ , (7)where (cid:15) µνρσ is the Levi-Civita tensor and u ν := P ν /µ isthe specific four-momentum. The inverse relation of thisequation reads S ρσ = − (cid:15) ρσγδ S γ u δ . (8)After substituting Eq. (8) into Eq. (5), we can derivethe relation for the spin magnitude in terms of the spinfour-vector S = S µ S µ . (9)The spin four-vector is from the definition (7) orthogonalto the four-momentum P µ S µ = 0, while from Eq. (8) onesees it is orthogonal also to the spin tensor S µν S µ = 0.Finally, from Eq. (10) it can be shown that it is orthog-onal to the four-velocity v µ S µ = 0 as well.Since the MPD equtions do not provide an evolutionequation for the four-velocity, it is convenient that for theTD SSC exists an explicit relation of the four-velocity interms of the four-momentum and the spin tensor [27].This relation reads v µ = m µ (cid:18) u µ + 2 S µν R νρκλ u ρ S κλ µ + R αβγδ S αβ S γδ (cid:19) , (10)where m = − P µ v µ is the rest mass with respect to thefour-velocity v µ . This mass m is not conserved underthe TD SSC, but it is used to conserve the normalization v µ v µ = − m = A µ (cid:112) A µ − B S , (11)where A = 4 µ + R αβγδ S αβ S γδ , (12) B = 4 h κη R κιλµ P ι S λµ R ηνωπ P ν S ωπ , (13) h κη = 1 S S κρ S ηρ . (14) A. The Kerr spacetime background
Since our work deals with the motion of a spinning inthe Kerr spacetime, let us briefly introduce this space-time. The Kerr geometry in BL coordinates ( t, r, θ, φ ) is described by the metricd s = g tt d t + 2 g tφ d t d φ + g φφ d φ + g rr d r + g θθ d θ , (15)where the metric coefficients are g tt = − (cid:18) − M r Σ (cid:19) ,g tφ = − aM r sin θ Σ ,g φφ = ( (cid:36) − a ∆ sin θ ) sin θ Σ , (16) g rr = Σ∆ ,g θθ = Σwith Σ = r + a cos θ , ∆ = (cid:36) − M r ,(cid:36) = r + a . (17)The Kerr spacetime is stationary and axisymmetric.This provides two Killing vector fields, the time-like one ξ µ ( t ) and the space-like one ξ µ ( φ ) . Due to these Killingvector fields, Eq. (6) provides two constants of motion.In particular, thanks to the time-like field, the energy E = − P t + 12 g tµ,ν S µν (18)is conserved, and thanks to the space-like field, the com-ponent of the total angular momentum parallel to therotational axis of Kerr ( z axis) J z = P φ − g φµ,ν S µν (19)is conserved. These two conserved quantities can be usedto parametrize the spinning particles orbits as discussedin Sec. II C. B. Equatorial Orbits
We are interested in equatorial orbits, where θ = π/ v θ component of the four-velocity must be always zero. Theorthogonality of the spin four-vector and the four-velocity v µ S µ = 0 implies that in order to achieve v θ = 0 forarbitrary equatorial orbit all the components of the spinfour vector should be zero except from S θ , i.e. S µ = S θ δ θµ . (20)The spin is, therefore, parallel to the z axis. Fromthe orthogonality of the spin four-vector and the four-momentum P µ S µ = 0, it holds that P θ = 0.From Eqs. (9) and (20) it can be shown that S θ = −√ g θθ S where the sign is chosen such that the spin mag-nitude is positive (negative) when the spin is parallel (an-tiparallel) to the z axis. Then, from Eq. (8) the onlynonzero components of the spin tensor are S tr = − S rt = − S u φ (cid:114) − g θθ g = − S u φ r ,S tφ = − S φt = S u r (cid:114) − g θθ g = S u r r ,S rφ = − S φr = − S u t (cid:114) − g θθ g = − S u t r , (21)where g is determinant of the metric. For Kerr spacetimeon equatorial plane, it holds (cid:112) − g θθ /g = 1 /r .Let us recheck the setup for equatorial orbits in a Kerrbackground. The total derivative with respect to propertime of the θ component of four-momentum can be ex-pressed from Eq. (1)d P θ d τ = − R θνρσ v ν S ρσ − Γ θνρ P ν v ρ . (22)The rhs of this equation is equal to zero on the equatorialplane. Furthermore, Eq. (10) reduces on the equatorialplane to v θ = ( m /µ ) P θ . This implies that when v θ = 0then P θ remains zero as well. Thus, the particle stays onthe equatorial plane by just demanding that v θ = 0.From Eqs. (18), (19) and (21), P t and P φ can be ex-pressed as functions of E and J z . These expressions indimensionless quantities read u t = − ˆ E − σ ˆ r (ˆ a ˆ E − ˆ J z )1 − σ ˆ r ,u φ = M ˆ J z − σ ˆ r (cid:104)(cid:0) − ˆ a + ˆ r (cid:1) ˆ E + ˆ a ˆ J z (cid:105) − σ ˆ r . (23)When we restrict the motion to the equatorial plane,it is possible to reproduce the equations of motion for thespinning particle from Eqs. (10) and (4). In particular,we can express u r from the normalization (4) as functionof E and J z and thanks to the fact that it holds2 S rν R νρκλ u ρ S κλ = 12 µ ˆ∆ σ x ˆ r Σ σ u r (24)we can write the equations of motion asΣ σ Λ σ dˆ t dˆ τ = m µ V t (ˆ r ) , (25a)Σ σ Λ σ dˆ r dˆ τ = m µ V r (ˆ r ) = ± m µ (cid:112) R σ (ˆ r ) , (25b)Σ σ Λ σ d φ dˆ τ = m µ V φ (ˆ r ) , (25c) where Σ σ = ˆ r (cid:18) − σ ˆ r (cid:19) , (25d)Λ σ = 1 − σ ˆ rx Σ σ , (25e) V t = ˆ a (cid:18) σ ˆ r Σ σ (cid:19) x + (cid:36) ∆ P σ , (25f) R σ = P σ − ˆ∆ (cid:18) Σ σ ˆ r + x (cid:19) , (25g) V φ = (cid:18) σ ˆ r Σ σ (cid:19) x + ˆ a ˆ∆ P σ , (25h) P σ = Σ σ ˆ E − (cid:16) ˆ a + σ ˆ r (cid:17) x , (25i) x = ˆ J z − (ˆ a + σ ) ˆ E . (25j)The rest mass with respect to v µ can be expressedfrom (11) as m µ = Λ σ (cid:115) − σ ˆ r − σ − (2 − Λ σ ) σ ˆ r . (26)This expression is identical to Eq. (49) in [29]. Eqs. (25)are identical to the equations (2.19)–(2.21) in [23] up tothe parametrization with d˜ τ / d τ = m /µ where ˜ τ is theparametrization used in [23]. By dividing Eqs. (25b) and(25c) we obtain Eq. (19) in [30]. Hence, we have checkedthe validity of the above equations.To simplify the equations of motion, it is useful toreparametrize Eqs. (25) with a time parameter λ whichis similar to the Mino time [24]. Eqs. (25) and (26) implythat the relation between ˆ τ and λ isdˆ τ d λ = ˆ r (cid:115)(cid:18) − σ ˆ r (cid:19) (cid:18) − σ − (2 − Λ σ ) σ ˆ r (cid:19) . (27)Then it holds dˆ x µ / d λ = V µ where ˆ x µ = (ˆ t, ˆ r, θ, φ ) with V θ = 0. V µ can be interpreted as dimensionless four-velocity with respect to λ . C. Constants of motion as orbital parameters
Let us see how we can use the constants of motion
E, J z to parametrize bounded equatorial orbits. Forbound orbits it holds that ˆ E < r and ˆ r , i.e. R σ (ˆ r ) = 0 , R σ (ˆ r ) = 0 . (28)These points are the turning points of an equatorial ec-centric orbit. If ˆ r ≤ ˆ r , then ˆ r is the pericenter and ˆ r the apocenter. For ˆ r < ˆ r < ˆ r , Eq. (25b) implies that R σ (ˆ r ) >
0. Moreover, it holds that R (cid:48) σ (ˆ r ) ≥ , R (cid:48) σ (ˆ r ) < . (29)We can parametrize each eccentric equatorial orbit byits semi-latus rectum p and its eccentricity e , which relateto the turning points as followsˆ r = p e , ˆ r = p − e . (30)The inverse relations read p = 2ˆ r ˆ r ˆ r + ˆ r , e = ˆ r − ˆ r ˆ r + ˆ r . (31)Eq. (28) can be written as two quadratic equationsin terms of ˆ E and ˆ J z . Using the same method as inAppendix B of [31] we can rearrange the formulas (28)for energy and angular momentum to arrive at f i ˆ E − g i ˆ E ˆ J z − h i ˆ J z − d i = 0 i = 1 , f = f (ˆ r ), f = f (ˆ r ) etc. and f (ˆ r ) =ˆ a (ˆ r + 2)ˆ r + ˆ r ++ σ (cid:18) ˆ a σ ˆ r + 2ˆ a (ˆ a + σ )ˆ r + 6ˆ a ˆ r − (ˆ r − rσ (cid:19) g (ˆ r ) =2ˆ a ˆ r + σ (cid:18) ˆ aσ ˆ r + ˆ a (2ˆ a + σ )ˆ r − (ˆ r − r (cid:19) (33) h (ˆ r ) = ˆ∆ − (cid:16) ˆ a + σ ˆ r (cid:17) d (ˆ r ) = ˆ∆(ˆ r − σ ) ˆ r These functions for σ = 0 are identical to the functions(B.6) – (B.9) in [31] with z − = 0. By manipulatingEq. (32) properly, we arrive atˆ E = κρ + 2 (cid:15) ˜ σ ± (cid:112) ˜ σ (˜ σ(cid:15) + ρ(cid:15)κ − ηκ ) ρ + 4 η ˜ σ , (34)ˆ J z = ρ ˆ E − κ σ ˆ E , (35)where κ = d h − d h ,(cid:15) = d g − d g ,ρ = f h − f h , (36) η = f g − f g , ˜ σ = g h − g h are the determinants appearing in [31]. Thanks to theidentity (cid:15)ρ − κη = ˜ σζ , where ζ = d f − d f , (37)we can rearrange Eq. (34) asˆ E = κρ + 2 (cid:15) ˜ σ − J z ) ˜ σ (cid:112) (cid:15) + κζρ + 4 η ˜ σ . (38) Since for ˆ a = σ = 0 the determinant ˜ σ = 0 and theEq. (35) is singular, it is better to substitute ˆ E into Eq.(35) and rearrange it as followsˆ J z = (cid:15)ρ − κη − sgn ( ˆ J z ) ρ (cid:112) (cid:15) + κζ ( ρ + 4 η ˜ σ ) ˆ E . (39)The signs of ˆ J z appearing in Eqs. (38) and (39) have beennumerically verified for spin values | σ | ≤ E and J z for given p and e have two solutions corresponding to the corotating or-bit and the counterrotating orbit. We can choose thecoordinates such that the z axis is parallel to the totalangular momentum, i.e. ˆ J z >
0. This convention impliesthat ˆ a > a < a σ > a σ < e = 0, both the numerator and the denominator ofEq. (38) become zero. This inconvenience can be avoidedby noticing that a coefficient e can be factored out fromthe determinants (36) and canceled out in Eq. (38). Inthis fashion, the solution (38) is valid even for e = 0.Actually, this allows us to verify that for e = 0 Eqs. (38)and (39) are identical to Eqs. (59) and (60) given in [29].There is a limit between the bounded and unboundedequatorial orbits defined by a separatrix. The term un-bounded orbits includes orbits escaping to infinity andorbits plunging to the central black hole. In the case theseparatrix splits plunging and bounded orbits, it holdsthat R (cid:48) σ (ˆ r ) = 0 and R (cid:48) σ (ˆ r ) <
0. The orbit with R (cid:48) σ (ˆ r ) = 0 is an unstable circular orbit, while a tra-jectory originating from ˆ r with energy and angular mo-mentum satisfying Eqs. (38) and (39) will asymptoticallyapproach the circular orbit at ˆ r either evolved forwardor backward in time. For a given Kerr parameter ˆ a andspin σ the effective potential R σ depends on ˆ E ( p, e ) andˆ J z ( p, e ), therefore the separatrices can be plotted on the p − e plane splitting it into two parts. In one part of theplane lie the bounded orbits, while in the other part lieunbounded orbits or initial conditions, which do not cor-respond to an orbit (Fig. 1). We can see that for given e the semi-latus rectum p of the separatrix decreases withincreasing spin.Fig. 2 shows two cases of a separatrix on the ˆ J z − ˆ E plane along with a grid of constant p and e lines. Notethat the intersection point between the separatrix andthe line e = 0 lying at the left lower corner of both panelsof Fig. 2 represents ISCO. In the limiting case that ˆ r = ˆ r the orbit is circular ( e = 0)and marginally stable, since it holds that R σ (ˆ r ) = R (cid:48) σ (ˆ r ) = R (cid:48)(cid:48) σ (ˆ r ) = 0. This orbit is often called the innermost stable cir-cular orbit (ISCO). p e a - p e a σ - σ - σ σ σ p e a FIG. 1. Separatrices for different Kerr parameters and spins.Points ( p, e ) on the depicted lines correspond to orbits asymp-totically approaching the unstable circular orbit lying atˆ r = p/ (1 + e ). For given e the semi-latus rectum p of theseparatrix decreases with increasing spin. Therefore, for aspinning particle it is possible to approach the horizon closerthan a non-spinning particle. All plots are for ˆ J z > D. Frequencies of eccentric equatorial orbits.
The radial motion of a particle in the equatorial planeparametrized by the time parameter λ has a period Λ r .This period can be defined as the time needed to go fromthe apocenter to the pericenter and back. Hence, Λ r canbe found by integrating the inversion of Eq. (25b), i.e.d λ dˆ r = 1 (cid:112) R σ (ˆ r ) , (40) J z E J z E FIG. 2. Separatrice (black thick solid) in the ˆ J z − ˆ E planealong with lines of constant semi-latus rectum (grey solid)and eccentricity (grey dashed) for Kerr parameter ˆ a = − . a = 0 . σ = 0 .
5. The eccentricity lines start at e = 0 for lower energies and reach e = 1 when ˆ E = 1 withstep 0 .
1. The semi-latus rectum ranges from p = 10 to p = 20for ˆ a = − . p = 3 to p = 20 for ˆ a = 0 . J z . over the above two branches (first from ˆ r to ˆ r and thenfrom ˆ r to ˆ r ) with respect to the radius ˆ r . However, theintegration over one branch is equal to the integrationover the other. Hence, we can find the Λ r by integratingEq. (25b) over the first branch to obtain the time elapsedduring the first branch and multiply the result by two[32], i.e. Λ r = 2 (cid:90) ˆ r ˆ r dˆ r (cid:112) R σ (ˆ r ) . (41)The radial frequency can be defined as Υ r = 2 π/ Λ r . Ifwe set the initial radius to r ( λ = 0) = r , then the radius r ( λ ) is an even function and can be written as r ( λ ) = r (0) + ∞ (cid:88) n =1 r ( n ) cos( n Υ r λ ) . (42)After substituting Eq. (42) to Eqs. (25a) and (25c)and integrating them, we obtainˆ t ( λ ) = Γ λ + ∆ˆ t ( λ ) ,φ ( λ ) = Υ φ λ + ∆ φ ( λ ) , (43)where Γ and Υ are frequencies with respect to λ andfunctions ∆ˆ t ( λ ) and ∆ φ ( λ ) can be written as series ofsines.The average rate of change of the azimuthal coordinateand time with respect to λ isΥ φ = 2Λ r (cid:90) ˆ r ˆ r V φ (ˆ r ) (cid:112) R σ (ˆ r ) dˆ r , (44)Γ = 2Λ r (cid:90) ˆ r ˆ r V t (ˆ r ) (cid:112) R σ (ˆ r ) dˆ r . (45)These integrals can be solved in terms of Lauricella’shypergeometric functions [33]. However, for achievingthis, the exact values of the roots of the radial poten-tial ˆ r R σ (ˆ r ), which is eighth order polynomial in ˆ r , mustbe found. This task can be only performed numerically.Thus, instead the integrals (41), (44) and (45) were calcu-lated directly numerically. These integrals have singularpoints at ˆ r and ˆ r , but this difficulty can be overcome.Namely, first we factor out the roots R s (ˆ r ) = (ˆ r − ˆ r )(ˆ r − ˆ r ) Q (ˆ r ) , (46)where Q (ˆ r ) is a polynomial with positive and negativepowers of ˆ r . To remove the singularities, an angle likecoordinate χ ∈ [0 , π ) is used by applying the transforma-tion ˆ r = p e cos χ . (47)Then, the integrals take the formΛ r = 2 √ − e p (cid:90) π (cid:112) J ( χ ) d χ , (48)Υ φ = 2 √ − e Λ r p (cid:90) π V φ (cid:18) p e cos χ (cid:19) (cid:112) J ( χ ) d χ , (49)Γ = 2 √ − e Λ r p (cid:90) π V t (cid:18) p e cos χ (cid:19) (cid:112) J ( χ ) d χ , (50)where J ( χ ) = (cid:88) k =0 (1 + e cos χ ) k k (cid:88) l =0 j ( p ) l j ( e ) k − l (1 − e ) k − l p l (51)is a polynomial in cos χ with coefficients j ( p )0 = 1 − ˆ E ,j ( p )1 = − ,j ( p )2 = ˆ a + 2ˆ a ˆ Ex + x ,j ( p )3 = − − ˆ E ) σ − ˆ Eσx + x ) ,j ( p )4 = 4 σ ,j ( p )5 = − aσ (ˆ aσ + x ( ˆ Eσ + x )) ,j ( p )6 = σ ((1 − ˆ E ) σ − x )((1 + ˆ E ) σ + x )) and j ( e )0 = 1 ,j ( e )1 = 2 ,j ( e )2 = e + 3 ,j ( e )3 = 4( e + 1) ,j ( e )4 = e + 10 e + 5 ,j ( e )5 = 2( e + 3)(3 e + 1) ,j ( e )6 = e + 21 e + 35 e + 7 . The polynomial J ( χ ) for σ = 0 is identical to the poly-nomial (40) in [31] with Carter constant Q = 0 up to thefactor 1 − e due to a different definition of J ( χ ) used in[31].We can define the frequencies with respect to the co-ordinate time asˆΩ r = Υ r Γ = πp √ − e (cid:82) π V t (ˆ r ( χ )) / (cid:112) J ( χ )d χ , (52)ˆΩ φ = Υ φ Γ = (cid:82) π V φ (ˆ r ( χ )) / (cid:112) J ( χ )d χ (cid:82) π V t (ˆ r ( χ )) / (cid:112) J ( χ )d χ . (53)We have numerically verified the above frequency formu-las by comparing them with frequencies obtained by adirect integration of the MPD equations for the respec-tive eccentric orbits. To integrate the MPD equationsan implicit Gauss-Runge-Kutta integrator was used asdescribed in [34].The equatorial plane equations of motion (25) given in t , r and φ can be rewritten in λ , t and φ parametrizedby χ , i.e. d λ d χ = (cid:115) − e p J ( χ ) (54)dˆ t d χ = V t (cid:18) p e cos χ (cid:19) (cid:115) − e p J ( χ ) (55)d φ d χ = V φ (cid:18) p e cos χ (cid:19) (cid:115) − e p J ( χ ) (56)These equations will be used later on, when the energyand angular momentum fluxes are calculated. III. GRAVITATIONAL WAVE FLUXESA. Teukolsky formalism
To calculate the GW fluxes we employ the Teukolskyformalism. The GWs are described perturbatively usingthe Weyl curvature scalarΨ = − C αβγδ n α m β n γ m δ , (57)where n µ and m µ are components of the Kinnersleytetrad n µ = 12Σ (cid:0) (cid:36) , − ∆ , , a (cid:1) , (58) m µ = ρ √ ia sin θ, , − , i csc θ ) , (59)where ρ = − ( r − ia cos θ ) − . The Weyl scalar Ψ isgoverned by the TE s O s ψ ( t, r, θ, φ ) = 4 π Σ T (60)with spin weight s = − − ψ = ρ − Ψ in the case ofthe GWs [35].
1. Frequency domain approach
This partial differential equation can be separated intoordinary differential equations after a Fourier transformin t and φ − ψ = ∞ (cid:88) l,m π (cid:90) ∞−∞ d ωe − iωt ψ lmω ( r ) − S aωlm ( θ, φ ) , (61)where − S aωlm ( θ, φ ) is spin weighted spheroidal harmonicfunction with spin weight − (cid:90) dΩ | − S aωlm ( θ, φ ) | = 1 . (62)For simplicity we use the notation S aωlm ( θ ) = − S aωlm ( θ, D ψ lmω ( r ) = T lmω (63)is obtained for the radial part ψ lmω ( r ), where D is a dif-ferential operator that can be found e.g. in [35] and T lmω is a source term discussed below. The asymptotic behav-ior of the homogeneous solutions R lmω ( r ) of Eq. (63) isdiscussed in [4, 25]. To satisfy physical boundary con-ditions, the solution must be purely outgoing at infinityand purely ingoing at the horizon; in other words, we aredealing with a retarded solution. We will denote a ho-mogeneous solution satisfying the first condition as R + lmω and a solution satisfying the second condition as R − lmω . An inhomogeneous solution satisfying boundary condi-tions can be found using the Green function formalismas ψ lmω ( r ) = C + lmω ( r ) R + lmω ( r ) + C − lmω ( r ) R − lmω ( r ) , (64) These functions are often denoted R ∞ lmω and R H lmω or R Up lmω and R In lmω . where the amplitudes are C ± lmω ( r ) = 1 W (cid:90) ∞ r + Θ ± ( r, r (cid:48) ) R ∓ lmω ( r (cid:48) ) T lmω ( r (cid:48) )∆ ( r (cid:48) ) d r (cid:48) (65)with the invariant Wronskian W = R + lmω ( r ) ∂ r R − lmω ( r ) − ( ∂ r R + lmω ( r )) R − lmω ( r )∆( r ) (66)and the Heaviside step functions defined asΘ + ( r, r (cid:48) ) = Θ( r (cid:48) − r ) , Θ − ( r, r (cid:48) ) = Θ( r − r (cid:48) ) . (67)Since we are interested in GW fluxes at the horizon and atinfinity, we will denote the relevant amplitudes as C − lmω ≡ C − lmω ( r → r + ) and C + lmω ≡ C + lmω ( r → ∞ ) respectively.In fact, the amplitudes are constant for r < r and r > r .Te source term in (63) can be written as T lmω = (cid:90) d t d θ d φ ∆ ( T nn + T nm + T mm ) e iωt − imφ , (68)where T nn = f (0) nn ( r, θ ) √− gT nn , T nm = ∂ r ( f (1) nm ( r, θ ) √− gT nm )+ f (0) nm ( r, θ ) √− gT nm , (69) T mm = ∂ rr ( f (2) mm ( r, θ ) √− gT mm )+ ∂ r ( f (1) mm ( r, θ ) √− gT mm ) + f (0) mm ( r, θ ) √− gT mm . The functions f ( i ) ab ( r, θ ) can be found in [29]. Projectionsof the stress energy tensor onto a tetrad e ( a ) µ read T ab = 1 √− g (cid:0) C ab − C σab (cid:1) δ − √− g ∂ ρ (cid:16) ( v t ) − S ρ ( µ v ν ) δ (cid:17) e ( a ) µ e ( b ) ν , (70)where C ab = ( v t ) − P ( µ v ν ) e ( a ) µ e ( b ) ν ,C σab = ( v t ) − S ρ ( µ Γ ν ) ρλ v λ e ( a ) µ e ( b ) ν . (71)The four-vectors P µ and v µ as well as the spin tensor S µν are functions of time, the Christoffel symbols areevaluated at the coordinates of the particle r p ( t ) , θ p ( t ),the delta functions are functions of both the space coor-dinates r, θ, φ and the coordinate time t and the squareroot of the determinant √− g , the functions f ( i ) ab and thetetrad legs e ( a ) µ are functions of r and θ . In our case, e ( a ) µ , e ( b ) µ are the Kinnersley tetrad components n µ and m µ .After integrating Eq. (68) over θ and φ and Eq. (65)over r using rules for integrating delta function, we obtaina relation for the amplitudes C ± lmω = (cid:90) ∞−∞ d te iωt − imφ p ( t ) I ± lmω ( r p ( t ) , θ p ( t )) (72)where I ± lmω ( r, θ ) = 1 W (cid:18) A − ( A + B ) dd r +( A + B ) d d r − B d d r (cid:19) R ∓ lmω ( r ) . (73)The coefficients A i in their general form can be found inAppendix B.Up to this point the derivation of GW fluxes holds for ageneric orbit of a spinning particle. In the following part,we confine it to equatorial orbits with the spin parallelto the z axis as described in Sec. II B.Thanks to the fact that the quantity I ± lmω ( r p ( t ) , π/ e im (Ω φ t − φ p ( t )) is periodic in timewith frequency Ω r (see eg. [37] for details), we can writethe amplitude as a sum over discrete frequencies C ± lmω = ∞ (cid:88) n = −∞ C ± lmn δ ( ω − ω mn ) , (74) ω mn = m Ω φ + n Ω r . (75)The partial amplitudes can be calculated as Fourier co-efficients by integrating over one period T r = 2 π/ Ω r C ± lmn = Ω r (cid:90) T r d tI ± lmω mn ( r p ( t ) , π/ × exp( iω mn t − imφ p ( t )) . (76)However, it is more convenient to integrate over the timeparameter λC ± lmn = Ω r (cid:90) Λ r d λ d t d λ I ± lmω mn ( r p ( t ( λ )) , π/ × exp( iω mn t ( λ ) − imφ p ( λ )) . (77)The integration over the two branches of the motion(from r to r which correspond to λ from 0 to Λ r / r to r which correspond to λ from Λ r / r )differs only by the sign of the radial velocity. Therefore,we can break the integral to two integrals, the first from 0to Λ r / r to Λ r / ω mn t ( λ ) − mφ ( λ ) = n Υ r λ + ω mn ∆ t − m ∆ φ . (78)From the fact that ∆ t and ∆ φ are series of sines it holds∆ t (Λ r − λ ) = − ∆ t ( λ ) and ∆ φ (Λ r − λ ) = − ∆ φ ( λ ). Afterchanging the integration variable to χ , we can write theintegral as a sum over the sign D r = ± of the radialvelocity, on which the coefficients A i depend, i.e. C ± lmn = Ω r (cid:90) π d χ d λ d χ (cid:88) D r = ± d t d λ I ± lmω mn ( r ( χ ) , π/ , D r ) × exp( iD r ( ω mn t ( χ ) − mφ ( χ ))) , (79) where d λ/ d χ comes from Eq. (54), I ± lmω mn comesfrom Eq. (73) and t ( χ ), φ ( χ ) are calculated fromEqs. (55), (56).The metric perturbation h µν = O ( q ) which can be de-fined as g exact µν = g µν + h µν + O ( q ), can be calculatedfrom the Weyl scalar Ψ [38]. GWs consist of two po-larizations and the metric perturbation can be writtenas h µν = h + e + µν + h × e × µν where e + µν and e × µν are the po-larization tensors. At infinity, the relation between thestrain h = h + − ih × and the Weyl scalar isΨ ( r → ∞ ) = ¨ h/ , (80)where the dots denote derivative with respect to the BLcoordinate time t . From Eqs. (61), (64), (74) and theasymptotic behavior of R + lmω it holds h = − r (cid:88) lmn C + lmn ω mn S aω mn lm ( θ ) e − iω mn ( t − r ∗ )+ imφ , (81)where r ∗ is tortoise coordinate defined as d r ∗ / d r = (cid:36) / ∆. The stress-energy tensor of the GW can be re-constructed from the strain which yields the energy andangular momentum fluxes at infinity (cid:28) d E ∞ d t (cid:29) = ∞ (cid:88) l =2 l (cid:88) m = − l ∞ (cid:88) n = −∞ (cid:12)(cid:12) C + lmn (cid:12)(cid:12) πω mn , (82) (cid:28) d J ∞ z d t (cid:29) = ∞ (cid:88) l =2 l (cid:88) m = − l ∞ (cid:88) n = −∞ m (cid:12)(cid:12) C + lmn (cid:12)(cid:12) πω mn (83)where the brackets denote time averaging. In the equa-torial case, the average can be calculated over one period T r . Similar derivation can be made for the fluxes at thehorizon [39] (cid:28) d E H d t (cid:29) = ∞ (cid:88) l =2 l (cid:88) m = − l ∞ (cid:88) n = −∞ α lmn (cid:12)(cid:12) C − lmn (cid:12)(cid:12) πω mn , (84) (cid:28) d J H z d t (cid:29) = ∞ (cid:88) l =2 l (cid:88) m = − l ∞ (cid:88) n = −∞ α lmn m (cid:12)(cid:12) C − lmn (cid:12)(cid:12) πω mn , (85)where α lmn = 256(2 M r + ) P ( P + 3 (cid:15) )( P + 16 (cid:15) ) ω mn | C lmω mn | (86)with (cid:15) = √ M − a / (4 M r + ), P = ω mn − ma/ (2 M r + )and the Teukolsky-Starobinsky constant is0 | C lmω | = (cid:16) ( λ lmω + 2) + 4 aω ( m − aω ) (cid:17) (cid:0) λ lmω + 36 aω ( m − aω ) (cid:1) − (2 λ lmω + 3) (48 aω ( m − aω )) + 144 ω (cid:0) M − a (cid:1) . (87)The partial amplitudes C ± lmn are proportional to thesecondary mass µ and therefore, if we use dimensionlessquantities on the rhs, we obtain (cid:28) d E ∞ d t (cid:29) = q (cid:88) l,m,n (cid:12)(cid:12)(cid:12) ˆ C + lmn (cid:12)(cid:12)(cid:12) π ˆ ω mn ≡ q (cid:88) l,m,n F E ∞ lmn , (88) (cid:28) d J ∞ z d t (cid:29) = M q (cid:88) l,m,n m (cid:12)(cid:12)(cid:12) ˆ C + lmn (cid:12)(cid:12)(cid:12) π ˆ ω mn ≡ M q (cid:88) l,m,n F J z ∞ lmn , (89)where we have defined the dimensionless fluxes F E ∞ lmn and F J z ∞ lmn that do not depend on the mass ratio q . The hori-zon fluxes F E H lmn and F J z H lmn can be defined in a similarfashion. We can write the dimensionless energy and an-gular momentum loss as (cid:42) d ˆ E ∞ dˆ t (cid:43) = q (cid:88) l,m,n F E ∞ lmn , (90) (cid:42) d ˆ J ∞ z dˆ t (cid:43) = q (cid:88) l,m,n F J z ∞ lmn . (91)These fluxes can be used for calculating the evolutionof the orbital parameters p and e during an adiabaticapproximation of an inspiral.
2. Time domain approach
To verify the frequency domain calculations, we nu-merically solved the TE (60) in the time domain. Forthis, we have employed the time domain solver
Teukode which is described in [40–42].
Teukode uses the methodof lines, i.e. finite differences in space and Runge-Kuttafor evolution in time. Instead of using Kinnersley tetradand BL coordinates, it solves TE using Campanelli tetrad[43] and hyperboloidal horizon-penetrating (HH) coor-dinates ( τ, ρ, θ, ϕ ) (for their definition see Eq. (10) in[41]). These coordinates reach future null infinity I + (“scri”) and horizon at finite radial coordinate ρ S so noextrapolation is needed to extract GW fluxes at infinity.Another advantage is that the coordinate light speed at In this section ρ denotes the radial HH-coordinate. the boundaries vanishes, therefore, no numerical bound-ary condition must be imposed. After the decompositioninto azimuthal m -modes ψ = (cid:80) m ψ m e imϕ the equationin (2 + 1)-dimensional form reads (cid:0) C ττ ∂ τ + C τρ ∂ τ ∂ ρ + C ρρ ∂ ρ + C θθ ∂ θ + C τ ∂ τ + C ρ ∂ ρ + C θ ∂ θ + C ) ψ m = S s , (92)where the coefficients C ττ , C τρ , . . . are functions of ρ and θ and S s is the source term for spinning particle discussedin [30].The source term consists of derivatives of delta func-tions up to third order. For accurate results proper rep-resentation of delta functions must be used. Approxi-mation as Gaussian function and piecewise polynomialsas described in [44] were implemented to the Teukode .According to [41], piecewise polynomial approximationis more accurate for circular equatorial orbits and fasterto calculate than Gaussian approximation, whereas cal-culations with Gaussian approximation are more stablewhen the particle is moving in ρ or θ direction. Thethird derivative of the delta function, which is neededfor spinning particle, was implemented only as Gaussianapproximation in the previous works. In our work weintroduced to Teukode an approach suggested in [45],which describes slightly different formulas for piecewisepolynomial approximation to construct delta functionand its derivatives.
Teukode has been tested extensivelyon circular equatorial orbits of a spinning particle in[29, 30, 46–48], but in this work it is tested for the firsttime on eccentric equatorial orbits of a spinning particle.
B. Numerical results
This Section discusses our numerical calculations ofGW fluxes in the frequency domain (as described inSec. III A 1) and compare them with time domain resultsobtained from the
Teukode (Sec. III A 2).First we present our approach to calculate quantitiesrelated to an orbit for given parameters ˆ a , σ , p and e .These quantities include the energy and the angular mo-mentum from Eqs. (38) and (39) respectively, the orbitalfrequencies ˆΩ r and ˆΩ φ from Eqs. (52) and (53) respec-tively and the functions ˆ t ( χ ) and φ ( χ ) from Eqs. (55)and (56) respectively. The integrals (52) and (53) werecalculated numerically using methods built-in to Mathe-matica. We used extended precision to 48 places, becausehigh precision of the parameters a and ω = m Ω φ + n Ω r is needed for the calculation of the radial function R ± lmω .1 - - π π π π - - χ FIG. 3. The real part of exp( iD r ( ω mn t ( χ ) − mφ ( χ ))) fororbital parameters ˆ a = 0 . σ = − . p = 12, e = 0 . m = 2, n = 15 (top panel) and for orbital parameters ˆ a = 0 . σ = − . p = 12, e = 0 . m = 2, n = 4 (bottompanel). The red dots indicate the values at which the functionis calculated during the numerical integration. To calculate the energy and angular momentum fluxesand the strain at infinity, one has to find the partialamplitudes ˆ C ± lmn and Eq. (79) implies integration over χ . The numerical integration errors depend on the em-ployed integration method and the number of points atwhich the function is enumerated. For our purposes,a fractional accuracy of the order of 10 − is sufficient.Therefore, we used the midpoint rule inducing an er-ror of the order O ( N − ) to the integration, where N is the number of points. The advantage of the mid-point rule is that for given accuracy, this method mini-mizes the number of points N needed for the calculation.However, more complex method can be implemented inthe future to improve the accuracy of this integration.The main oscillatory part of Eq. (79) is contained in theexponential term exp( iD r ( ω mn t ( χ ) − mφ ( χ ))). Fig. 3shows the behavior of this oscillatory part for certainsetups. The higher the value of n is, the more theexponential function oscillates. These oscillations withhigh frequency are present especially around χ = π inhigh eccentricity cases. The number of the points N needed for the integration is calculated from the max-imum of the derivative of the function ω mn t ( χ ) − mφ ( χ )with respect to χ , which in dimensionless quantities reads(ˆ ω mn V t (ˆ r ( χ )) − mV φ (ˆ r ( χ ))) (cid:112) (1 − e ) /J ( χ ) /p .The radial functions R ± lmn were calculated usingthe BHPToolkit [36], which employs the Mano-Suzuki-Takasugi (MST) method [49] or a numerical integra-tion of the radial TE. The angular functions S ˆ a ˆ ωlm werealso calculated using the BHPToolkit which employs theLeaver’s method [50].The strain is calculated from Eq. (81) and the fluxesare calculated from Eqs. (88). The range of l and n for l l l l - - - - - n C lmn + FIG. 4. Absolute values of the partial amplitudes (cid:12)(cid:12)(cid:12) ˆ C + lmn (cid:12)(cid:12)(cid:12) fororbital parameters ˆ a = 0 . σ = − . p = 12, e = 0 . m = 2 given m -mode was found in the following way. First wecalculate the coefficient ˆ C + lmn for l = max( | m | ,
2) for arange of n to find the mode with the maximal (cid:12)(cid:12)(cid:12) ˆ C + lmn (cid:12)(cid:12)(cid:12) .Then, we calculate other l and n modes until the abso-lute value is less than a chosen accuracy times the maxi-mal mode. In our calculations, we have chosen accuracy10 − . However, in some cases the absolute value (cid:12)(cid:12)(cid:12) ˆ C + lmn (cid:12)(cid:12)(cid:12) is not monotonous in n and it drops suddenly for some n .Because of this, after such a sudden decrease, amplitudesfor more n must be calculated. In Fig. 4, the absolutevalues of the coefficients (cid:12)(cid:12)(cid:12) ˆ C + lmn (cid:12)(cid:12)(cid:12) are plotted for an orbitwith ˆ a = 0 . σ = − . p = 12, e = 0 . m = 2 for different l and n . We can see that, forgiven accuracy, only limited number of modes is needed(for l = m = 2 it is 21 n -modes) and the absolute valueof the amplitudes is decreasing exponentially with | n | forsufficiently high | n | . Note that although the astrophys-ical relevant value of the spin σ is of the same order asthe mass ratio q (cid:28)
1, it is possible to calculate the GWfluxes for higher spins and then linearize the result in σ to find the contribution of spin σ (cid:28)
1. We use also theselarge values to make any deficiencies in our calculationsprominent.In Appendix C we compare our coefficients ˆ C ± lmn andfluxes F E ∞ lmn and F E H lmn with that of [25]. A simplifiedversion of our code calculating GW fluxes from circu-lar equatorial orbit of a spinning particle around a KerrBH was used to independently verify the results of [29].These results are discussed in detail in [51]. Tables ofthe values of the partial amplitudes ˆ C ± lmn for future ref-erences are in Appendix D.
1. Comparison of frequency domain and time domain
To compare the time domain and the frequency domainresults, we have calculated the coefficients ˆ C + lmn for somerange of l and n in the frequency domain for different2values of the spin σ and of the eccentricity e . We haveused these coefficients to find the respective strains andenergy fluxes at infinity. Then, these results have servedas reference values in our comparison with the azimuthal m -mode of the strain at infinity multiplied by the radialcoordinate ˆ r h m and the energy fluxes at infinity F E ∞ m obtained in the time domain. Because of the fact thatthe space discretization applied in Teukode induces nu-merical errors to the time domain calculations, we haverun the time domain calculations for several resolutionsand tested the convergence of the code.To calculate the strains and the fluxes with in the timedomain with
Teukode , we need to approximate the deltafunctions representing the secondary body in the ρ and θ directions. To do that we have used different combina-tions of Gaussian functions and piecewise polynomials inthese directions. The accuracy appears to be higher whenthe piecewise polynomial are used in both ρ and θ direc-tion or Gaussian function in ρ direction and piecewisepolynomial in θ direction, than in the other two possiblesettings, i.e. Gaussian in both directions and Gaussianin θ direction with piecewise polynomial in ρ direction.When the piecewise polynomial is used in both directions,calculations are faster and, therefore, we have used thisapproximation in most cases. In our calculations, thestrain has been extracted at r = ∞ and θ = π/ T r starting at the retarded coordinate around u = 350 M ,where u = t − r ∗ .In order to provide a first comparison of the frequencyand the time domain results, we use the relative differ-ence of the azimuthal mode m of the strain at r = ∞ and θ = π/ δh m = (cid:12)(cid:12)(cid:12)(cid:12) − h td m h fd m (cid:12)(cid:12)(cid:12)(cid:12) , (93)where h td m is the strain calculated using Teukode and h fd m is m -mode of the strain calculated in frequency domainusing Eq. (81) without the sum over m . Fig. 5 shows therelative difference of the azimuthal modes m = 1 , , , u . In this plot, the strain calculated in the frequencydomain (the denominator of δh m ) remains fixed, whileeach time domain calculated evolution of the strain isperformed for different number of points in the ρ di-rection N ρ (resolution). The delta function is approx-imated by a piecewise polynomial for five resolutions( N ρ =1200,1704,2400,3384,4800), while in one case is ap-proximated by a Gaussian function for N ρ = 4800. Wecan see that the relative difference δh m tends to decreaseas the resolution increases, but for the highest resolu-tion N ρ = 4800 the numerical noise becomes significant.Though the Gaussian approximation is less accurate, theamplitude of its noise is relatively smaller than the am-plitude of the noise for the piecewise polynomial approx-imation with the same resolution. We speculate that thecause of this numerical noise comes from the fact that asthe resolution increases, the approximation becomes less × - δ h × - × - × - δ h × - × - × - δ h
200 400 600 800 1000 1200 14005. × - × - × - u δ h FIG. 5. The relative difference of the strain δh m from the m = 1 mode (top panel) to the m = 4 mode (bottom panel) asa function of the retarded coordinate ˆ u at r = ∞ and θ = π/ ρ direction N ρ . The piecewise polynomialapproximation of the delta function was used for all casesapart from one, for which the Gaussian approximation withresolution 4800 was employed. The parameters of the orbitare ˆ a = 0 . σ = − . p = 12, e = 0 .
2. The initial noise iscaused by zero initial data in time domain. smooth. Namely, we have used a 12th order approxima-tion of the delta function, which is 12 points wide, for3each resolution; therefore, the higher the resolution is,the narrower and higher is the delta function. Note thatthe m = 1-mode has very small value and the noise hasrelatively higher amplitude than in m = 2 , , m = 0-mode, which is not shown here, althoughnon-zero, has extremely small value allowing the numer-ical noise to be dominant. σ - σ - σ - σ σ σ × - Δρδ F E ∞ N ρ FIG. 6. The relative difference of the energy flux δ F E ∞ lm ofthe l = m = 2 mode as function of the grid length in the ρ direction of the time domain calculations. Note that the timedomain calculations have been projected on the Y lm basis,while the frequency domain ones on the S aωlm . Each curve rep-resents a different value of the secondary spin, while the Kerrparameter ˆ a = 0 .
9, semi-latus rectum p = 12 and eccentricity e = 0 . To further check our results, we have calculated therelative difference of the energy fluxes δ F E ∞ lm = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − F E ∞ lm, td F E ∞ lm, fd (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (94)where F E ∞ lm, td is the value calculated using Teukode and F E ∞ lm, fd is the value calculated with the frequency domainapproach summed over n . Fig. 6 shows how the timedomain calculations of the dominant l = m = 2 modeof the energy fluxes converges to the frequency ones asthe resolution increases. For this plot we have keptfixed the Kerr parameter ˆ a = 0 .
9, the semi-latus rec-tum p = 12 and the eccentricity e = 0 .
2, while we haveused for each curve a different value of the secondaryspin σ spanning from − . .
5. The relative differ-ence in the fluxes should converge to zero as the gridlength ∆ ρ = ( ρ S − ρ + ) /N ρ decreases (increasing resolu-tion). However, the relative differences do not convergeto zero, because in the frequency domain calculations weuse the projection to spin-weighted spheroidal harmonics S aωlm and Teukode projects the strain to the spin-weightedspherical harmonics Y lm = S lm . For the dominant modethe difference between the projections to these functionsis low because for low aω , the spheroidal functions S aωlm can be approximated by the spherical functions Y lm . × - δ F E ∞ N ρ σ - σ - σ - σ σ σ × - × - × - δ F E ∞ × - × - × - Δρδ F E ∞ FIG. 7. Comparison of frequency domain and time domainresults. The relative difference δ F E ∞ m =1 (top panel), δ F E ∞ m =2 (middle panel) and δ F E ∞ m =3 (bottom panel) is plotted for dif-ferent values of the secondary spins σ spanning from − . .
5. The Kerr parameter ˆ a = 0 .
9, the semi-latus rectum p = 12 and the eccentricity 0 . Because of the aforementioned projection issue, for aproper comparison of the time and frequency domain re-sults, we must calculate the sum of the fluxes over l . Therelative difference δ F E ∞ m = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − F E ∞ m, td F E ∞ m, fd (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (95)for m = 1 , , σ in the frequency domain and intime domain we used different resolutions ( N ρ =1200 , , , , ρ corresponding to the highest resolu-tion N ρ = 4800 shows variance in the relative differences.4This is caused by the fact that the noise amplitude is thehighest for the highest resolution, which can be seen inFig. 5. Especially in the case m = 1 where the energyflux is significantly lower than for m = 2, the variancein the relative differences is clearly visible. For the high-est resolution, the relative difference is higher for highervalues of spin | σ | . This can be caused by the numericalnoise in time domain calculations induced by the non-smoothness of the piecewise polynomial approximationof the third derivative of the delta function. Namely, theterm with the third derivative is proportional to the spin σ . To check the dependence of our calculations on thevalue of eccentricity, we have calculated the energy fluxesfor fixed Kerr parameter ˆ a = 0 .
9, secondary spin σ = 0 . p = 12, while the eccentricity e value span from 0 . .
8. For each eccentricity wehave calculated the relative difference in the energy fluxes δ F E ∞ m for m = 1 , ,
3. Then, we have compared thedependence of the relative difference on the resolutionfor different eccentricities as in the case with the chang-ing secondary spin. This comparison is shown in Fig. 8.First, we have calculated the dominant m = 2 mode intime domain with piecewise polynomial approximation ofthe delta function in both ρ and θ direction (p-p), butfor e = 0 . δ F E ∞ does not converge to zero (purple line in themiddle panel of Fig. 8). Therefore, for m = 2 and othermodes, we performed the time domain calculations for e = 0 . ρ directionand the piecewise polynomial approximation in θ direc-tion (G-p, red line in all panels of Fig. 8). However, the m = 1 mode has low amplitude and the noise is thereforemore significant and the p-p approximation for e = 0 . e = 0 . m = 1 mode we repeated the calculation for e = 0 . e = 0 . δ F E ∞ m converges tozero for all m -modes and eccentricities e . It is also appar-ent that even the piecewise polynomial approximation isin general more efficient, it has its own limitation thatin our example became prominent for high eccentricitiesand modes with small values, i.e. in modes that the nu-merical noise is dominant. IV. SUMMARY
In this work, we have studied the motion of a spinningparticle in the equatorial plane of a Kerr black hole andthe GW fluxes from these orbits. The only possible con-figuration of the spins in this setup is the spins to be par-allel or antiparallel. In this framework, we have deriveda reduced set of equations of motion equivalent to theMPD equations with TD SSC. Taking advantage of thefact that an orbit can be characterized by its constants e = - p e = - G - δ F E ∞ N ρ e e e e = - p e = - p × - × - δ F E ∞ × - × - × - Δρδ F E ∞ FIG. 8. Comparison of frequency domain and time domainresults. The relative difference δ F E ∞ m =1 (top panel), δ F E ∞ m =2 (middle panel) and δ F E ∞ m =3 (bottom panel) is plotted for dif-ferent values of the eccentricity e spanning from 0 . . a = 0 .
9, the secondary spin σ = 0 . p = 12 are kept fixed for all thecases. If not specified, the delta function is approximated bya piecewise polynomial in both ρ and θ direction. For m = 1, e = 0 . m = 2, e = 0 . ρ direction and piecewise polynomialin θ direction. For m = 1, e = 0 . ρ and θ directions. of motion, namely the energy E and the z componentof the total angular momentum J z , we have providedexplicit formulas for the energy and the angular momen-tum in terms of the eccentricity e and semi-latus rectum p . Furthermore, through the reduced equations of mo-tion and by introducing a Mino-like time parameter λ ,we were able to find expressions allowing the numericalcalculation of the frequencies of the radial and azimuthalmotion. These expressions provide the frequencies with5respect to λ or the BL time.The orbital findings were then implemented in the cal-culation of the GW fluxes from the equatorial orbits inthe frequency domain. Namely, this work introduces theformulas giving the strain h , the energy fluxes and theangular momentum fluxes at infinity and at the hori-zon from a spinning secondary moving on the equato-rial plane of a Kerr black hole. For this purpose, wehave developed a Mathematica code calculating the am-plitudes C ± lmn on which the frequency domain GW fluxesdepend. We plan to make this code publicly availablethrough the Black Hole Perturbation Toolkit repository.The frequency domain results were, then, compared withtime domain results obtained from a TE solver called Teukode . To improve the efficiency of
Teukode , we haveimplemented a piecewise polynomial to approximate thedelta functions and its derivatives in the spinning-particlesource term. The comparison has shown good agreementbetween the frequency domain results with the time do-main ones.To check the discretization error in the time domaincalculations introduced by the piecewise polynomial, wehave calculated the fluxes in time domain for differentresolutions and compared them with the respective fre-quency domain results. The difference between the re-sults from these two approaches tend to consistently de-crease with increasing resolution. However, for the high-est resolution, which we have implemented, the numeri-cal noise in the time domain calculations becomes signif-icant. This behavior occurs for different calculation se-tups. Namely, we have checked our calculations by vary-ing the secondary spin while keeping the other parame-ters fixed and by varying the eccentricities while keepingthe other parameters fixed.These calculations are part of the on-going effort tobuild post-adiabatic gravitational waveforms modellinggravitational waves emitted by extreme mass ratio inspi-rals. In a future work, the frequency domain fluxes willbe used to find the adiabatic evolution of the orbit onthe equatorial plane under the influence of radiation re-action. The influence of the secondary spin on the changeof the orbital parameters and phase of the GW will bestudied.
ACKNOWLEDGMENTS
The authors have been supported by the fellowshipLumina Quaeruntur No. LQ100032102 of the CzechAcademy of Sciences. The authors would like to acknowl-edge networking support by the GWverse COST ActionCA16104, “Black holes, gravitational waves and funda-mental physics. V.S. would also like to express gratitudefor the hospitality of the Theoretical Physics Instituteat the University of Jena. We would like to thank Se-bastiano Bernuzzi, Enno Harms, Vojtˇech Witzany andTom´aˇs Ledvinka for useful discussions and comments.This work makes use of the Black Hole Perturbation Toolkit. Computational resources were supplied by theproject ”e-Infrastruktura CZ” (e-INFRA LM2018140)provided within the program Projects of Large Research,Development and Innovations Infrastructures.
Appendix A: List of dimensionless quantities
Throughout this work, we use several quantities bothin their full form and dimensionless form. The dimen-sionless form is denoted by a hat. Their list with relationbetween the full and dimensionless form is in Table I.Some quantities such as the time parameter λ or x aredefined only as dimensionless whereas other quantitiesare used only in their full form. ˆ t = t/M BL timeˆ r = r/M BL radial coordinateˆ a = a/M Kerr parameter σ = S/ ( µM ) Secondary spinˆ E = E/µ
Energyˆ J z = J z / ( µM ) Angular momentumˆ τ = τ /M Proper timeˆ∆ = ∆ /M ˆ (cid:36) = (cid:36) /M ˆΩ r = Ω r M Radial BL frequencyˆΩ φ = Ω φ M Orbital BL frequencyˆ ω = ωM Frequencyˆ C ab = C ab /µ ˆ C σab = C σab /µ ˆ C ± lmn = C ± lmn M /µ Partial amplitudesˆ u = u/M Retarded coordinateTABLE I. List of dimensionless quantities
Appendix B: Formulas for GW fluxes
In this Appendix we derive the coefficients A i = A i ( r, θ ) and B i +1 = B i +1 ( r, θ ), i = 0 , ,
2, in Eq. (73)for calculation of partial amplitudes of GWs from generalbound orbits of a spinning particle around a Kerr blackhole. Then we list explicit formulas for equatorial orbitswith secondary spin parallel to the z axis.To find the form of the coefficients A i and B i +1 inEq. (73), the integrals (68) and (65) must be evaluatedusing rules for integrating delta functions. We can clas-sify the parts of the coefficients A i according to termfrom which they originate: A = (cid:88) ab = nn,nm,mm ( A ab + A tφab + A rab + A θab ) , (B1) A = (cid:88) ab = nm,mm ( A ab + A tφab + A rab + A θab ) , (B2) A = A mm + A tφmm + A rmm + A θmm . (B3)The terms A abi originate from the first term of the stress-energy tensor (70) containing the non-spinning part of6 T µν and parts containing Christoffel symbols. The terms A tφabi originate from the second term of (70) containing t and φ derivative. Similarly, the terms A rabi or A θabi originate from the second term of Eq. (70) containing r or θ derivative respectively. The subscripts ab denote thetetrad legs in Eq. (69). A abi can be found by integrating θ and φ after substi-tuting the first term of Eq. (70) into Eq. (68) by replacing θ → θ p ( t ), φ → φ p ( t ) and then using integration by partsin Eq. (65), where the derivatives with respect to r in(69) are shifted to the radial function R ± lmn to obtain A abi = (cid:0) C ab − C σab (cid:1) f ( i ) ab , (B4)where C ab and C σab are defined in Eq. (71).To find the form of A tφabi , we must perform integrationby parts in Eq. (68) where the t or φ derivative in thesecond term of Eq. (70) are shifted to exp( iω − imφ )because no other functions depend on t and φ . Fromthis, we get terms multiplied by iω and − imφ . Afterthat, an integration over r of Eq. (65) is done similarlyas in the previous case and we obtain A tφabi = d τ d t ( iωS tµ − imS φµ ) v ν e ( a )( µ e ( b ) ν ) f ( i ) ab . (B5)The term A θabi is derived in similar way. The derivativewith respect to θ in the second term in (70) is shifted tothe functions f ( i ) ab and the tetrad legs. The boundaryterm vanishes because f ( i ) ab ( r,
0) = f ( i ) ab ( r, π ) = 0. Thefinal term has the form A θabi = d τ d t S θ ( µ v ν ) f ( i ) ab ∂ θ ( e ( a ) µ e ( b ) ν )+ d τ d t S θ ( µ v ν ) e ( a ) µ e ( b ) ν ∂ θ f ( i ) ab . (B6)Now let us focus on the term containing the r deriva-tive in Eq. (70). After substituting the stress-energy ten-sor (70) into Eq. (69), the derivative of the delta functioncan be shifted to the function f ( i ) ab and the tetrad legs. Forexample, from the first term of T nm we obtain ∂ r (cid:16) f (1) nm ( r, θ ) n µ m ν ∂ r (cid:16) ( v t ) − S r ( µ v ν ) δ (cid:17)(cid:17) = ∂ r (cid:16) f (1) nm ( r, θ ) n µ m ν ( v t ) − S r ( µ v ν ) δ (cid:17) − ∂ r (cid:16) ∂ r (cid:16) f (1) nm ( r, θ ) n µ m ν (cid:17) ( v t ) − S r ( µ v ν ) δ (cid:17) (B7)After substituting Eq. (68) into Eq. (65) we can changethe order of the t and r integrals and integrate by parts.From the second term in Eq. (B7) we obtain a term withderivatives with respect to r of f ( i ) ab and the tetrad legs A rabi = d τ d t S r ( µ v ν ) f ( i ) ab ∂ r ( e ( a ) µ e ( b ) ν )+ d τ d t S r ( µ v ν ) e ( a ) µ e ( b ) ν ∂ r f ( i ) ab . (B8) From the second term in Eq. (B7) we obtain terms withone order higher derivatives of the radial function R ± lmω ,of which the integration by parts we can perform to ob-tain the coefficients B = (cid:88) ab = nn,nm,mm B ab , (B9) B = (cid:88) ab = nm,mm B ab , (B10) B = B mm , (B11)where B ab ( i +1) = − d τ d t S r ( µ v ν ) e ( a ) µ e ( b ) ν f ( i ) ab . (B12)The functions f ( i ) ab = f ( i ) ab ( r, θ ) in the equatorial planeare given by f (0) nn (cid:16) r, π (cid:17) = − r ∆ (cid:18) L † L † − iar L † (cid:19) S aωlm ( θ ) (cid:12)(cid:12)(cid:12) θ → π , (B13) f (0) nm (cid:16) r, π (cid:17) = 2 √ r ∆ (cid:18) iK ∆ + 2 r (cid:19) L † S aωlm ( θ ) (cid:12)(cid:12)(cid:12) θ → π , (B14) f (1) nm (cid:16) r, π (cid:17) = 2 √ r ∆ L † S aωlm ( θ ) (cid:12)(cid:12)(cid:12) θ → π , (B15) f (0) mm (cid:16) r, π (cid:17) = (cid:18) i∂ r (cid:18) K ∆ (cid:19) − i K ∆ r + K ∆ (cid:19) S aωlm (cid:16) π (cid:17) , (B16) f (1) mm (cid:16) r, π (cid:17) = − (cid:18) r + i K ∆ (cid:19) S aωlm (cid:16) π (cid:17) , (B17) f (2) mm (cid:16) r, π (cid:17) = − S aωlm (cid:16) π (cid:17) , (B18)where K = ( r + a ) ω − am , (B19) L † n = ∂ θ − m csc θ + aω sin θ + n cot θ . (B20)Up to this point the analysis holds for generic orbits ofa spinning particle. When we constrain the particle onequatorial orbits with its spin set parallel to the z axis,then S θµ = 0 for all µ and, therefore, A θabi = 0. For thepresentation of the equatorial case, we prefer to use thedimensionless quantities.In the definition of C ab and C σab (71) we can replacethe derivative with respect to τ in v µ with derivativewith respect to λ and use the fact that V t = dˆ t/ d λ , V r = dˆ r/ d λ and V φ = d φ/ d λ . From Eqs. (23) and (25)we obtain ˆ C nn = d λ dˆ t V n Σ σ , (B21)ˆ C nm = d λ dˆ t V m V n (cid:0) r + σ (cid:1) σ (ˆ r + 2 σ ) , (B22)ˆ C mm = d λ dˆ t ˆ rV m (ˆ r + 2 σ ) , (B23)7ˆ C σnn = d λ dˆ t σ r Σ σ (cid:32) aV n − ˆ∆(ˆ r + 2 σ )ˆ r − σ V n x + (ˆ a − ˆ r ) V r x (cid:33) , (B24)ˆ C σnm = d λ dˆ t iσ √ r Σ σ (cid:32) − ˆ a − ˆ r ∆ (cid:16) V n + ( V r ) (cid:17) − ˆ a − ˆ r ˆ∆ V n V r + 3ˆ aσ ˆ r Σ σ V n x − ˆ aV r x + ˆ∆(ˆ r + 2 σ )2ˆ r Σ σ x (cid:33) , (B25)ˆ C σmm = d λ dˆ t σ Σ σ (cid:18) − ˆ a ˆ∆ (cid:16) V n + 2 V n V r + ( V r ) (cid:17) + i ˆ r √ V n + V r ) V m (cid:19) , (B26)where V n = V t n t + V r n r + V φ n φ M = − P σ (ˆ r ) + V r , (B27) V m = V t m t + V φ m φ M = − ix (cid:0) ˆ r + 2 σ (cid:1) √ σ . (B28)We can rewrite the expressions for A tφabi , A rabi and B ab ( i +1) into dimensionless quantities asˆ A tφabi = d λ dˆ t (cid:16) i ˆ ω ˆ S t ( a − im ˆ S φ ( a (cid:17) V b ) ˆ f ( i ) ab (cid:16) ˆ r, π (cid:17) , (B29)ˆ A rabi = d λ dˆ t ( ˆ S r ( ∂ ˆ r a V b ) + ˆ S r ( b V ∂ ˆ r a ) ) f ( i ) ab (cid:16) ˆ r, π (cid:17) + d λ dˆ t ˆ S r ( a V b ) ∂ ˆ r ˆ f ( i ) ab (cid:16) ˆ r, π (cid:17) , (B30)ˆ B ab ( i +1) = − d λ dˆ t ˆ S r ( a V b ) ˆ f ( i ) ab (cid:16) ˆ r, π (cid:17) , (B31)where we used the dimensionless projections of S µν intothe tetradˆ S tn = 1 µM (cid:0) S tr n r + S rφ n φ (cid:1) = σ (cid:0) x ˆ (cid:36) − aV n (cid:1) r Σ σ , (B32)ˆ S rn = 1 µM (cid:0) − S tr n t + S rφ n φ (cid:1) = − σx ˆ∆2ˆ r Σ σ , (B33)ˆ S φn = 1 µ (cid:0) − S tφ n t − S rφ n r (cid:1) = σ (ˆ ax − V n )2ˆ r Σ σ , (B34)ˆ S tm = 1 µM S tφ m φ = − iσ ˆ (cid:36) V r √ σ , (B35)ˆ S rm = 1 µM (cid:0) − S tr m t + S rφ m φ (cid:1) = − iσP σ (ˆ r ) √ σ , (B36)ˆ S φm = − µ S tφ m t = − iσ ˆ aV r √ σ . (B37)The quantities V ∂ ˆ r a and ˆ S r∂ ˆ r a can be understood asdimensionless projections on the differentiated tetrad ∂ r e µ ( a ) V ∂ ˆ r n = M (cid:18) V t ∂ r n t + V r ∂ r n r + V φ ∂ r n φ M (cid:19) = (ˆ a − ˆ r ) P σ (ˆ r )ˆ r ˆ∆ , (B38) V ∂ ˆ r m = M (cid:18) V t ∂ r m t + V φ ∂ r m φ M (cid:19) = V m ˆ r − i √ aP σ (ˆ r )ˆ∆ , (B39)ˆ S r∂ ˆ r n = 1 µ (cid:0) − S tr ∂ r n t + S rφ ∂ r n φ (cid:1) = σ (ˆ a − ˆ r ) x ˆ r Σ σ , (B40)ˆ S r∂ ˆ r m = 1 µ (cid:0) − S tr ∂ r m t + S rφ ∂ r m φ (cid:1) = − iσ (2ˆ ax + P σ (ˆ r )) √ r Σ σ , (B41)where the covariant components of the tetrad are n µ = 12Σ (cid:0) − ∆ , − Σ , , a ∆ sin θ (cid:1) , (B42) m µ = − ρ √ (cid:0) ia sin θ, , Σ , − i(cid:36) sin θ (cid:1) . (B43) Appendix C: Comparison with [Phys. Rev. D 73,024027 (2006)]
This section compares our frequency domain calcula-tions for a non-spinning particle with results obtained in[25]. In that work, the GW fluxes were calculated fromgeneric orbits of a non-sinning particle moving arounda Kerr black hole using Teukolsky formalism with thefractional accuracy of the energy flux l, m -modes set to10 − .We have compared our data with theirs for an equa-torial orbit around a Kerr black hole with a = 0 . p = 8 . . r ISCO and e = 0 .
3. In particular, wehave compared our energy fluxes F E ∞ lmn , F E H lmn and ampli-tudes ˆ C ± lmn with their data. In the top panel of Fig. 9,we plot the difference between our calculated fluxes F E ∞ lmn and the fluxes F E ∞ lmn DH calculated in [25] normalized by8 - - - - - δ F l m nE ∞ l m l m l m - - - - - - - - n δ C l m n + FIG. 9. Differences between our frequency domain resultsand the results obtained in [25]. Top panel: The difference δ F E ∞ lmn between the fluxes normalized by max n min ≤ n ≤ n max (cid:12)(cid:12) F E ∞ lmn (cid:12)(cid:12) .Bottom panel: The difference δ ˆ C + lmn between the coefficientsnormalized by max n min ≤ n ≤ n max (cid:12)(cid:12)(cid:12) ˆ C + lmn (cid:12)(cid:12)(cid:12) . the maximum of F E ∞ lmn over n for each lm -mode δ F E ∞ lmn = (cid:12)(cid:12) F E ∞ lmn − F E ∞ lmn DH (cid:12)(cid:12) max n min ≤ n ≤ n max F E ∞ lmn . (C1)We can see that for each lmn -mode the error is less than10 − of the maximal value for given l and m . In a similarway, we have compared the coefficients ˆ C + lmn using thequantity δC + lmn = (cid:12)(cid:12)(cid:12) ˆ C + lmn − ˆ C + lmn DH (cid:12)(cid:12)(cid:12) max n min ≤ n ≤ n max (cid:12)(cid:12)(cid:12) ˆ C + lmn (cid:12)(cid:12)(cid:12) . (C2)The result of this comparison is shown in the bottompanel of Fig. 9. The normalized difference for the coeffi-cients ˆ C + lmn is higher than in the flux comparison, becausethe flux is calculated from the second power of ˆ C + lmn andthe error is thus relatively smaller. Similar comparisonwas calculated for the horizon fluxes F E H lmn and ˆ C ± lmn .The result is shown in Fig. 10. Although the accuracy isless than 10 − for some modes, the contribution from thehorizon fluxes is smaller than from the fluxes to infinityand the overall accuracy remains higher. - - - - - - δ F l m nEH l m l m l m - - - - - - - - n δ C l m n - FIG. 10. Differences between our frequency domain resultsand the results obtained in [25]. Top panel: The difference δ F E H lmn between the fluxes normalized by max n min ≤ n ≤ n max (cid:12)(cid:12) F E H lmn (cid:12)(cid:12) .Bottom panel: The difference δ ˆ C − lmn between the coefficientsnormalized by max n min ≤ n ≤ n max (cid:12)(cid:12)(cid:12) ˆ C − lmn (cid:12)(cid:12)(cid:12) . Appendix D: Data tables
In this appendix we present data tables of the partialamplitudes C ± lmn (Tables II to V) for an orbit with or-bital parameters ˆ a = 0 . σ = − . p = 12, e = 0 . E = 0 . . . . ˆ J z = 3 . . . . ˆΩ φ = 0 . r = 0 . (cid:12)(cid:12) C + lmn (cid:12)(cid:12) > − are listed for 1 ≤ m ≤
4. The accuracy of the dominant modes should be atsix significant digits, but for lower modes, the accuracydrops. This accuracy depends mostly on the accuracy ofthe radial function R ± lmn and the coordinates t ( χ ) and φ ( χ ).9 TABLE II. List of partial amplitudes C ± l n for an orbit with orbital parameters ˆ a = 0 . σ = − . p = 12, e = 0 . l m n Re { C + lmn } Im { C + lmn } Re { C − lmn } Im { C − lmn } × − -1.467715 × − × − × − × − -4.271458 × − × − × − × − -1.043309 × − -5.859804 × − × − × − -1.824696 × − -1.435506 × − × − × − -1.244299 × − -1.084857 × − × − × − × − -6.110857 × − × − × − -5.100353 × − -4.159829 × − × − × − -6.201949 × − -3.035062 × − × − × − -3.735099 × − -1.503776 × − × − × − -1.633602 × − -6.182375 × − × − × − -5.891306 × − -2.272133 × − × − × − -1.849669 × − -7.750262 × − × − × − -5.168866 × − -2.508722 × − × − × − -1.283265 × − -7.819832 × − × − × − -2.742169 × − -2.371637 × − × − × − -9.135246 × − × − × − × − -2.513096 × − × − × − × − -5.375419 × − × − × − × − -6.372160 × − × − × − × − × − × − × − × − -6.663590 × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − -3.354890 × − × − -9.269538 × − × − × − × − -7.573044 × − × − × − × − -4.275336 × − × − × − × − -1.991900 × − × − × − × − -8.225436 × − × − × − × − -3.122972 × − × − × − -1.963261 × − -1.023901 × − × − -2.593480 × − -8.203217 × − -2.319941 × − × − -1.783264 × − -3.940505 × − -8.863485 × − [1] Pau Amaro-Seoane, Heather Audley, Stanislav Babak,John Baker, Enrico Barausse, Peter Bender, EmanueleBerti, Pierre Binetruy, Michael Born, Daniele Bortoluzzi,Jordan Camp, Chiara Caprini, et al. Laser InterferometerSpace Antenna. arXiv e-prints , page arXiv:1702.00786,February 2017.[2] Stanislav Babak, Jonathan Gair, Alberto Sesana, EnricoBarausse, Carlos F. Sopuerta, Christopher P. L. Berry,Emanuele Berti, Pau Amaro-Seoane, Antoine Petiteau,and Antoine Klein. Science with the space-based inter-ferometer LISA. V. Extreme mass-ratio inspirals. Phys.Rev. D , 95:103012, May 2017.[3] Tanja Hinderer and ´Eanna ´E. Flanagan. Two-timescaleanalysis of extreme mass ratio inspirals in Kerr space-time: Orbital motion.
Phys. Rev. D , 78:064028, Sep 2008. [4] Yasushi Mino, Misao Sasaki, Masaru Shibata, HideyukiTagoshi, and Takahiro Tanaka. Black hole perturbation:Chapter 1.
Prog.Theor.Phys.Suppl. , 128:1–121, 1997.[5] Sarp Akcay, Sam R. Dolan, Chris Kavanagh, JordanMoxon, Niels Warburton, and Barry Wardell. Dissipationin extreme mass-ratio binaries with a spinning secondary.
Phys. Rev. D , 102(6):064013, September 2020.[6] Michael D. Hartl. Dynamics of spinning test particles inKerr space-time.
Phys. Rev. , D67:024005, 2003.[7] J. Frenkel. Die elektrodynamik des rotierenden elektrons.
Zeitschrift fuer Physik , 37(4-5):243–262, 1926.[8] Cornel Lanczos. Ueber eine invariante formulierung dererhaltungssaetze in der allgemeinen relativitaetstheorie.
Zeitschrift fuer Physik , 59(7-8):514–539, 1930. TABLE III. List of partial amplitudes C ± l n for the same orbit as in Table II. l m n Re { C + lmn } Im { C + lmn } Re { C − lmn } Im { C − lmn } × − -2.613368 × − × − × − × − -3.499305 × − × − × − × − -1.253110 × − × − × − × − × − -9.742181 × − -2.076711 × − × − -8.824264 × − × − × − × − -1.089118 × − × − × − × − -7.404687 × − × − × − × − -3.764006 × − × − × − × − -1.599176 × − × − × − × − -5.997531 × − × − × − × − -2.047506 × − × − × − × − -6.480880 × − × − × − × − -1.923652 × − × − × − × − -5.403652 × − × − × − × − -1.460369 × − × − × − × − -3.880442 × − × − × − × − -9.455299 × − × − × − × − -1.713718 × − × − × − × − × − × − -9.847625 × − × − -1.111405 × − × − -1.653879 × − × − -1.883465 × − × − -1.566972 × − × − -1.615521 × − × − -9.383274 × − × − -9.843758 × − × − -4.500460 × − × − -4.845333 × − × − -1.888329 × − × − -2.055320 × − × − -7.241796 × − × − -7.798475 × − × − -2.603261 × − × − -2.708873 × − × − -8.912589 × − × − -8.743846 × − × − -2.937305 × − × − -2.649605 × − × − -9.389021 × − × − -7.609938 × − × − -2.926757 × − × − × − × − -4.205354 × − × − -3.437785 × − × − -4.056764 × − × − -1.313240 × − × − -3.881701 × − × − × − × − -2.357344 × − × − × − × − -1.142665 × − × − × − × − -4.822614 × − × − × − × − -1.851227 × − × − × − × − -6.628891 × − × − × − × − -2.250158 × − × − × − × − -7.319933 × − × − -2.652971 × − -9.489133 × − -1.219765 × − × − -8.273581 × − -1.010641 × − -1.194663 × − × − × − -7.026352 × − -7.612572 × − × − × − -3.875359 × − -3.833430 × − × − × − -1.841952 × − -1.655902 × − × − × − -7.882070 × − -6.404217 × − [9] Myron Mathisson. Neue mechanik materieller systemes. Acta Phys.Polon. , 6:163–2900, 1937.[10] Myron Mathisson. Republication of: New mechanics ofmaterial systems.
General Relativity and Gravitation ,42(4):1011–1048, 2010.[11] Achille Papapetrou. Spinning test particles in generalrelativity. 1.
Proc.Roy.Soc.Lond. , A209:248–258, 1951.[12] E. Corinaldesi and Achille Papapetrou. Spinning testparticles in general relativity. 2.
Proc. Roy. Soc. Lond. , A209:259–268, 1951.[13] W Tulczyjew. Motion of multipole particles in generalrelativity theory.
Acta Phys. Pol , 18:393, 1959.[14] W.G. Dixon. A covariant multipole formalism for ex-tended test bodies in general relativity.
Il Nuovo Ci-mento , 34(2):317–339, 1964.[15] W.G. Dixon. Dynamics of extended bodies in gen-eral relativity. II. Moments of the charge-current vector.
Proc.Roy.Soc.Lond. , A319:509–547, 1970. TABLE IV. List of partial amplitudes C ± l n for the same orbit as in Table II. l m n Re { C + lmn } Im { C + lmn } Re { C − lmn } Im { C − lmn } × − -7.760638 × − × − × − × − × − -3.882412 × − -1.922022 × − × − -4.030524 × − × − × − × − × − -9.633994 × − -5.865711 × − × − × − -1.135235 × − -7.650771 × − × − × − -7.199807 × − -5.358218 × − × − × − -3.482019 × − -2.851315 × − × − × − -1.438478 × − -1.289482 × − × − × − -5.363338 × − -5.227212 × − × − × − -1.863735 × − -1.957612 × − × − × − -6.164127 × − -6.902951 × − × − × − -1.969217 × − -2.321446 × − × − × − -6.141824 × − -7.513658 × − × − × − -1.884840 × − -2.356231 × − × − × − -5.716907 × − -7.198380 × − × − × − -1.717137 × − -2.149742 × − × − × − -5.246898 × − -6.287407 × − × − × − -2.488432 × − -2.373154 × − × − -3.145041 × − × − × − × − × − -2.124982 × − -2.014665 × − × − × − -2.891155 × − -2.717792 × − × − × − -2.242899 × − -2.082785 × − × − × − -1.324312 × − -1.210354 × − × − × − -6.609932 × − -5.923846 × − × − × − -2.939162 × − -2.573399 × − × − × − -1.200513 × − -1.023073 × − × − × − -4.593515 × − -3.795772 × − × − × − -1.668689 × − -1.331898 × − × − × − -5.810556 × − -4.462024 × − × − × − -1.953240 × − -1.437145 × − × − × − -6.478361 × − -3.587094 × − × − -1.123248 × − × − × − × − × − -3.389273 × − -7.785256 × − × − × − -4.536955 × − -1.859315 × − × − × − -3.512562 × − × − × − × − -2.071569 × − × − × − -1.206677 × − -1.031213 × − × − × − -1.709111 × − -4.563110 × − × − × − -1.026788 × − -1.850203 × − × − × − -4.714943 × − -7.009564 × − × − × − -1.860621 × − -2.514576 × − × − × − -6.710533 × − -8.623420 × − × − × − × − -8.984566 × − × − × − × − -1.151088 × − × − × − × − -9.081402 × − × − × − -7.155788 × − -2.769497 × − × − [16] W.G. Dixon. Dynamics of extended bodies in gen-eral relativity. I. Momentum and angular momentum. Proc.Roy.Soc.Lond. , A314:499–527, 1970.[17] William G Dixon. Dynamics of extended bodies ingeneral relativity. iii. equations of motion.
Philo-sophical Transactions of the Royal Society of LondonA: Mathematical, Physical and Engineering Sciences ,277(1264):59–119, 1974. [18] Robert M. Wald. Gravitational spin interaction.
Phys.Rev. , D6:406–413, 1972.[19] O. Semer´ak. Spinning test particles in a Kerr field. 1.
Mon.Not.Roy.Astron.Soc. , 308:863–875, 1999.[20] K Kyrian and O Semer´ak. Spinning test particles in aKerr field.
Mon.Not.Roy.Astron.Soc. , 382:1922, 2007.[21] D. Bini, G. Gemelli, and R. Ruffini. Spinning test par-ticles in general relativity: Nongeodesic motion in theReissner-Nordstrom space-time.
Phys.Rev. , D61:064013, TABLE V. List of partial amplitudes C ± l n for the same orbit as in Table II. l m n Re { C + lmn } Im { C + lmn } Re { C − lmn } Im { C − lmn } × − -2.987138 × − × − -1.364431 × − × − × − -5.564356 × − × − × − -1.203041 × − × − -1.718247 × − × − × − -4.314063 × − × − × − × − -9.355450 × − × − × − × − -8.650047 × − × − × − × − -5.704990 × − × − × − × − -3.087320 × − × − × − × − -1.461746 × − × − × − × − -6.277339 × − × − × − × − -2.501089 × − × − × − × − -9.388208 × − × − × − × − -3.356380 × − × − × − × − -1.152115 × − × − × − × − -3.820558 × − × − × − × − -1.229874 × − × − × − × − -3.860830 × − × − × − × − -1.182908 × − × − × − × − -3.023775 × − × − × − -1.021033 × − × − -5.960966 × − × − × − -1.071155 × − × − × − × − -2.097613 × − × − × − × − -2.025645 × − × − × − × − -1.413811 × − × − × − × − -8.077376 × − × − × − × − -4.013179 × − × − × − × − -1.795737 × − × − × − × − -7.399921 × − × − × − × − -2.851190 × − × − × − × − -1.038250 × − × − × − × − -3.601133 × − × − × − × − -1.196530 × − × − × − × − × − × − × − -5.945050 × − -9.799602 × − -5.966157 × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − -9.583384 × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − Phys. Rev. , D67:104023, 2003.[23] M. Shibata and Y. Mino. Gravitational waves from aspinning particle plunging into a Kerr black hole.
Phys-ical Review D - Particles, Fields, Gravitation and Cos-mology , 58(6):3–4, 1998.[24] Yasushi Mino. Perturbative approach to an orbital evo-lution around a supermassive black hole.
Phys. Rev. D , 67:084027, Apr 2003.[25] Steve Drasco and Scott A. Hughes. Gravitationalwave snapshots of generic extreme mass ratio inspirals.
Phys.Rev. , D73:024027, 2006.[26] Guillaume Faye, Luc Blanchet, and Alessandra Buo-nanno. Higher-order spin effects in the dynamics ofcompact binaries. I. Equations of motion.
Phys.Rev. ,D74:104033, 2006. [27] Juergen Ehlers and Ekkart Rudolph. Dynamics of ex-tended bodies in general relativity center-of-mass de-scription and quasirigidity. General Relativity and Grav-itation , 8(3):197–217, 1977.[28] Vojtˇech Witzany, Jan Steinhoff, and Georgios Lukes-Gerakopoulos. Hamiltonians and canonical coordinatesfor spinning particles in curved space-time.
Classical andQuantum Gravity , 36(7):075003, April 2019.[29] Gabriel Andres Piovano, Andrea Maselli, and Paolo Pani.Extreme mass ratio inspirals with spinning secondary: Adetailed study of equatorial circular motion.
PhysicalReview D , 102(2), apr 2020.[30] Enno Harms, Georgios Lukes-Gerakopoulos, SebastianoBernuzzi, and Alessandro Nagar. Asymptotic gravita-tional wave fluxes from a spinning particle in circularequatorial orbits around a rotating black hole.
Phys.Rev. , D93(4):044015, 2016.[31] W. Schmidt. Celestial mechanics in Kerr spacetime.
Clas-sical and Quantum Gravity , 19(10):2743–2764, 2002.[32] Ryuichi Fujita and Wataru Hikida. Analytical solutionsof bound timelike geodesic orbits in Kerr spacetime.
Clas-sical and Quantum Gravity , 26(13):135002, July 2009.[33] Eva Hackmann, Claus Laemmerzahl, Yuri N. Obukhov,Dirk Puetzfeld, and Isabell Schaffer. Motion of spinningtest bodies in Kerr spacetime.
Phys.Rev. , D90(6):064035,2014.[34] Georgios Lukes-Gerakopoulos, Jonathan Seyrich, andDaniela Kunst. Investigating spinning test particles: spinsupplementary conditions and the Hamiltonian formal-ism.
Phys.Rev. , D90(10):104019, 2014.[35] Saul A. Teukolsky. Perturbations of a rotating black hole.1. Fundamental equations for gravitational electromag-netic and neutrino field perturbations.
Astrophys. J. ,185:635–647, 1973.[36] Black Hole Perturbation Toolkit. (bhptoolkit.org).[37] Kostas Glampedakis and Daniel Kennefick. Zoom andwhirl: Eccentric equatorial orbits around spinning blackholes and their evolution under gravitational radiationreaction.
Phys.Rev. , D66:044002, 2002.[38] Paul L. Chrzanowski. Vector potential and metric per-turbations of a rotating black hole.
Physical Review D ,11(8):2042–2062, 1975.[39] S.A. Teukolsky and W.H. Press. Perturbations of a ro-tating black hole. III - Interaction of the hole with grav-itational and electromagnetic radiation.
Astrophys.J. ,193:443–461, 1974.[40] Enno Harms, Sebastiano Bernuzzi, and BerndBr¨ugmann. Numerical solution of the 2+1 Teukol-sky equation on a hyperboloidal and horizon penetratingfoliation of Kerr and application to late-time decays.
Class.Quant.Grav. , 30:115013, 2013. [41] Enno Harms, Sebastiano Bernuzzi, Alessandro Nagar,and An Zenginoglu. A new gravitational wave generationalgorithm for particle perturbations of the Kerr space-time.
Class. Quant. Grav. , 31(24):245004, 2014.[42] Alessandro Nagar, Enno Harms, Sebastiano Bernuzzi,and Anil Zenginoglu. The antikick strikes back: recoilvelocities for nearly-extremal binary black hole mergersin the test-mass limit.
Phys.Rev. , D90(12):124086, 2014.[43] Manuela Campanelli, Gaurav Khanna, Pablo Laguna,Jorge Pullin, and Michael P. Ryan. Perturbations ofthe Kerr space-time in horizon penetrating coordinates.
Class.Quant.Grav. , 18:1543–1554, 2001.[44] Pranesh A. Sundararajan, Gaurav Khanna, and Scott A.Hughes. Towards adiabatic waveforms for inspiral intoKerr black holes: I. A new model of the source forthe time domain perturbation equation.
Phys. Rev. ,D76:104005, 2007.[45] Johan Wald´en. On the approximation of singular sourceterms in differential equations.
Numerical Methods forPartial Differential Equations , 15(4):503–520, 1999.[46] Enno Harms, Georgios Lukes-Gerakopoulos, SebastianoBernuzzi, and Alessandro Nagar. Spinning test bodyorbiting around a Schwarzschild black hole: Circulardynamics and gravitational-wave fluxes.
Phys. Rev. ,D94(10):104010, 2016.[47] Georgios Lukes-Gerakopoulos, Enno Harms, SebastianoBernuzzi, and Alessandro Nagar. Spinning test body or-biting around a Kerr black hole: Circular dynamics andgravitational-wave fluxes.
Phys.Rev. D , 96(6):064051,September 2017.[48] Alessandro Nagar, Francesco Messina, Chris Kavanagh,Georgios Lukes-Gerakopoulos, Niels Warburton, Sebas-tiano Bernuzzi, and Enno Harms. Factorization and re-summation: A new paradigm to improve gravitationalwave amplitudes. III. The spinning test-body terms.
Phys. Rev. D , 100:104056, Nov 2019.[49] Misao Sasaki and Hideyuki Tagoshi. Analytic black holeperturbation approach to gravitational radiation.
LivingRev.Rel. , 6:6, 2003.[50] E. W. Leaver. An Analytic representation for the quasinormal modes of Kerr black holes.
Proc. Roy. Soc. Lond. ,A402:285–298, 1985.[51] Viktor Skoup´y and Georgios Lukes-Gerakopoulos. Grav-itational wave templates from Extreme Mass Ratio In-spirals. arXiv e-prints , page arXiv:2101.04533, January2021.[52] Enrico Barausse, Etienne Racine, and Alessandra Buo-nanno. Hamiltonian of a spinning test-particle in curvedspacetime.