Evolution of the Carter constant for a resonant inspiral into a Kerr black hole: I. The scalar case

We discuss the inspiral of a small body around a Kerr black hole. When the time scale of the radiation reaction is sufficiently longer than its orbital period, the leading order orbital evolution is described only by the knowledge of the averaged evolution of the constants of motion, i.e., the energy, azimuthal angular momentum and the Carter constant. Although there is no conserved current composed of the perturbation field corresponding to the Carter constant, it has been shown that the averaged rate of change of the Carter constant can be given by a simple formula, when there exists a simultaneous turning point of the radial and polar oscillations. However, an inspiralling orbit may cross a"resonance"point, where the frequencies of the radial and polar orbital oscillations are in a rational ratio. At the resonance point, one cannot find a simultaneous turning point in general. Hence, even for the averaged rate of change of the Carter constant, a direct computation of the 'self-force", which is still challenging especially in the case of the Kerr background, seems to be necessary. In this paper, we present a method of computing the averaged rate of change of the Carter constant in a relatively simple manner at the"resonance"point.


I. INTRODUCTION
The study of the two-body problem in general relativity has attracted much attention in recent years in order to obtain a better understanding of the gravitational wave emission from a binary system. If one body is much less massive than the other, i.e., when we consider extreme mass ratio inspirals of a compact object into a super massive black hole, the compact object can be treated as a point particle and the two-body problem is handled within the self-force program (see Refs. [1,2] and references therein).
When a nonspinning particle with mass µ is moving along a bound orbit around a Kerr black hole with mass M (≫ µ) and the spin parameter a, the long-term evolution of the orbit due to the radiation reaction of the gravitational wave emission is dictated by the long-time averaged rates of change of three constants of motion: the energy E, azimuthal angular momentum L and Carter constant Q [3] since the characteristic time scale of the secular orbital evolution is sufficiently longer than the orbital period. The orbital phase errors originating from the other effects scale as O((M/µ) 0 ). This particular approximation of the orbital evolution is referred to as the adiabatic approximation [4][5][6][7][8].
On average, the energy and azimuthal angular momentum losses of the particle are balanced with those carried by gravitational waves to the infinity or into the Kerr horizon due to their global conservation laws. This is nothing but the well established balance argument. Although we do not have any balance argument for the Carter constant, based on Mino's argument [9], Sago et al. [10,11] (see also Drasco et al. [12] for the scalar case) have succeeded in deriving a formula for the averaged rate of change of the Carter constant dQ/dt similar to that obtained by the balance argument for the energy and azimuthal angular momentum.
However, the above-mentioned work assumes that the particle's radial and polar orbital frequencies, Ω r and Ω θ , are in an irrational ratio. When Ω r and Ω θ cross a resonance point, at which Ω r /Ω θ = j r /j θ holds with coprime integers j r and j θ , the above formula is not correct, and the phase error due to ignorance of the correct formula can be as large as O((M/µ) 1/2 ). This was first pointed out by Mino [4] (see also Ref. [5]). Recently, Flanagan and Hinderer [13] investigated this issue quantitatively within the post-Newtonian approximation. For the inspiral in the extreme mass ratio regime M ≫ µ, this phase error is much larger than that from the first-order conservative self-force or the long-time averaged second-order dissipative self-force, which scales as O((M/µ) 0 ) [5][6][7][14][15][16]. For this reason, it is very desirable to derive a formula for dQ/dt in the adiabatic approximation for resonance orbits.
The main differences from the off-resonance case in the derivation of dQ/dt are the following two points. First, the orbit is no longer characterized by the three constants of motion alone. Instead, as was pointed out by Flanagan et al. [17], the orbit depends on the relative offset phase specified by ∆λ, the separation between the instances when the rand θ-oscillations reach their minima. The offset phase also slowly evolves compared to the orbital period near the resonance point. Thus, the resonance orbit must be characterized by four parameters: {E, L z , Q, ∆λ}. Second, we cannot use Mino's argument [4], which proves that dQ/dt can be evaluated by using the radiative field defined by "the half-retarded minus half-advanced field", which is by definition regular at the location of the particle. The key ingredient of Mino's argument is an invariance of the geodesic under the transformation (t, r, θ, φ) → (−t, r, θ, −φ) that is assured if and only if the rand θ-components of the body's four velocity vanish simultaneously at a certain time, i.e., when we can choose ∆λ = 0. Clearly, this condition is not satisfied for resonance geodesics.
Our final goal is to obtain the evolution of resonant inspirals with a phase accuracy better than O((M/µ) 0 ). Toward this ambitious goal, we here begin with deriving a practical formula of dQ/dt for the resonant inspirals. In the resonance case, we find that dQ/dt acquires contributions from both the radiative and symmetric fields, the latter defined by "the half-retarded plus half-advanced field". In general, the symmetric field diverges at the particle's location, thus requires some regularization. In a later section, we will propose a new regularization scheme for the symmetric part of dQ/dt . The key point of our regularization method is that it is compatible with the Teukolsky formalism [18]. This makes the calculation much simpler than the self-force calculation, which has already proved challenging numerically for a particle orbiting around a Schwarzschild black hole (see Ref. [1] for recent progress).
In this paper, to avoid complications, we focus on a point particle with mass µ coupled to a massless scalar field with the charge q. Although in principle there seems to be no essential obstacle to extending our scalar formalism to the gravitational case, the latter is more involved since it requires the reconstruction of metric perturbations, which is a little complicated in the Kerr case [19,20].
The remainder of this manuscript is organized as follows. In Sec. II, we begin with describing the model of a point scalar particle coupled to a massless scalar field. We review the generic bounded geodesic motion around a Kerr black hole, particularly focusing on the resonance case. The retarded solution of the scalar field equation expressed by the mode decomposition is also summarized there. We derive a general formula for the long-time averaged rate of change of the Carter constant in Sec. III. We separate the formula into two parts, the radiative and symmetric parts that originate from the radiative and symmetric fields, respectively. The result for the radiative part is shown in Sec. III B, and that for the symmetric part is presented in Sec. III C, where we also discuss the mode sum regularization of the symmetric part since it diverges at the particle's location. Finally, we summarize this paper in Sec. IV. To clarify the impact of our work, we also show how the long-term evolution of the resonance orbit is described in the adiabatic approximation. In Appendix A, we discuss a possibility that the resonance effect on the phase error might be larger than O((µ/M ) −1/2 ). Some detailed derivations of formulas related to the regularized symmetric part are summarized in Appendix B.
Throughout this paper, we use geometrical units G = c = 1 and the sign convention of the metric is (−, +, +, +). We basically follow the notation used by Drasco et al. [12].

