"Flux-balance formulae"for extreme mass-ratio inspirals

The"flux-balance formulae"that determine the averaged evolution of energy, azimuthal angular momentum, and Carter constant in terms of the averaged asymptotic gravitational-wave fluxes for inspirals of small bodies into Kerr black holes were first derived about 15 years ago. However, this derivation is restricted to the case that the background Kerr geodesics are non-resonant (i.e., the radial and angular motions are always incommensurate), and excludes the resonant case that can be important for the radiative dynamics of extreme mass-ratio inspirals. We give here a new derivation of the flux formulae based on Hamiltonian dynamics of a self-forced particle motion, which is a valuable tool for analyzing self-force effects on generic (eccentric, inclined) bound orbits in the Kerr spacetime. This Hamiltonian derivation using action-angle variables is much simpler than the previous one, applies to resonant inspirals without any complication, and can be straightforwardly implemented by using analytical/numerical Teukolsky-based flux codes.

one to rewrite the time-averaged radiative self-force in terms of the asymptotic amplitudes of gravitational radiation at infinity and the horizon. The flux formulae then follow from the time-averaged self-forced equations of motion. The time-averaged rates of change ofÊ and L calculated by the flux formulae are matched precisely by the total time-averaged fluxes of energy and angular momentum out to infinity and down to the horizon of a Kerr black hole, as expected [22,40]. Although there is no known gravitational-wave fluxes of Carter constantQ (see, e.g., Ref. [43]), somewhat surprisingly, this self-force derivation can provide the desired flux formulae for the time-averaged rates of change ofQ.
However, the generalization of the flux formulae to the resonant case (ii) did not appear straightforward 3 . The technical difficulties and subtleties trace back to the dependence on initial phases of the orbits in the resonant case [28], and in fact Refs. [35,36,41] take full advantage of non-resonant assumptions to prove the absence of this initial phase dependence. The non-resonant assumptions in the analysis by Mino [41] and Sago et al. [35,36] have been partially removed by Hinderer and Flanagan [16,26] as well as us [45,46], and Flanagan, Hughes, and Ruangsri [47], respectively. The purpose of this paper is to eliminate remaining non-resonant restrictions, and to establish the complete set of flux formulae for radiative inspirals of small bodies into Kerr black holes at O(η), including the case of resonant orbits.