II. PRELIMINARIES
We consider a point scalar particle with the bare rest mass µ = 1 1 and scalar charge q, moving in the Kerr spacetime. The Kerr metric in the Boyer-Lindquist coordinates is given by where M and a are the mass and the Kerr parameter, respectively, and Σ := r 2 + a 2 cos 2 θ, There are two Killing vectors associated with the stationarity and axisymmetry of the Kerr spacetime: The Kerr spacetime also has an irreducible Killing tensor given by with the Kinnersley null tetrad defined by ℓ µ := (r 2 + a 2 , ∆, 0, a)/∆, n µ := (r 2 + a 2 , −∆, 0, a)/(2Σ) and m µ := (ia sin θ, 0, 1, i/ sin θ)/( √ 2(r + ia cos θ)), where ℓ µ n µ = −1 and m µm µ = 1 while all the other inner products vanish.

A. Bounded geodesic motion around a Kerr black hole in resonance
Neglecting the effect of self-interaction, a particle can be regarded as a test particle that moves along a geodesic on the background spacetime, z µ (τ ) := (t(τ ), r(τ ), θ(τ ), φ(τ )) parametrized by the proper time τ . There are three integrals of motion from the symmetries of the Kerr spacetime: the energy E, azimuthal angular momentum L and Carter constant Q [3], respectively. These are expressed as where u α := dz α /dτ is the particle's four velocity (note that E and L are O(µ), and Q is O(µ 2 ) if we recover µ).
For the Kerr geodesics, it is convenient to parametrize the particle's orbit as z(λ) := (t(λ), r(λ), θ(λ), φ(λ)) with the "Mino time" λ, which is related to τ as dλ := dτ /Σ [4]. Using λ, the equations of motion are written as where Here it should be noted that the equations of r and θ are completely decoupled, and those of t and φ are divided into r and θ -dependent parts. First, the rand θ -oscillations are independently periodic for a bound orbit: from Eqs (5). Here, n and k are integers and the periods Λ r and Λ θ with respect to the Mino time λ are given by with the turning points r min , r max and θ min (< π/2). For later convenience, we introducer andθ which denote the fiducial solutions of Eqs. (5) such thatr(nΛ r ) = r min andθ(kΛ r ) = θ min .
The precise meaning of the resonance is that the orbital frequencies where j r and j θ are coprime integers. We define the resonance frequency, Υ, and the corresponding period, Λ, as In the resonance case, there is a difference between the Mino times of reaching the minima of the rand θ -oscillations in general. We call this difference the "offset phase" and denote it by ∆λ. We note that ∆λ can be shifted by depending on which minima we focus on. Therefore, we choose ∆λ such that |∆λ| ≤ Λ ′ /2. Setting the origin of the λ coordinate to satisfy θ(λ) =θ(λ), the orbit is written as without loss of generality. Second, the motion in the t-direction is also separated into rand θ-dependent parts. We introduce the averaged values of V tr (r) and V tθ (θ) defined by In integrating the t-component of the equation of motion, the choice of the origin of the t-coordinate can be changed freely because of the time-translation symmetry of the Kerr spacetime. Then, without loss of generality, t(λ) is given by [12] t(λ) = Υ t λ + ∆t r (λ − ∆λ) + ∆t θ (λ), where Υ t := V tr + V tθ , and The function φ(λ) is deduced in the same way as t(λ). Again considering the axial symmetry of the Kerr spacetime, we have where we define Υ φ := V φr + V φθ and with From the above discussion, the orbit is expressed as with the two-variable functionz where the oscillating part of the orbit is denoted by

B. The retarded solution of the scalar field equation
The equation of the scalar field induced by a charged particle is given by where ∇ µ and g are the covariant differentiation and the determinant for the background Kerr metric, g µν , in Eq. (1), respectively (see Ref. [12] for further details). Assuming that there is no incoming scalar wave, we take the retarded solution, Φ (ret) . To obtain the solution, we here construct the retarded Green function in terms of mode functions. From the fact that the scalar field equation is separable in the Kerr spacetime [18,21], we can write the mode functions in the form of where r * is the tortoise coordinate defined by dr * := ((r 2 + a 2 )/∆)dr, and S ωlm (θ, φ) := Θ ωlm (θ)e imφ is the spheroidal harmonics normalized as The superscript ♭ represents one of the four distinct boundary conditions: "up", "down", "in" and "out". The corresponding radial functions are defined by [12,22] respectively. Here p mω := ω − ma/(2M r + ), r + := M + √ M 2 − a 2 , and the superscript * denotes the complex conjugate. The coefficients α ωlm and β ωlm are normalization constants. τ ωlm and σ ωlm are complex transmission and reflection coefficients, respectively. Both the complex coefficients µ ωlm and ν ωlm are also defined by Eqs. (23). In this work, we choose the normalizations of the "in" and "up" solutions as so that the invariance of the radial equation under the transformation (ω, m) → (−ω, −m) implies the identity From this identity (25), we find It should also be noted that we choose the phase of the normalization of the angular function Θ ωlm (θ) as The retarded Green function for the scalar field equation (20) is expressed in the factorized form of the mode functions [12]: In showing the second equality of Eq. (28), we used the relations between the radial functions u ♭ ωlm (r * ) in Eq. (23). This Green function leads the retarded solution as For the later analysis, it is helpful to decompose the retarded field into a linear combination of the radiative (anti-symmetric) and symmetric fields to be where with Here the advanced Green function, G (adv) (x, x ′ ), is derived from the retarded Green function as The radiative Green function can be rewritten as [12] and the Heaviside functions are eliminated. On the other hand, we rewrite the symmetric Green function in the following way. Writing G (adv) (x, x ′ ) in terms of the mode functions in a manner similar to Eq. (28), replacing (ω, m) with (−ω, −m), and rearranging it with Eqs. (24), (25), (26), and (27), we find Substituting Eqs. (28) and (36) into Eq. (33), we find that the symmetric Green function is rewritten as where g (sym) and ℑ[X] represents the imaginary part of X.