Outline and summary of this paper
Our approach to the flux formulae is based on a Hamiltonian approach in the selfforce theory, originally developed for the conservative self-force dynamics in the Kerr spacetime [46,48,49]. The Hamiltonian method allows us to conveniently formulate the (first-order) self-forced equations of motion as the Hamilton's canonical equations in terms of action-angle variables on the phase space [12,16], which manifestly respect the integrability and tri-periodicity of the bound Kerr geodesic orbits in the test-mass limit. Furthermore, the triplets of {Ê,L,Q} are also promoted to canonical variables on the phase space, which are just invertible functions of action variables. As a consequence, the analysis has a great advantage over the previous spacetime approach of Refs. [35,36,39,41,45], allowing a better control in the resonant case.
We begin our derivation in Sec. 2 with a review of our Hamiltonian formalism of selfforce dynamics in the Kerr geometry [46,48]. Based on the generalized equivalence principle (i.e., the self-forced motion in the background Kerr spacetime is equivalent to the geodesic motion in a certainly perturbed Kerr spacetime) [50,51], the Hamiltonian method describes the motion of a nonspinning point particle by the geodesic Hamiltonian, H ≡ H (0) + H (1) on the 8D phase space spanned by the canonical variables (x µ , u µ ) with the effective metric g µν ≡ g µν , and H (1) (∝ η) is an interaction Hamiltonian that accounts for the (first-order) self-force generated by the "regular" part of the metric perturbation h (R) µν (∝ η) associated with the orbit. We work out the Hamiltonian H in terms of the actionangle variables (w α , J α ) defined in Eqs (20) and (22), and the Hamilton's equation for the actions J α is given byJ where the overdot denotes the derivative with respect to the proper time τ along the orbit (measured in g µν ). Notice thatJ α involves only the interaction Hamiltonian H (1) because the background Hamiltonian is H (0) = H (0) (J). In Sec. 3, we simplify the Hamilton's equation of Eq. (1) by the long-time average. We split the interaction Hamiltonian H (1) into the time-antisymmetric ("radiative") pieces H rad and time-symmetric piece H sym 4 , i.e., H (1) = H rad + H sym , and show that the general (quasi gauge-invariant) expression for the averaged Hamilton's equation ofJ α is 5 where · τ is the average over the proper time τ , and δ/δw α 0 is the total variation with respect to the initial values of angle variables w α 0 . This equation applies to both the non-resonant and resonant orbit, but we will show that the non-resonant orbit yields δ H sym τ /δw α 0 = 0 because H sym τ is guaranteed to be independent of w α 0 in this case. In Sec. 4, we establish a connection between the τ -averaged Hamilton's equation of Eq. (2) and the amplitudes of the gravitational radiation at future null infinity and the horizon. The basic idea is to derive the mode expression of H rad based on the Teukolsky formalism; the (half-retarded-minus-half-advanced) radiative metric perturbation contained in H rad is a global vacuum perturbation, whose asymptotic forms can be expressed in terms of mode solutions to the homogeneous Teukolsky equation. We then show our main results of flux formulae for J α .
In the case of the non-resonant orbit, the flux formulae are expressed as (see Eq. (60)) where · t is the average over an asymptotic time coordinate t, ℓ is the (spheroidal) angular number of gravitational wave modes, ω mkn and p mkn are the functions of the linear combination of the orbital frequencies Ω α in terms of mode integers {m, k, n}, ε α ≡ (−ω mkn , n, k, m), andZ out/dowm ℓmkn are asymptotic amplitudes of gravitational waves modes 'out (down)' to the infinity (horizon).
Here, the on-resonance 'out/down' amplitudes of gravitational wave modesZ out/down ℓmN explicitly depend on the initial value w α 0 of the orbit (see Eqs. (62)). Importantly, to evaluate dJ r /dt t or dJ θ /dt t , we need the additional symmetric contribution of δ H sym τ /δw α 0 as well as dJ /dt t : recall Eq. (2).
In Sec. 5 we derive the flux formulae for the t-averaged rates of change of specific energŷ E, azimuthal angular momentumL and Carter constantQ, built on the flux formulae of Eqs. (3) and (4). In the Hamiltonian formalism, the triplets of {Ê,L,Q} are canonical variables on the phase space, which are invertible functions of J α . After a simple canonical transformation, we obtain dÊ/dt t = − dJ t /dt t , dL/dt t = dJ φ /dt t and (see Eqs. (75) the factors (r 2 + a 2 )P /∆ λ , aP /∆ λ and Υ r , which are all local to the orbit, are defined in Sec. 5.1. Substituting the flux formulae of Eq. (3) into these equations, the end results then reduce to those for dÊ/dt t , dL/dt t and dQ/dt t in the non-resonant case that were first presented by Sago et al. [35,36]. When instead substituting Eqs. (4), we arrive at precisely the same expressions for dÊ/dt t and dL/dt t in the resonant case that were first reported by Grossman, Levin, and Perez-Giz [44], and Flanagan, Hughes and Ruangsri [47]. However, we find that our on-resonance expression for dQ/dt t only partially agrees with that presented in Ref. [47]; the difference is attributed to the fact that Ref. [47] only considers the radiative sector of Eq. (5) obtained from H rad , and discards all contributions to dJ r /dt t coming from H sym by construction (see Eqs. (63)) 6 . We conclude this introduction by briefly discussing some open issues and prospects about our results. First, it seems that the flux formulae of Eqs. (3) and (4) would suggest the "balance laws" for the actions J α of the radiative inspiral of a small body into a Kerr black hole (in a time-averaged sense). At a conceptual level, our above interpretation sounds natural at least for dJ t /dt t and dJ φ /dt t because the corresponding fluxes on their righthand sides in Eqs. (3) and (4) coincide precisely with the fluxes of energy and axial angular momentum that gravitational radiation carries to infinity and down to the horizon [22,40,44,47]. However, the similar interpretation for dJ r /dt t and dJ θ /dt t is subtle. While the forms of Eqs. (3) and (4) would resemble more closely the "balance laws" than that of Eq. (5) for the Carter constant, still, the precise physical meaning of the 'out/down' gravitationalwave fluxes encoded in dJ r /dt t and dJ θ /dt t are not known to us. For example, how does the background Kerr geometry evolve due to the radiative losses of such the "fluxes" of 6 The conservative self-force contribution to dQ/dt t in the resonant case may be explained in a heuristic way, without invoking technical Hamiltonian analysis. Consider a purely conservative selfforced dynamics in the corotating frame with the angular velocity of the periastron advance; the orbital plane does not precess in this frame. Introduce the axes (x,z) in the corotating frame to point the direction connecting the center of the background Kerr black hole and the periastron, and that normal to the orbital plane, respectively. In the non-resonant case, the orbit can ergodically fill up the orbital plane (bounded by two radii r min ≤ r ≤ r max ) that has ax-z plane symmetry, and so there should be no "torque of the conservative self-force" in thex direction (in a certain time-averaged sense). In contrast, the resonant orbit with a different initial value of w α 0 traces out a different "shape" in the orbital plane and breaks the plane symmetry, giving rise to the "torque" in thex direction. Because this torque induces the change in the orbital inclination angle related with the Carter constant Q (see, e.g., Refs. [37,40]), we can see that dQ/dt t on resonance would have the conservative contribution of the self-force (i.e., H sym ).

5/27
J r and J θ ? This question remains an open issue, and perhaps the answer may be provided through the on-going development of the second-order formalism for the perturbed Einstein field equation [52,53].
Second, the flux formulae of Eqs. (3) and (4) can be straightforwardly implemented, making use of the analytical/numerical Teukolsky-based flux codes developed by "B.H.P.C." [23,36,37,[54][55][56][57] with some mild adaptations (or any of the existing Teukolsky platforms such as the ones of Refs. [20,24,40,47,[58][59][60]). Because the flux formulae are built on the radiative sector of the averaged Hamilton's equation of Eq. (2), it is completely equivalent to the effect of standard time-averaged dissipative self-force for generic (resonant) orbits in the Kerr geometry. Furthermore, the flux formulae are (quasi) gauge-invariant characterization of the radiative dynamics for such generic orbits. Therefore, we expect the concrete calculation of the flux formulae will clarify an interesting possibility of the "sustained" resonance in the extreme mass-ratio inspirals (the orbit "stuck" on the resonance) [45,61], help in refining known practical schemes for simulating the radiative extreme-mass-ratio inspirals [6,25,[62][63][64], and provide a (yet another) accurate strong-field benchmark for the extreme mass-ratio regime of a more generic two-body problem in general relativity (see, e.g., Ref. [65]). We shall leave the task of actually producing the numerical data or analytical approximations to our forthcoming publication.

Hamiltonian dynamics in the perturbed Kerr geometry with radiation
To set the stage, we begin with an extension of Hamiltonian formulation of the conservative dynamics in perturbed Kerr geometry [46] to incorporate the full, physical, retarded metric perturbation sourced by a point particle. For the most part we shall import many helpful results from Sec. II of this reference, but the reader should bear in mind that the discussion presented here allows for the gravitation radiation from the particle. Our development below is limited to the self-force theory at first order in the mass ratio for the sub-extremal Kerr spacetime: we shall not be concerned with the second-order perturbation theory [51,52,66,67] and the (perturbed) motion in the extremal Kerr geometry.
Throughout this work the Kerr metric g (0) µν of mass M and spin S ≡ aM (< M 2 ) will be written in terms of Boyer-Lindquist coordinates (t, r, θ, φ). It is given by where The mass of the particle is µ (≪ M ) and the mass ratio is defined by η ≡ µ/M . We will assume that particle's orbits in the limit η → 0 recover generic bound geodesic orbits in Kerr geometry, restricted in r min ≤ r ≤ r max and θ min ≤ θ ≤ π − θ min , respectively.

4D geodesic Hamiltonian and Hamilton's equations
The Hamiltonian formulation in Ref. [46] begins with the observation that the motion of a small particle is the geodesic motion in the smooth vacuum effective metric [50,68]: µν is the background Kerr geometry, and h (R) µν = O(η) is the regular-part of the metric perturbation (≡ the R field) that is defined by subtracting an appropriate singular-part metric perturbation (≡ the S field) h (S) µν (x; γ) from a physical, retarded metric perturbation h + µν (x; γ) generated by the source orbit γ [69]. At first order in the mass ratio, the R field can be always decomposed into the time-antisymmetric ("radiative") field h (rad) µν µν with the advanced metric perturbation h − µν (x; γ), and Ref. [46] discussed only the symmetric field instead of the full R field to define conservative dynamics. Despite this difference in the metric perturbation, the essential picture remains unchanged: the particle moves on a geodesic of the effective metric. The generalization of the Hamiltonian formulation [46] to the full effective metric g µν of Eq. (8) is therefore immediate.
In our framework, the geodesic Hamiltonian in the effective metric is a standard 4D Hamiltonian. We write this as or, equivalently, expand it as where H (0) (x, u) is the unperturbed background Hamiltonian, simply defined by the expression (9) with the substitution g µν → g µν (0) , and is the perturbed interaction Hamiltonian (∝ η). There is no need to display a more explicit form of h (R) µν in this section, and we shall defer it to Sec. 3. The Hamiltonian of Eq. (9) leads to with Hamilton's equations for the canonical position x µ and momentum u µ :ẋ where the overdot stands for the derivative with respect to the proper time τ along the orbit.
Here, it is important to recognize that τ is measured in the effective metric of Eq. (8), not in the background Kerr metric g (0) µν . As a consequence, the canonical variables (x µ , u µ ) have to be normalized according to for any physical orbits (i.e., on-shell solutions of Eqs. (12)).
For these canonical transformations, we first need a set of canonical variables (X α , P α ) such that the canonical momenta P α recover the constants of motion for Kerr geodesics in the test-mass limit η → 0 [12]. Thanks to the symmetry of Kerr geometry associated with 7/27 the Killing vectors t µ and φ µ , and the Killing tensor K µν [70] (see Appendix A for their explicit expressions), the Kerr geodesics admit the three non-trivial constants of motion: the specific energyÊ, azimuthal angular momentumL and Carter constantQ [11]. In the Hamiltonian formalism of Ref. [46],Ê,L andQ are all promoted to canonical momenta P α on the phase space to define where we introduceμ to complete 4D canonical momenta: it is important to distinguish the canonical variableμ from the physical mass of the particle µ.
We next use Eqs. (19) to introduce action variables J α . They are defined by [12,16] where denotes twice the integral over the allowed region of motion described by R(r, P ) ≥ 0 and Θ(cos θ, P ) ≥ 0. By definition, we have and J α (P ) is an invertible function to give P µ = P µ (J), as is shown in Ref. [12]. This allows us to rewrite the generating function of Eq. (16) as W(x, J) ≡ W (x, P (J)) and it then generates the desired canonical transformation between (x µ , u µ ) and action-angle variables (w α , J α ). The results are u µ = (∂W /∂x µ ) J , and With these action-angle variables, Hamilton's canonical equations reaḋ where we use the fact that (∂H (0) (J)/∂w α ) J = 0 in the second equality of the latter equation: recall Eqs. (14). It is important to recognize that (w α , J α ) evaluated along the physical orbit γ must satisfyẇ due to the normalization condition of Eq. (13). While this relation seems to have been unnoticed in many papers, it is now recognized as an algebraic formula to be interpreted as the "first law of binary mechanics" (in the self-force theory) [46,72] 9 , and plays a crucial role in reading out the physical effects of the self-force with the Hamiltonian formalisms [48].

Orbital resonance in Kerr geodesics
The manipulation of the Hamilton's equations of Eq. (23) will require as the input (source) orbit γ for the interaction Hamiltonian H (1) (w, J; γ), and it is approximated by generic bound geodesic orbits in background Kerr spacetime. For future reference, we now take the test-mass limit η → 0 in Eq. (23) and obtain the solutions (w α (τ ), J α (τ )) for the bound Kerr geodesics. The spatial components of angle variables {w r , w θ , w φ } given by Eqs. (22) are 2π-periodic coordinates on the phase space. This observation then implies that the Hamilton's equation 9 "The first law of binary mechanics" is a conjectured variational formula for (point-particle) binary systems that relates the local property of each body to either global conserved charges of the binary system (the Bondi binding energy, etc.) or local conserved quantities along the orbit (the mechanical energyÊ, etc). The first law has been formulated in a number of context of two-body dynamics in general relativity, including the post-Newtonian method, self-force theory and numerical-relativity simulation; the recent development is reviewed in Sec. 3 forẇ α can define the angular frequencies of Kerr geodesics with respect to the proper time τ [12,16,46,72], where we have introduced the (background) "redshift" variable z ≡ (ω t ) −1 in the phase space and the associated fundamental frequencies Ω α measured with respect to the Boyer-Lindquist coordinate time t (i.e., the proper time of a static observer in the asymptotically far region) 10 . With Eq. (25), the solutions to the Hamilton's equations of Eqs. (23) can be expressed as where τ 0 and w α 0 are some initial values, and the constants of motion J α (τ ) = J α given by Eqs. (20).
When the radial and angular motions become commensurable, there exists only one fundamental frequencyΩ for their motions. This is the (r-θ) orbital resonance, which will require a separate treatment in our analysis. For the fundamental frequencies Ω r and Ω θ of Eqs. (25), we define the resonant relation for Kerr geodesics bỹ with a pair of coprime integers {β r , β θ }, and the associated resonant phase by The initial values w α 0 in Eq. (26) will play a very important role in the next two sections. Not all of them, however, are relevant to our discussion. Without loss of generality, we have whether or not the orbit is resonant because the Kerr geometry is stationary and axially symmetric: recall Eqs. (22). While w r 0 and w θ 0 cannot be set to zero making use of the symmetry of the Kerr geometry in general, we are still allowed to have modulo 2π by selecting a suitable value of τ 0 . If the orbit is non-resonant, Eq. (30) further yields w r 0 = 0 = w θ 0 thanks to the fact that the orbital periods 2π/ω r and 2π/ω θ are incommensurable; this is essentially the same argument as firstly produced by Mino [41]. The above dependence of the initial phase of the resonant orbit was stressed in a number of works [28,44,45,47,61], and our observation here agrees with their analysis.

Green's-function-based definition of the interaction Hamiltonian
In this section, we simplify the Hamilton's equations ofJ α in Eqs. (23) by taking a long-time average, defining for various functions f (s) = O(η) along the background orbits parameterized by s (e.g., the coordinate time t or proper time τ ). This is essentially equivalent to working in the leadingorder two-timescale ("adiabatic") approximation of Eqs. (23) (see, e.g., Refs. [15,16,26]) 11 . We begin with the Green's-function-based expression for the interaction Hamiltonian of Eq. (11) to split it into the radiative and symmetric portions, and then summarize some of pertinent properties of the averaged Hamilton's equations based on that split: they are extensively studied in Sec. III and IV of Ref. [46] to which the discussion here refers the reader for details.

Interaction Hamiltonian: the radiative-symmetric decomposition
A key property of the interaction Hamiltonian of Eq. (11) is the functional dependence of the source orbit γ with which the R-field h (R) µν (x; γ) is generated [10]. This aspect becomes especially clear in view of Green's-function-based definition, which is written as where G µν ρσ is the R-part of the Green's function for the linearized Einstein equation defined by the retarded Green's function G µν ρσ + (x; x ′ ) and the S-part of the Green's function G µν ρσ (S) (x; x ′ ) [50]: details are reviewed in Sec. 16 of Ref. [7]. The variables with a prime here correspond to the (source) orbit γ approximated by the generic bound geodesics in Kerr spacetime (at leading order in the two-timescale approximation).
We next turn to decompose H (1) into the anti-symmetric ("radiative") pieces H rad and "symmetric" piece H sym , adopting the split of G µν ρσ (R) into the anti-symmetric and symmetric portions 12 . Following Gal'tsov [42] and Mino [41], we define the anti-symmetric radiative Green's function G µν ρσ (rad) and the (regularized) symmetric Green's function G µν ρσ (sym−S) by where we have introduced the advanced Green's function G µν ρσ . Armed with such definitions, we simply define H rad and H sym by the substitution of G µν ρσ 11 It should be noted that the long-time average of Eq. (43) used for our analysis is not always equivalent to the phase-space average over the angle variables w α in the literature (see, e.g., Ref. [16] for the precise definition). The two averaging procedures are reconciled only when the orbits are non-resonant [44], and we need to be mindful of their difference in the resonant case. 12 We have omitted the superscript '(1)' from H rad and H sym as we shall focus only on the first-order effects.

Averaged Hamilton's equation
The radiative-symmetric split of H (1) is convenient for a number of reasons, and its most important advantage is that the τ -averaged Hamilton's equation J α τ is simplified to for the non-resonant case, and for the resonant case (recall Eq. (27)). Details of the proof of Eqs. (34) and (35) are provided in Sec. III of Ref. [46], but the key ideas behind the proof are that (i) the symmetric interaction Hamiltonian H sym , which is defined by Eq. (32) with the substitution G µν ρσ is symmetric under the exchange of the "field variables" (x µ , u µ ) and "source variables" (x ′µ , u ′ µ ), and that (ii) the τ -averaged symmetric interaction Hamiltonian H sym τ depends only on the actions J α and the initial values of the angle variables w α 0 in Eq. (26). These two observations allow us to write where δ/δw α 0 is the total variation with respect to w α 0 , and we have used an identity 1 2 with x, u). A factor 1/2 in the right hand side of Eq. (36) accounts for δ/δw α 0 acting on both the field and source orbits parameterized by w α (τ ) and w ′α (τ ), respectively, while the partial derivative in the left hand side acts only on the field variable of w α .
Without loss of generality, we can assume Eq. (29), i.e., H sym is independent of w t 0 and whether or not the orbits experience resonance. Similarly, one can further show that because we are always allowed to have w r 0 = 0 = w θ 0 , i.e., H sym τ is independent of w r 0 and w θ 0 : recall the discussion in Sec. 2.3. Equation (34) then follows. For the resonant orbits, there is no known argument to guarantee (∂H sym /∂w r ) J τ = 0 = (∂H sym /∂w θ ) J τ . It is possible, however, to make some progress by taking the linear 12/27 combination β r (∂H sym /∂w r ) J τ + β θ (∂H sym /∂w θ ) J τ with integers β r and β θ of Eq. (27), which characterize resonance. We recall that the τ -derivative of H sym = O(η) implieṡ where we have used Hamilton's equation of Eq. (23). At the linear order in η, the τ -average of its left-hand side of Eq. (40) gives because the orbit in the test-mass limit η → 0 is the bound Kerr geodesics, and H sym = O(η) does not secularly grow (in the class of gauge in which the effective metric of Eq. (8) is well-defined [69]). Inserting this and Eq. (38) into Eq. (40) with Eq. (27), we arrive at With Eq. (36) and the resonant phase w ⊥ 0 of Eq. (28), this relation can be translated to 1 2 Equation (43) confirms that the presence of (∂H sym /∂w r ) J τ and (∂H sym /∂w θ ) J τ is the consequence of the initial-phase dependence of the resonant orbit. Thus, the general expressions for J α τ should be given by Eq. (35) 13 .

Quasi gauge-invariance of J α τ
Before proceeding, we remind the reader that our Hamiltonian formalism assumes on a certain restricted class of gauges in which the R-field perturbation h (R) µν is well defined everywhere around the orbit [69]. Although Eqs. (23) are altered by a gauge transformation, J α τ are quasi gauge-invariant within that class of gauges 14 .
In the context of the Hamiltonian formalism of Ref. [46], the gauge freedom corresponds to the infinitesimal canonical transformation associated with a generating function Ξ ≡ ξ µ u µ with a gauge vector ξ µ = O(η), which describes a standard infinitesimal coordinate transformation, x µ → x µ + ξ µ + O(η 2 ). Assuming that ξ µ = O(η) holds everywhere in the spacetime (to avoid any spurious secular growth in the metric perturbation), the gauge transformation induces J α → J α +δ ξ J α , whereδ ξ J α do not contain any secularly growing terms (see Eq. (4.7) of Ref. [46]). By taking the τ -average of τ -derivative of this relation, we haveδ and this result immediately establishes the quasi gauge-invariance of J α τ . 13 If we further average Eq. (43) with respect to w ⊥ 0 , which is essentially equivalent to the phase space average of H sym over the 2-torus parameterized by w r and w θ , it identically vanishes. This agrees with the conclusion of Ref. [16] (i.e., the phase-space averaged rates change of J α do not have the contribution from the conservative self-force whether or not the orbit experiences resonance).
14 The quantity J α τ is only locally defined along the orbit, and is not gauge invariant in the strict mathematical sense in general relativity. In fact, the gauge transformation here should be (at least) restricted to respect the (tri)periodicity of the orbit obtained from Eq. (23): see Sec. IV of Ref. [46] and Sec. 7.6 of Ref. [4] for more details.

Flux formulae from the averaged Hamilton's equation
In this section, we specialize the discussion to the radiative sector of Eqs. (34) and (35), and translate them into a more practical language of flux formulae. These formulae involve only asymptotic amplitudes of gravitational waves, which are readily computable in welldeveloped Teukolsky's framework of the black-hole perturbation theory (see, e.g., Refs. [21,22] for reviews). We shall leave for future the more difficult task of actual evaluation of the symmetric piece δ H sym τ /δw ⊥ 0 in Eq. (43) for the resonant orbits 15 . Throughout this section, we will use overbar and "c.c." to mean complex conjugate. Our strategy here closely follows a number of techniques developed by Sago et al. [35][36][37], Drasco et al. [39,40,85], Grossman et al. [44] and Flanagan et al. [47], and our presentation is largely patterned after Ref. [47].

Radiative interaction Hamiltonian in terms of Teukolsky mode functions
A natural starting point of our computation may be the Green's-function-based expression for H rad given by Eq. (32). The key notion here is that the tensorial radiative Green's function G µν ρσ (rad) in a particular, traceless "radiation" gauges can be reconstructed from the scalar radiative Green's function for the Teukolsky equation of perturbed Weyl curvature scalars, which was first derived by Gal'tsov [42] (and later corrected in Ref. [39]). The scalar radiative Green's function can be expressed in terms of only its homogeneous solutions to the Teukolsky equations (see Eq. (B12)), which are separated in all of the variables, and the reconstruction procedure is straightforward, relying on the classical method developed long ago by Chrzanowski [76], and Cohen and Kegeles [77,78]. As we have seen, since (the right-hand side of) Eqs. (34) and (35) are quasi gauge invariant, the reconstruction approach provides an extremely efficient route to evaluate them 16 .
The derivation of H rad based on the Teukolsky equations of Eqs. (B6) is provided in Refs. [42] as well as in Appendix A of Ref. [36], and we shall not need the technical details here. Thus, we simply import the final result from Eq. (3.9) of Ref. [36] 17 : where ω is a continuous frequency, ℓm ≡ ∞ ℓ=2 ℓ m=−ℓ with a pair of integers (ℓ, m), p ωm ≡ ω − ma/(2M r + ) (where r + ≡ M + √ M 2 − a 2 ) is the superradiant factor, and 'down' and 'out' modes are defined by specifying the boundary conditions for the Teukolsky functions imposed at the event horizon and infinity (see Eqs. (B9)). Here, we have defined a mode 15 A preliminary work in the toy problem of the scalar-field self-force can be found in Ref. [45]. 16 We should mention that the mass and angular momentum of the background Kerr spacetime are not altered (in the Abott-Deser sense [79]) by adding the metric perturbation associated with the reconstructed G µν ρσ (rad) through the method of Ref. [76][77][78]80]. Refs. [81][82][83] discuss more details of this statement. 17 The retarded and advanced metric perturbations share the same static ω = 0 modes, and they are not contained in H rad . Their contribution is classified here as a part of H sym . 14/27 scalar Φ out/down ωℓm on the phase space (with the spin-weight s = −2) by 18 in terms of a pure mode function (for the Heltz potential) given by (see Eq. (B5a)) where −2 τ † µν is a certain second-order differential operator (B3), s S ωℓm is the spin-weighted spheroidal harmonics, −s R out/down ωℓm is a (spin-flipped) 'out' or 'down' mode of the radial Teukolsky function, and s N ωℓm is a normalization factor. The amplitude of the mode scalar Φ out/down ωℓm , Z out/down ωℓm , is given by integrating the complex conjugate Φ out/down ωℓm The explicit expressions for −2 τ † µν and Φ

Harmonic decomposition of the mode scalars and amplitudes
When inserting Eq. (45) into the τ -averaged Hamilton's equation of Eqs. (34) and (35), we find In this subsection, we shall simplify this expression, first for the non-resonant case and next for the resonant case, respectively. We will consider only the 'out' mode of Eq. (49) below and the 'down' mode is precisely analogous.
The mode scalar Φ out ωℓm defined on the stationary and axially symmetric Kerr background is proportional to e −iωt e imφ , which implies Φ out ωℓm ∝ e −iωw t e imw φ from Eq. (22). Using this result, Φ out ωℓm allows the Fourier expansions in w r and w θ given by Φ out where the sum k,n = +∞ k=−∞ +∞ n=−∞ is over pairs of integers (k, n), and the Fourier coefficients are It is important to recognize that Φ out ωℓm is a function defined on the phase space, not restricted to the orbit accounting for Eq. (24). Inserting Eqs. (26) and (50) into Eq. (48), a straightforward computation returns The choice of spin weights corresponds to that of one of two radiation gauges. For s = −2, the radiative metric perturbation (re)constructed from G µν ρσ (rad) satisfies h Here, we have introduced discretized fundamental frequencies, initial phases, and amplitudes as respectively, replacing ω with ω mkn because of a delta function in Eq. (52). Notice that the amplitude Z out ωℓm depends on the initial phases of the orbit only through the overall phase e iχmkn [40,85].
We may now simplify the τ -averaged derivative of Φ out ωℓm with respect to w α in Eq. (49). Invoking Eqs. (26), (50) and (52), we have where ε α ≡ −ω mkn , n ′ , k ′ , m , and we have used Eqs. (53). The expression of Eq. (54) applies to both non-resonant and resonant orbits. In the non-resonant case, the τ -average of Eq. (54) forces k ′ = k and n ′ = n in the sum, and we just arrive at and the initial phase w α 0 of the orbit does not appear here. In the resonant case, however, all pairs (k, n) satisfy kΩ θ + nΩ r = NΩ with an integer N ≡ kβ θ + nβ r owing to the resonant condition of Eq. (27). In terms of the angle variables w r and w θ , this implies withω ≡Ω/z and Eq. (28). Then, the τ -average of Eq. (54) enforces N ′ = N (where N ′ ≡ k ′ β θ + n ′ β r ) in the sum, and a straightforward computation, making use of the second equality of Eq.(57), returns 16/27 where ε α = (−ω mN , n ′ , k ′ , m) now, and respectively. The notation (k, n) N means the summation over all pairs of (k, n) which satisfy the relation of N = kβ θ + nβ r . Notice that the cross term between pairs (k, n) N and (k ′ , n ′ ) N remains even if k = k ′ and n = n ′ . It is important to recognize that, unlike Eq. (56), the expression of Eq. (58) can generically depend on the initial phase through the resonant phase w ⊥ 0 . Indeed, the resonant phase term (k − k ′ )β θ w ⊥ 0 of Eq. (58) vanishes only after the average over w ⊥ 0 . These key observations agree with the results found in Refs. [44,47].

Flux formulae
Our final task is to assemble the results in previous subsections, and derive flux formulae from the τ -averaged Hamilton's equations of Eq. (49). For non-resonant case, Eq. (56) can now be substituted into the right-hand side of Eq. (49). We simply arrive at where we have used the relation about the various long-time average (see, e.g., Sec. 9 of Ref. [39]) with the redshift variable z (recall Eq. (25)). Equation (60) is the final form of the flux formulae for non-resonant orbits. For the resonant case, recalling Eq. (58), it is convenient to introduce the initial-phase dependent amplitudes of 'out'-mode defined by [44,47] and those of 'down'-modes as well. Then, Eq. (58) can be substituted into the right-hand side of Eq. (49). By using Eqs. (35) and (43), a simple computation with the relation of Eq. (61) gives (now recovering the symmetric contribution to Eq. (49) for completeness) 17/27 which are the final forms of the flux formulae for resonant orbits. In general, the flux formulae of dJ r /dt t and dJ θ /dt t in Eqs. (63) are invalid unless one includes δ H sym τ /δw ⊥ 0 . Nevertheless, we see that the special linear combination is valid because it does not involve any contribution from δ H sym τ /δw ⊥ 0 at all, thanks to Eqs. (42) and N = kβ θ + nβ r . This expression is also useful since it is expressed in terms of Z out,down ℓmN only, in the same manner as Eqs. (63a) and (63b). We expect that a most advanced self-force code (such as the one of Ref. [64]) will soon be able to test this relation directly.

Flux formulae for the energy, angular momentum and Carter constant
In this section, we calculate the canonical transformation between the evolution of the specific energy, azimuthal angular momentum and Carter constantṖ α (recall Eq. (14)) and that of the action variablesJ α using Eq. (21), and produce the flux formulae for dP α /dt t in terms of dJ α /dt t obtained in Sec. 4. To prepare the way for the discussion, we introduce the (Carter-)Mino time [11,41] λ related to the proper time τ of Eq. (13) by and the associated long-time averaging f (λ) λ for various functions f (λ): recall Eq. (31).

Identities for Kerr geodesics
We establish here a number of identities satisfied by the partial derivatives of the action variables J α with respect to the canonical momenta P α in the test-mass limit η → 0. The notation for the expression of Kerr geodesics in terms of λ is adopted from Drasco and Hughes [85], and Fujita and Hikida [86]. We begin with computing the partial derivatives (∂J r /∂Q) and (∂J θ /∂Ĉ). From the definitions of actions of Eqs. (20), we immediately obtain (see, e.g., Eqs.(3), (5) and (7) of Ref. [86]) where r min/max (r min ≤ r max ) are two largest roots of R(r, P ) = 0, cos θ min > 0 is the smallest positive root of Θ(cos θ, P ) = 0 (recall Eqs. (17)), and are the angular frequencies of Kerr geodesics defined with respect to the Mino time λ. Note that we have (∂r min /∂P α ) √ R = 0 = (∂r max /∂P α ) √ R because R(r, P ) = 0 at r = r min/max , (and similarly for (∂ cos θ min /∂P α ) √ Θ = 0). Equations (66) are clearly separated in r and θ. This comes from the separation property of the generating function W (x, P ) in Eq. (16), 18/27 and it reflects the fact that the Carter constantQ (orĈ) is the separation constant for the Hamilton-Jacobi equation of H (0) in Eq. (10) [11,12].

Averaged evolution of the energy, angular momentum and Carter constant
First, we derive the expressions forṖ α in terms ofJ α using the canonical transformation of Eq. (21). The evolution of the specific energy and angular momentum,Ė = −J t andL =J φ trivially follow from their definitions in Eqs. (20). The expressions for the evolution of the specific Carter constants,Q andĊ are easily produced from the proper-time derivative of J r = J r (H (0) ,Ê,L,Q) and J θ = J θ (H (0) ,Ê,L,Ĉ). After some simple algebra, making use of the identities Eqs. (66), (70) and (71), we arrive at We will now simplifyṖ α averaged over τ to derive the associated flux formulae. Again, the averaged rate of change ofÊ andL are trivially given by making use of the relation of Eq. (61). To simplify Eqs. (72), recall that Eqs. (9) and (13) imply thatḢ =Ḣ (0) +Ḣ (1) = 0, which means (recall the computation of Eq. (41)) for the τ -averaging. With this relation, the τ -average of Eqs. (72) then read where we have used Eqs. (61) and (73). We can clearly see the r-θ split in the averaged rate of change of the Carter constants: the first version of the expression dQ/dt t in Eq. (75a) is described by only the 'r-components' (i.e., r, Υ r , dJ r /dt t etc.), and the second version of the expression dĈ/dt t in Eq. (75b) involves only the 'θ-components' (i.e., θ, Υ θ , dJ θ /dt t etc.). However, they are not independent formulae because of the relation of Eq. (18), and indeed, Eq. (75a) is equivalent to Eq. (75b). This statement is easily understood, making use of the relation obtained from "the first-law of binary mechanics" [46,87,88]. Importing, for example, from Sec. 5 of Ref. [72] and keeping in mind that the fundamental frequency Ω α of Eq. (25) is related to the Mino-time frequency Υ α by the relation Ω α = Υ α /Γ, we have 20 with (see, e.g., Eqs.(7) of Ref. [86]) Equations (75) and (76) can now be substituted into the long-time average of the t-derivative of Eq. (18) given by dĈ/dt t = dQ/dt t − 2(aÊ −L){a dÊ/dt t − dL/dt t } which easily reveals the equivalence between Eqs. (75a) and (75b). Equations (73) and (75) are the final form of the flux formulae for the energy, azimuthal angular momentum and Carter constants (in terms of those for J α ). When the flux formulae of Eq. (60) for the non-resonant orbits are inserted into Eqs. (73) and (75), we have the same results obtained by Sago et al. [36] and displayed in their Eqs. (3.13), (3.15), (3.24) and (3.26), respectively. Similarly, for the resonant orbits substitution of Eqs. (63) in Eqs. (73) and (75) gives the same results obtained by Flanagan et al. [47] and displayed in their Eqs. (3.35) -(3.40), respectively. We note that, however, the flux formulae of dQ/dt t in Eq. (3.40) of Ref. [47] must be supplemented by the additional contribution from the symmetric interaction Hamiltonian H sym . This is clear especially in view of Eqs. (75), i.e., in the resonant case dJ r /dt t and dJ θ /dt t in Eqs. (75) generally require the flux formulae of Eq. (63), involving δ H sym τ /δw ⊥ 0 .

A. Killing vectors and tensors for Kerr geometry
In this Appendix, we collect a few key results of Killing vectors and tensors of the Kerr spacetime that play a central role when constructing constants of motion for Kerr geodesic orbits: recall Sec. 2.2. The stationary and axisymmetric Kerr geometry admits the two Killing vectors Besides, it is now well-known that there is another "hidden symmetry" of the Kerr geometry associated with the (rank-2 irreducible) Killing tensor K µν = K (µν) [11,70,89] that satisfies the Killing's equation ∇ µν , and parenthesis embracing indices are the total symmetrization of a given tensor. Explicit expressions for K µν may need a standard Kinnersley null tetrad (recall Eq. (7)), and or, the associated basis 1-forms and m α = 1 √ 2(r + ia cos θ) −ia sin θ, 0, Σ, i(r 2 + a 2 ) sin θ (A5) that satisfy ℓ µ n µ = −1 and m µ m µ = 1 (the overbar denotes complex conjugate). Making use of these basis 1-forms, we then write the Killing tensor as [39,90] K µν ≡ 2a 2 cos 2 θ ℓ (µ n ν) + 2r 2 m (µ m ν) .
(A6) 21/27 The Kerr metric of Eq. (6) can also be written in terms of the basis 1-forms {ℓ µ , n µ , m µ , m µ } as With Eqs. (7), this can be substituted into the right hand side of Eq. (A6) to establish the "duality" of the Killing tensor:
The master variables in the Teukolsky formalism are essentially Weyl scalars ψ 0 and ψ 4 defined by (recall the Kinnersley null tetrad of Eq. (A2)) where C αβγδ is the (perturbed) Weyl tensor, and ρ ≡ (r − ia cos θ) −1 . The master variables s Ψ satisfy the Teukolsky equation where s O and s τ αβ are differential operators for the spin-weight s, and T αβ is the energymomentum tensor of the matter source (point-particle source etc). In the bulk of our paper, we may need the explicit expression for the adjoint of −2 τ αβ which is given by (see, e.g., Eq. (A.28) of Ref. [36]) 21 , where the 'plus (+)'-operators are with integers n and s. The explicit expressions for other operators are displayed in, e.g., Appendix A of Ref. [36]. 21 Suppose a linear differential operator M acting on an n-index tensor T , and taking it to another k-index tensor M • T . The adjoint of M is defined by M † so that the relation ( is satisfied for any pair of such two linear operators M 1 and M 2 [94]. Here, we use a dagger ( †) to denote the adjoint, and a 'plus (+)' to imply the transformation (ω, m) → (−ω, −m). Notice that these symbols would differ from common choices; e.g., † and + operators correspond to a 'star ( * )' and 'dagger ( †)' in Ref. [36], respectively.

22/27
The Teukolsky equation of Eq. (B2) admits a full separation of variables in the frequency domain. We may write the solution s Ψ and the source term s T ≡ s τ αβ T αβ as where ω is a continuous frequency, ℓ and m are integers, and we have introduced the spinweighted spheroidal harmonics s S ωℓm ≡ (1/ √ 2π) s Θ ωℓm (θ)e imφ . Substituting Eqs. (B5) into Eq. (B2), we obtain the angular and radial Teukolsky equations with potentials and K ≡ ω(r 2 + a 2 ) − am.
The differential equation of Eq. (B6a) defines the (polar part of) spin-weighted spheroidal harmonics s Θ ωℓm (θ) normalized as π 0 dθ sin θ| s Θ ωℓm | 2 = 1, and the associated eigen values s λ ωℓm in Eqs. (B7). At the same time, the homogeneous solution of Eq. (B6b) defines the four independent "mode" functions, depending on their boundary conditions at the infinity and horizon. In keeping with common nomenclature, they are defined by (see, e.g., Ref. [42]) for 'up' and 'in' modes, and for 'out' and 'down' modes, where p mω ≡ ω − ma/(2M r + ) and r * is the tortoise coordinate that satisfies dr * /dr = (r 2 + a 2 )/∆. The complex-valued coefficients s B inc ωℓm , s B ref ωℓm and s B trans ωℓm ( s C inc ωℓm , s C ref ωℓm and s C trans ωℓm ) are, respectively, incidence, reflection and transmission coefficients of the 'in'-mode ('up'-mode) solutions. 23/27 For given boundary conditions, the Green's function of the Teukolsky equation of Eq. (B2) is defined by the solution of the differential equation where δ (4) (x − x ′ ) is a 4D coordinate delta function, and its explicit expressions can be constructed in terms of the mode functions of Eqs. (B8) and (B9). For example, the retarded Green's function that satisfies the retarded boundary condition G + (x, x ′ ) = 0 for t < t ′ is given by (see, e.g., Eq.(A.38) of Ref. [36]) where we have introduced a step function H(x) ≡ x −∞ δ(y)dy. The advanced Green's function G − (x, x ′ ) that satisfies the advanced boundary condition G − (x, x ′ ) = 0 for t > t ′ is similarly obtained by Eq. (B11) × s A s R down ωℓm (r) −s R down ωℓm (r ′ ) + ω p ωm s B s R out ωℓm (r) −s R out ωℓm (r ′ ) .(B12) We do not need the explicit expressions for the normalization factors s A and s B here, but they can be straightforwardly computed from the results in Appendix A of Ref. [36]. Notice that, unlike Eq. (B11), the radiative Green's function of Eq. (B12) contains no step function. The step functions in s G + (x, x ′ ) and s G − (x, x ′ ) are exactly cancel each other out, and s G (rad) (x, x ′ ) indeed satisfies the homogeneous Teukolsky equation, s O s G (rad) = 0. This is the key property of s G (rad) (x, x ′ ) that allows us to construct the radiative interaction Hamiltonian in the simple "homogeneous" form of Eq. (45), leading to the flux formulae displayed in Eqs. (60) and (63); for the full details to obtain Eq. (45) from Eq. (B12), once again, we refer readers to Appendix A of Ref. [36].

C. More on evolution of constants of motion
Unlike the flux formulae of Eqs. (60) and (63), it is not so difficult to express the evolution of constants of motionṖ α (see Eq. (14)) directly in terms of the local self-force f µ (∼ O(η)) [95,96]. It is given in our notation by using Refs. [17,36] 22 To be precise, the proper time of those literature isτ normalized with respect to the background Kerr metric g (0) µν , which differs from ours τ of Eq. (13). However, this difference is negligible here because dτ = dτ (1 + H (1) ) [97], and the self-force is already f α ∼ O(η).

24/27
and this is a so-called "forcing terms" in Refs. [16,61,62,64]. The objective of this appendix is to examine how this expression can be derived in the Hamiltonian formulation.
We write the Hamilton's equations in the "mixed" canonical variables of (x µ , P α ) ≡ (x µ (X, P ), P α ). From the Hamilton's equation for P α , we quickly obtaiṅ where ∂P α /∂X β = 0, and the last equality follows from (∂H (0) (P )/∂X α ) P = 0 as well as an identity with the standard Poisson bracket {x µ , P α } 23 . This gives a derivation ofṖ α expressed in the form of Eq. (C1) in the language of Hamiltonian formalism. In particular, the comparison of Eq. (C1) with Eq. (C2) reveals an interesting relation,