III. ADIABATIC EVOLUTION OF THE CARTER CONSTANT AT RESONANCE
We denote the three constants of motion as I i := {E, L, Q}. The time averaged rates of change of I i can be re-expressed as with the definition of the λ average, 2 One key property of the adiabatic approximation developed in Ref. [4] for off-resonance orbits is that dI i /dλ λ can be evaluated by using the radiative solution of the field equation (20) as where f α is a differential operator defined by [23][24][25] Since the radiative field has no divergent behavior at the location of the particle, this relation gives dI i /dλ λ without any regularization procedure. The final expressions for dI i /dλ λ are given in terms of the asymptotic amplitudes of the scalar field at infinity and on the horizon [12]. The derivation of Eq.(41) relies on the assumption that there exists a simultaneous turning point of the rand θ -oscillations of the particle's motion. In the case of a resonance orbit, however, the frequencies of the rand θ -oscillations are in a rational ratio, and hence one cannot find a simultaneous turning point in general. This means that the orbit depends on the non-vanishing offset phase, ∆λ. Therefore, the expressions for dI i /dt t should be derived without using the relation (41).
Incidentally, as for dE/dt t and dL/dt t , the balance argument can be used to relate them to the fluxes of the energy and angular momentum to infinity or into the black hole horizon. Therefore, one can conclude that the expressions for dE/dt t and dL/dt t in the resonance case are the same as those in the non-resonance case [17]. On the other hand, there is no known flux composed of the perturbation field corresponding to the Carter constant, and hence we cannot use the balance argument to derive the averaged rate of change, dQ/dt t . We need to reformulate the calculation of dQ/dt t without using the relation (41). In this section, we derive a formula for computing dQ/dt t in the resonance case, by starting from the basic formula of the self-force acting on a particle.

A. Formal reduction of dQ/dt t for a resonance orbit
In contrast to dE/dt t and dL/dt t , we need to directly evaluate the time derivative of the Carter constant defined in Eq. (4). This is given by the self-force acting on the particle as with the aid of the Killing tensor equation ∇ (γ K αβ) = 0. Here, a α is the four acceleration due to the scalar self-force with the retarded scalar field given by [26] where Φ (R) is the so-called R-part of the retarded solution, which satisfies the homogeneous equation of Eq. (20). Φ (R) is schematically expressed as where Φ (S) (x) is the S-part of the retarded field [26], which satisfies the same equation as the retarded (or symmetric) solution.
It should be noted that, the difference between Φ (sym) and Φ (S) is regular at the location of the particle while each of them is singular. Now plugging Eq. (44) into Eq. (43), the time derivative of the Carter constant is rewritten by the R-part of the retarded field Φ (R) (x): Since the first term is a total derivative up to the leading order in q, it does not contribute after taking a long-time average. Hence, we omit the first term in Eq. (46). Substituting the explicit expression for K αβ into Eq. (46), and dropping the first term, we obtain whereÔ x is a differential operator defined bŷ Following the decomposition of the retarded field in Eq. (45), we also decompose dQ/dt t into two parts as dQ dt where G (S) (x, x ′ ) is the S-part of the retarded Green function (see Eq. (81) below). We reduce Eqs. (50) and (51) into more tractable expressions. For this purpose, we consider a function of two variables, we find Here we insert the δ-function, δ(λ r − λ θ + ∆λ), on the second line to deal with the rand θ -oscillations separately, keeping the information on the offset phase. In the last equality, we perform integration by parts with respect to λ r and replace dδ(λ r − λ θ + ∆λ)/dλ r with dδ(λ r − λ θ + ∆λ)/d(∆λ). The differentiation d/d(∆λ) can be extracted to the outside of the integral because the integrand depends on ∆λ only through this δ-function.
As for dQ/dt (rad) t , we arrive at the following expression with the aid of Eq.(53): We call attention to the fact that the ∆λ-differentiation in the second term operates only through z(λ), not through z(λ ′ ).
For dQ/dt we obtain a similar formula to the radiative part as dQ dt where Using the following properties of G (sym−S) (x, x ′ ) (see Appendix B for details), we find that the first term in Eq. (55) vanishes, and then obtain dQ dt In Eq. (59), the differentiation with respect to ∆λ operates on both z(λ) and z(λ ′ ), while, in the second term of Eq. (55), it operates only on z(λ). Since the part in the square brackets in Eq. (59) is symmetric under the exchange of z(λ) and z(λ ′ ), however, the contribution from the ∆λ-differentiation through z(λ ′ ) is the same as that through z(λ). This is the reason why the factor 2 in the expression of Eq. (55) disappears in Eq. (59). The detailed derivation is also presented in Appendix B. The result that the tand φ-derivative parts in Eq. (55) have no contribution, can be expected from the fact that dE/dt t and dL/dt t have no contribution from the symmetric part of the Green function, which is shown by the balance argument of the corresponding conserved currents. We emphasize that, from the above expression, dQ/dt (sym−S) t will not vanish in general. This means that the secular change of the Carter constant would also be induced by the contribution from the symmetric part for the resonance orbit.
Our first task is to rewrite Eq. (54) in terms of the mode functions given in Eq. (21). With the aid of the radiative Green function (35) expressed as a sum of mode functions, we have dQ dt where, for brevity, we have introducedπ ♭ ωlm (x) := Σ(x)π ♭ ωlm (x), and also replaced the partial differentiations with respect to t and φ with −iω and im, respectively. We can writeπ ♭ ωlm [z(λ)] as with a function of two time variables, λ r and λ θ : where ∆z(λ r , λ θ ) is the oscillating part of the orbit defined in Eq. (19). Since the function J ♭ ωlm (λ r , λ θ ) is biperiodic with the frequencies Υ r and Υ θ for λ r and λ θ , respectively, it can be expanded in the Fourier series as where we introduce ω mnrn θ := Υ −1 t (mΥ φ + n r Υ r + n θ Υ θ ) = mΩ φ + (j r n r + j θ n θ )Ω, (69) and withN := {(n r , n θ )|j r n r + j θ n θ = N } .
Here it should be noted that N denotes an integer, and that we rearranged the order of summation over (n r , n θ ) in the second equality in Eq. (68).
Owing to the δ-function in Eq. (68), we can replace the Fourier integral in Eq. (64) with the discrete Fourier series. Hence, we can also replace ω in π ♭ ωmN lm [z(λ)] λ with ω mN , and then In the second equality, we used the periodicity of the integrand with the period Λ. Finally, using Eqs. (68) and (73), we can reduce Eq. (64) to dQ dt In the second equality, we used the expressions of dE/dt t and dL/dt t [12], The contribution from the radiative field in Eq. (74) is nothing but the scalar analogue of dQ/dt (rad) obtained by Flanagan et al. [17]. We note that the above expression for the radiative part is also valid for the off-resonance case if we interpret the summation over N as meaning the summation over independent frequencies in this case.

C. Practical formula of dQ/dt (sym−S)
t

: the regularized symmetric part
To obtain a more practical expression for dQ/dt (sym−S) t , we rewrite the potential functions of the symmetric and S-parts given in Eq. (63), respectively. The symmetric part can be treated in a similar manner to the radiative part discussed in Sec. III B, though we cannot completely separate the averages over λ and λ ′ at the level of each mode. Substituting Eq. (37) (with Eq. (38)) into the potential function of the symmetric part in Eq. (63), we obtain ] is a periodic function with period Λ, we find that the λ ′ -integral in Eq. (77) produces δ(ω − ω mN ), and therefore we can discretize the Fourier integral as where ǫ 1 = ǫ cos ζ, and ǫ 2 = ǫ sin ζ. In the above equation, the long-time average over λ can be replaced by an average over one period because the part in the angle brackets is periodic with period Λ. For later convenience, in subtracting the the potential function of the S-part, Ψ (S) (∆λ; ǫ, ζ), from Ψ (sym) (∆λ; ǫ, ζ), we rewrite Eq. (78) by rearranging the order of the summation with respect to l and m as Here we define the mode decomposition of Ψ (sym) (∆λ; ǫ, ζ) as Equations (79) and (80) can be interpreted as the Fourier expansion in the (ǫ 1 , ǫ 2 )-space and the Fourier coefficients, respectively. The point is that each (m, N ) Fourier mode is finite, which can be calculated numerically. Next, we derive Ψ (S) (∆λ; ǫ, ζ). Using half the squared geodesic distance between x and x ′ , σ(x, x ′ ), the S-part Green function defined in the local convex neighborhood of the particle's position is written as where L is the scale of curvature of the background spacetime, and V (x, x ′ ) is a regular function in the coincidence limit x → x ′ (see, e.g., Ref. [2] for more concrete definitions of these quantities). The integrand containing V (x, x ′ ) contributes as the finite value to Ψ (S) (∆λ; ǫ, ζ) since the λ ′ integral has its support only in a short interval around λ ′ ≈ λ. In the coincidence limit, this interval shrinks to zero, and thus this contribution also vanishes. Then, we need to consider only the term including U (x, x ′ ).
Abbreviating the argument ∆λ, the argument of the δ-function is expanded as where we have introduced the new variables and assumed δλ = O(ǫ) since the solution of σ = 0 for δλ has the same behavior. We stress that the O(ǫ 3 ) term is absent since σ is unchanged under the transformation ǫ → −ǫ and δλ → −δλ. With the aid of the above expansion for σ, we expand both Σ[z(λ)] and Σ[z(λ ′ )] in terms of ǫ as well, and eventually obtain where we replace the long-time average overλ by the average over one period because the integrand is periodic with period Λ. We note that ψ(ζ) depends on ζ through ξ µ . This is the only contribution of the S-part that remains in the limit ǫ → 0. Since (g µν + u µ u ν ) is a spacelike projection operator, the expression inside the square root is positive semi-definite.
To subtract the S-part of the potential function from the symmetric part mode by mode, we calculate the (m, N )-modes of Ψ (S) (∆λ; ǫ, ζ) as This expression is a double integral with respect toλ and ζ, and is finite. Since ψ(ζ) is independent of m and N , once we obtain ψ(ζ) as a function of ζ by performing theλ-integral, the remaining integral over ζ for each pair of m and N is just a single integral. We stress that subtracting Eq. (87) from Eq. (80) is a variant of the so-called mode sum regularization [1, 27] developed in the context of self-force, but our proposal is suitable for the decomposition in terms of the spheroidal harmonics. We do not have to re-expand the perturbation using the spherical harmonics, which is usually utilized. After subtracting this S-part contribution from Ψ (sym) mN (∆λ) mode by mode, one should be able to take the summations over m and N without any divergence. The convergence of these summations will be even accelerated by fitting the asymptotic form as usual. We will report the result of the explicit convergence test of Ψ

IV. SUMMARY AND FUTURE DIRECTIONS
In this paper, we have derived the long-time averaged rate of change of the Carter constant, dQ/dt t , for a self-interacting scalar charge within the adiabatic regime, which is also applicable to a resonant inspiral. The essential point is that, by contrast to off-resonance orbits, a resonance orbit depends on the offset phase, ∆λ, which is the difference between the times when the rand θ -oscillations reach their minimum values. While this offset phase is identically zero for an off-resonance orbit, it takes a non-zero value for a resonance orbit, and also evolves slowly compared to the time scale of the orbital motion. Therefore, we need to concern ourselves with the evolution of not only the three constants of motion but also ∆λ. The existence of a non-zero ∆λ requires a direct computation of the self-force acting on the charged particle even for calculating the averaged rate of change of the Carter constant. This means that dQ/dt t for a resonance orbit is no longer determined solely by the radiative field Φ (rad) (x) alone, and that we need to consider the contribution from the symmetric field Φ (sym) (x) as well.
As for the radiative part dQ/dt (rad) t , we have derived a practical expression in a similar manner to the off-resonance cases with slight modification, and found that our result is equivalent to the scalar analogue of the long-time averaged rate for the gravitational case shown in Ref. [17]. On the other hand, the symmetric part is more complicated than the radiative one because of the singular behavior of the symmetric field at the location of the particle. We have proposed a method of subtracting the singular behavior, the so-called S-part, from the symmetric part via the point splitting regularization in the Killing directions. The key point is that it is compatible with the Teukolsky formalism, which makes the calculation much simpler than the direct selfforce calculation. While there is an expectation that dQ/dt (sym−S) t might eventually become suppressed, which is suggested by post-Newtonian analysis [28], this conjecture should be tested in a fully general relativistic analysis. Since the numerical code that calculates the long-time average of the constants of motion for generic off-resonance orbits around the Kerr black hole is well developed [29,30], there is in principle no obstacle to numerically testing the conjecture. The code is now under development and the result will be reported in our next publication.
To avoid unnecessary confusion over various works, we comment on the following two points related to our current work. The first point is the difference between the long-time average and the phase space average for a resonance orbit, discussed in Refs. [6,13]. Unlike the long-time average, the phase space average for a resonance orbit means that the averaging is further performed with respect to the slowly evolving variable ∆λ, which is an important control variable of the resonance orbit. The phase space average eliminates all dependence of ∆λ both in dQ/dt (rad) t as well as dQ/dt (sym−S) t , and thus it is likely to miss the important contribution to the accumulated orbital phase correction that would exceed O(1) radians (see the discussion in Appendix A). For this reason, we think that only the long-time averaged rate of change of constants of motion can be a good candidate for obtaining the adiabatic evolution of a resonant inspiral with a phase accuracy better than the order of unity.
The other point is that our new method does not reduce the importance of developing the direct self-force calculation, which has recently undergone rapid progress. We think that these two methods are complementary. Relatively high precision will be required in evaluating the effect that may cause an orbital phase error enhanced by an inverse power of the mass ratio. Numerical evaluation of such quantities might be a tough problem if we cannot take advantage of the Teukolsky formalism that reduces the problem into one to solve one-dimensional radial functions. Even if achieving high accuracy numerically is not so difficult, our method will give reference numbers which are useful to check the accuracy of the results obtained by direct self-force calculations.
Before closing this paper, we would like to mention how to extend our method to the case of the gravitational radiation reaction. Since our method is based on the mode sum regularization in which each mode contribution is gauge invariant, we do not have to concern ourselves with the gauge adjustment between the symmetric part and the S-part. Even for the gravitational case, the contribution from the S-part evaluated in the Lorenz gauge is almost identical to the scalar case discussed in this paper. To calculate the contribution from the symmetric part, we can use the standard Teukolsky formalism with the aid of the established method of reconstructing the metric perturbation from the master variables [31][32][33][34]. The metric perturbation reconstructed in this manner is given in the so-called radiation gauge, in which a singularity runs from the particle's trajectory in the radial null direction. However, almost all regularized trajectories shifted in the Killing directions do not cross this singularity. Therefore we expect that the presence of this radial singularity will not cause any trouble. Although we need further study for practical implementation, we think that our method can be extended to the gravitational case, as briefly described above. the spin angular momentum aM . To apply the following discussion to the scalar model, we just need to replace µ with (q 2 /µ) in the results, where q is the scalar charge of a point scalar particle.
As was discussed in Refs. [5,9], as far as we consider the first-order self-force in the osculating orbit approximation, the long-term orbital evolution is basically determined by tracing the evolution of slowly changing variables, i.e., the constants of motion I i , and the orbital frequencies as their functions, Υ a (I i ). We only need to know the non-oscillating part of the rates of change of I i and Υ a (I i ) for a generic orbit as the corrections due to the self-force. In particular, we only need to know the former if it is allowable to neglect O((M/µ) 0 ) error in the orbital phase. This observation is essentially equivalent to the more systematic two-time scale analysis given in Ref. [6].
Near the resonance, however, the adiabatic evolution of the orbit cannot be described by tracing the constants of motion I i because the self-force also depends on the offset phase ∆λ, the relative phase between the rand θ -oscillations, as discussed in the main text. Neglecting the part that oscillates in the orbital time scale, the time derivatives of I i are determined as functions of I i and ∆λ, which we denote byİ i (I i , ∆λ). At the same time, however, ∆λ also evolves unless the exact resonance condition is satisfied, and the resonance orbit effectively comes back to the off-resonance orbit when the offset phase starts to change rapidly. We examine possible patterns of evolution of {I i , ∆λ} below.
Suppose that the resonance condition is only approximately satisfied, but its deviation is small. In this situation,İ i may still depend on ∆λ. Moreover, we are also allowed to replace I i inİ i (I i , ∆λ) by a fixed value at the exact resonance point, which we denote I i res . This can be understood in the following manner. The typical time scale for crossing the resonance point is estimated as τ res ∼ O(Ω φ M/µ). During the resonance crossing, the net change of the constants of motion from I i res is estimated as [5,6]. This means that the correction to I i res is always suppressed by a factor of O( µ/M ) in the resonance crossing. Based on this observation, we can estimate the error inİ i that is associated with the replacement of its argument from I i to I i res as After passing the resonance, this error term causes the shifts of the constants of motion, which will be evaluated as ∆I i err /I i ∼ O(µ/M ). However, the net phase error due to these shifts, which is simply given by ∆I i err times τ rad , is at most O((M/µ) 0 ). This phase error is small compared with what we are interested in here. Thus, we conclude that the evolution of the constants of motion during the resonance crossing is given byİ i (I i res , ∆λ). Next, we discuss the evolution of the offset phase, ∆λ. In the situation considered here, there is no longer an exact common resonance period Λ. Instead, we introduce the approximate common periodΛ bỹ During the periodΛ, the offset phase is shifted by where ∆Υ := j θ Υ r − j r Υ θ . When ∆Υ = 0, we have an exact resonance. In the adiabatic regime, we are also allowed to replace the time derivative of ∆λ with the averaged value over one periodΛ. As a result, the evolution of ∆λ is described as Now, we investigate the evolution of ∆λ around a resonance point in more detail. For this purpose, focusing on the time derivative of d∆λ/dt t , we consider the two cases separately, depending on whether or not there exists ∆λ c such that d dt (see also Ref. [35] for the classification of resonance orbits). The latter case has already been discussed in Ref. [5] and is essentially the same as the transient resonance analyzed in Refs. [13,17]. As ∆T typically scales like (∝ 1/µ), the phase corrections due to this effect can be significantly large (∝ 1/ √ µ) if the inspiral is in the extreme mass ratio regime, M ≫ µ.
A more interesting situation may happen if ∆λ c exists. This case corresponds to the so-called sustained resonance [13]. When ∆λ c is sufficiently close to ∆λ c , we obtain, with Taylor expansion around ∆λ c , d dt where Since the time scale 1/|̟| scales like (∝ 1/ √ µ), it is likely to be much shorter than the radiation reaction time scale, which scales like (∝ 1/µ). If ̟ 2 > 0, the above system behaves like a harmonic oscillator with its potential minimum varying adiabatically. In this case, ∆λ oscillates around ∆λ c with the amplitude ∝ 1/ √ ̟. This evolution may last as long as the condition ̟ 2 > 0 and the existence of ∆λ c are maintained. Therefore, in principle, the resonance may last for the radiation reaction time scale. This time scale is much longer than we naively expect when the possible presence of ∆λ c is neglected. If ̟ 2 < 0, the value of ∆λ close to ∆λ c is disfavored. The time scale for the variation of ∆λ is again ∝ 1/ √ µ as in the case of the absence of ∆λ c . Hence, the corrections to the phase evolution are also of the same order as before. Although the post-Newtonian analysis [13] and the gravitational wave fluxes from a resonance orbit [17] suggest that the sustained resonance may be absent in inspirals around the Kerr black hole, or may exist only in very restricted orbital parameter regions even if it is achieved, the existence of the sustained resonance in the inspirals into the Kerr black hole is still an open question.
The relation (B2) simply tells that the Green functions depend on t and t ′ through t − t ′ , and on φ and φ ′ through φ − φ ′ , respectively. Based on this understanding, we find that Eq. (58) holds. Next, we show that the first term on the right-hand side of Eq. (55) vanishes by using Eqs. (56), (57), and (58). We focus only on the t-derivative part, but one can show that the φ-derivative part vanishes in the same way. By using the property of Eq.(57), we can replace the long-time average over λ with the average over one period Λ as In the third equality, we transformed (λ, λ ′ ) to (λ − nΛ, λ ′ − nΛ), and used Eq. (57) and Σ[z(λ − nΛ)] = Σ[z(λ)] (note that Σ(x) is a function of r and θ). We also rewrite the integrals over λ and λ ′ as Further, for the Green function G (sym−S) (x, x ′ ), by exchanging the labels in the last line of Eq. (B4) as λ ↔ λ ′ and x ↔ x ′ , we obtain where we used Eqs. (56) and (58) in the last equality. Therefore, we find As for the ∆λ-derivative part, we transform where we replaced the long-time average with respect to λ with the average over one orbital period Λ as shown in Eq. (B3), and exchanged the labels as λ ↔ λ ′ and x ↔ x ′ , and the domains of integration with respect to λ and λ ′ in the same way as in Eq. (B4). Using Eq. (B7), we find (B8) Removing the tand φ-derivative parts from Eq. (55) and using Eq. (B8), we finally obtain Eq. (59).