{'fraction_non_alphanumeric': 0.04139736562218768, 'fraction_numerical': 0.014504271806079872, 'mean_word_length': 4.195652173913044, 'pattern_counts': {'":': 1, '<': 6, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Gyroscopes orbiting black holes: A frequency-domain approach to precession and spin-curvature coupling for spinning bodies on generic Kerr orbits

A small body orbiting a black hole follows a trajectory that, at leading order, is a geodesic of the black hole spacetime. Much effort has gone into computing"self force"corrections to this motion, arising from the small body's own contributions to the system's spacetime. Another correction to the motion arises from coupling of the small body's spin to the black hole's spacetime curvature. Spin-curvature coupling drives a precession of the small body, and introduces a"force"(relative to the geodesic) which shifts the small body's worldline. These effects scale with the small body's spin at leading order. In this paper, we show that the equations which govern spin-curvature coupling can be analyzed with a frequency-domain decomposition, at least to leading order in the small body's spin. We show how to compute the frequency of precession along generic orbits, and how to describe the small body's precession and motion in the frequency domain. We illustrate this approach with a number of examples. This approach is likely to be useful for understanding spin coupling effects in the extreme mass ratio limit, and may provide insight into modeling spin effects in the strong field for non-extreme mass ratios.


I. INTRODUCTION
A. The motion of spinning bodies in curved spacetime Perturbative techniques have proven to be very fruitful for understanding the two-body problem of general relativity. Typically one uses the mass ratio of the system as the small parameter defining the perturbation. An example of such a system is an extreme mass ratio inspiral, or "EMRI," a small compact object (mass µ ∼ 1 − 100 M ⊙ ) whose orbit about a massive black hole (mass M ∼ 10 6 M ⊙ ) evolves (at lowest order) due to gravitational-wave-driven backreaction. Such systems are regarded as important sources of gravitational waves for low-frequency, space-based gravitational-wave detectors such as eLISA [1]. Modeling EMRI systems has been a major motivator for much of the progress this field has seen. More fundamentally, such a system represents a clean limit of the general two-body problem. Many quantities which describe the evolution of binaries can be solved in this limit quite precisely. Although formally applicable only to systems with large mass ratios, perturbation-theory-based analyses have proven to be useful for understanding binary systems even for mass ratios of order unity [2][3][4].
At zeroth order in the mass ratio, the small body moves on a geodesic of the background spacetime: The operator D/dτ is a covariant derivative evaluated along the small body's trajectory, with τ proper time along that trajectory. Corrections to geodesic motion enter at higher order in the mass ratio. Schematically, we write Dp α dτ = f α . (1. 2) The force f α that pushes the smaller body away from a geodesic can arise from a number of physical effects. A particularly important and interesting example is the gravitational self force, f α sf . It arises from the small body's own contribution to the binary's spacetime, and enters at order (µ/M ) 2 . This force encodes the loss of energy and angular momentum from the binary due to gravitational waves, as well as smaller scale but observationally important shifts to observables like orbital frequencies. Many groups have put considerable effort into computing the gravitational self-force; see Refs. [5] and [6] for a very detailed overview and discussion of the selfforce problem in general relativity, and Refs. [7][8][9][10][11][12][13][14] for some examples of recent progress in the field.
Other forces arise when we consider that the small body is not a point mass, but has internal structure. The leading example of such structure is spin. Spin couples to the curvature of the background spacetime, modifying the smaller body's motion [15][16][17][18]: These are the Mathisson-Papapetrou equations, commonly called the Papapetrou equations. Here, R α µλσ is the Riemann curvature of the background spacetime (which we'll take to be a Kerr black hole), u α = dx α /dτ is the orbiting body's four-velocity (tangent to its orbit worldline), and p α is its four-momentum. The tensor S λσ describes the spin of the small body; we define it precisely in Sec. III. If the small body is itself a Kerr black hole, then S λσ ∝ sµ 2 , where s ≤ 1 is the small body's dimensionless spin parameter, and where the small body's mass µ will be defined precisely in Sec. III A. In this case, typical components of the spin will have magnitude S ≡ sµ 2 ≤ µ 2 , so the leading effect of the spin-curvature force will be of order µ 2 . A major motivation for our analysis is that nongeodesic effects associated with the self force and with spin-curvature coupling both lead to forces that scale at leading order as µ 2 (at least if the small body is considered to be a Kerr black hole, in the case of spin-curvature coupling). This very crudely suggests that the effects of these two forces might be comparable. However, it should be emphasized that this scaling argument can only motivate a more careful analysis which ultimately examines gauge invariant, observable quantities. It is only through such an analysis that one can ascertain whether these forces drive secularly accumulating phase shifts (for example), or produce periodic signatures that may be searched for in astrophysical data. The work in this paper is the first step in a program to develop such an analysis, laying the foundations for an examination of such effects for generic orbits of spinning bodies with arbitrary spin-orbit orientations.
Equations (1.3) and (1.4) must be supplemented by an additional constraint, reflecting an internal degree of freedom associated with the small body's structure. One can reference an extended body's motion to an infinite number of worldlines which pass through it. The additional constraint which must be imposed can be understood as specifying which of these worldlines one uses. This constraint is called the spin-supplementary condition, or SSC, in the case of a body described by its mass and spin. Figure 1 of Ref. [19] provides a very clear and simple illustration of why the SSC is needed. In non-relativistic dynamics, a natural choice for the worldline that one uses to describe a body's motion is the one which passes through its "center of mass," or COM. However, in relativistic dynamics, the COM is observer dependent. An observer who is at rest with respect to the axis of rotation of a uniform spinning sphere will place the COM at the sphere's geometric center. An observer moving relative to the rotation axis will not: thanks to the relativistic transformation of mass currents, and thanks also to the fact that different mass elements of the sphere have different relative velocities to this observer, the moving observer places the spinning sphere's COM away from the rotation axis.
In this language, the SSC can be regarded as choosing a particular observer and tying the body's worldline to their choice of COM. Costa and Natário [19] provide excellent discussion and many further details elaborating on this point. An SSC which is commonly used in the literature is the Tulczyjew SSC [20], (1.5) Several variations on Eq. (1.5) have been formulated, and are compared by Kyrian and Semerák in Ref. [21]. Many SSCs 1 differ from one another at O(S 2 ); if the smaller body is itself a Kerr black hole, this difference is O(µ 4 ). These differences will be extremely small for large mass ratio systems, suggesting that it would be useful to linearize our system in spin. Related to the internal degree of freedom is the fact that the four-momentum and four-velocity are in general not parallel to one another when the small body's spin is taken into account. Instead, we have (1.6) the momentum and velocity are only parallel if the second term vanishes (which in fact it does for some SSCs; see Sec. 3.4 of Ref. [21] for an example). For several cases (including the Tulczyjew SSC, which we will use), In such a case, p α and u α are parallel at order S. Given an initial position, velocity, and spin configuration, Eqs.  (1.5) can be ignored. The trajectories they describe reduce to Kerr geodesic orbits. For S small, these equations describe a nearly geodesic orbit, modified by a spin-curvature force of order S = sµ 2 . As already mentioned, the scaling of the spin-curvature force with mass is naively similar to the scaling for the gravitational self force, motivating a more careful analysis to ascertain whether its effects are of observational relevance. Examples of such analyses have been done for special cases. Burko [22] and Burko and Khanna [23] find that the spin-curvature force may in fact dominate over certain aspects of the self force for circular and equatorial orbits, at least when the small body's spin is large and aligned with the orbit; similar results were found by Steinhoff and Puetzfeld [24].
Although there has been a great deal of work studying motion under the influence of the Papapetrou equations (see Kyrian and Semarák [21] and references therein for an outstanding review), the impact that the small body's spin has for potentially relevant EMRI sources has not yet been carefully quantified. It appears that, to date, the most careful and rigorous analyses of spin-enhanced, large-mass ratio dynamics have focused upon simple orbits and simple spin-orbit configurations. Most work on this problem has focused upon circular orbits with spin parallel to the orbital angular momentum [24,[56][57][58][59]; at least one paper [60] considered the small body's spin to be slightly misaligned from parallel, leading to precession kinematics which influenced the system's emitted gravitational waves. A very recent analysis using a Hamiltonian formulation of spin-curvature coupling has the promise of analyzing motions for general orbits and spin orientations [61], though at present the authors confine their analysis to spherical background spacetimes and planar orbits. With one exception, analyses which focus on issues of gravitational-wave measurement and data analysis likewise use the circular, equatorial, parallel spin limit [22,62]. The one exception is Barack and Cutler [63], who developed a "kludge" model for the inspiral and waveforms by pushing the post-Newtonian expansion beyond its formal range of validity in order to assess parameter measurement accuracy with spacebased gravtiational-wave detectors. Barack and Cutler found that the small body's spin has a negligible impact on parameter accuracy, except when the small body has maximal spin. It will be important to revisit this analysis using techniques that accurately describe the small body's motion and precession deep into the strong field. C. Our approach to spin-enhanced orbital motion A major goal of the work we present here is to develop tools that can be used to efficiently compute orbits that include the impact of spin-curvature forces for EMRI systems with arbitrary spin-orbit configurations. In this paper, we focus on computation of the precession and forces for EMRIs; later work will examine how to evolve orbits using the resulting force. We assume µ ≪ M , where M is the mass of the large black hole, and µ is the mass of the orbiting body. We assume that the smaller body has a spin S, and define the dimensionless spin parameter s = S/µ 2 . If the small body's spin obeys the Kerr bound, then 0 ≤ s ≤ 1.
Our approach to modeling spin-enhanced orbital motion is based on the fact that, for large mass ratio systems, the motion will be "close to" a Kerr geodesic, with the spin-curvature force introducing corrections to the orbiting body's worldline of order µ/M ≪ 1. This suggests that a useful way to compute the force's effect is with the method of osculating geodesics [64,65], in which the body is regarded as evolving through a sequence of geodesic orbits. This sequence of geodesics shares the tangent to the body's true worldline at each moment. Osculating geodesics generalize the technique of osculating orbits that are commonly used in astrophysical dynamics, and are a powerful way to describe motion which is perturbed away from an exact, integrable trajectory.
It should be emphasized at this point that the equations of motion for a spinning body on a generic orbit about a Kerr black hole are not integrable 2 . One can in fact find configurations that lead to chaotic orbital evolution [66][67][68]. Our proposed method of treating the motion as a sequence of osculating geodesics presumes that the motion is "nearly" integrable, in the sense that the non-integrable spin-curvature coupling acts as a weak perturbation to the integrable geodesic motion. This is presumably reasonable provided that the small body's spin is "small enough," though it remains to be quantified what "small enough" really means.
In our formulation of this problem, we use the fact that functions arising from motion on a Kerr orbit can be expanded in a Fourier series [69]. A Kerr black hole orbit is triperiodic, with separate frequencies describing its radial, polar, and axial motions. Many functions which arise from motion along such an orbit, including those governing spin-curvature coupling, can be written where z is the worldline followed by the orbiting body, t is time as measured by a distant observer, and Ω θ,r are frequencies characterizing the orbit's motion in the θ and r directions. (A third frequency associated with axial motions, Ω φ , exists, but is not important in this analysis thanks to the axisymmetry of the Kerr spacetime.) Similar forms can be written down using different time parameterizations. Including the small body's spin introduces a new frequency Ω s which is associated with the small body's precession as it moves along its orbit. These frequencies all vary smoothly with the parameters that describe bound Kerr orbits, and so evolve in a simple way as the worldline moves through its sequence of osculating geodesics. We find that frequency-domain expansions of the quantities needed to model spin-enhanced motion tend to converge rapidly, at least for orbital eccentricities e 0.7 or so 3 . This approach should thus facilitate efficient calculation of the quantities needed to compute spin-enhanced orbits. Describing these quantities in the frequency domain can be useful for understanding how harmonic features of the orbital and precessional motion are imprinted on observationally important characteristics of the system. We expect that a frequency-domain description will prove useful for extending codes which compute radiation from generic Kerr orbits (e.g., [70]) to include the effects of spin. Work of this nature has been done so far for simpler orbit and spin configurations [56-58, 60, 71].

D. Organization of this paper
In this paper, we show how to compute the new frequency associated with the small body's spin and examine the Fourier decomposition of all the relevant quantities needed to compute spin-enhanced EMRI orbits. Since our analysis is based on expanding about geodesics of Kerr black holes, we begin in Sec. II by reviewing the properties of these geodesics. We introduce the notation that we use, and describe briefly the parameterization and properties of Kerr black hole orbits. Chief among these properties is the periodic structure of these orbits. We review the triperiodic structure of Kerr black hole orbits, and describe how any function which arises from motion along such an orbit can be expanded in a Fourier series of the orbit's fundamental frequencies. In Sec. III, we move beyond geodesics, summarizing how the small body's spin couples to spacetime curvature. We first review the general equations that govern spin-curvature coupling in Sec. III A, and then examine in Sec. III B how these equations simplify when we linearize in the small body's spin. 3 Based on experience studying gravitational waves from eccentric binary systems (e.g., Ref. [70]), we expect that the frequencydomain expansion will converge very slowly as e → 1. We have not, however, done a systematic study to determine the value of e at which the frequency-domain approach becomes less useful.
In Sec. IV, we then study the precession of the small body's spin in the frequency domain. We do this analysis in what is often called "Mino time," a time variable that fully separates the radial and polar motions of Kerr black hole orbits, and that has proven to be well suited for studies of strong-field motions. We begin by studying precession along a circular equatorial orbit (Sec. IV A). This limit turns out to be very clean when analyzed using Mino time. Changing to the Boyer-Lindquist time coordinate (appropriate for analyzing quantities that distant observers would measure) and taking the weak-field limit, we recover well-known results describing gyroscope precession in general relativity. We then generalize to more complicated orbit geometries in Secs. IV B and IV C. In all these cases, we find that the frequency characterizing precession is the eigenvalue of a matrix representation of the relevant equations. This representation can become fairly complicated for generic orbits, thanks to the complicated time-frequency structure of the underlying orbits. Nonetheless, we find it is not terribly difficult in this framework to accurately construct the solution for the small body's precession.
With the small body's precession in hand, we then examine the spin-curvature force (1.3) in Sec. V. We find that the frequency-domain expansion converges quickly, and is an excellent way to describe these quantities in the large mass ratio regime. Interestingly, there exist certain configurations in which the precession and orbit frequencies are "commensurate" with one another: some linear combination of the precession and orbit frequencies sums to zero. Fourier components of the force corresponding to such combinations would be constant in time, and one might imagine that the system's evolution would be substantially modified as compared to "nearby" orbit configurations (as has been seen in studies of the self force's behavior at black hole orbit resonances [72,73]). In Sec. VI, we show that, at least at linear order in spin, resonant spin-curvature coupling does not have any special impact on the system's evolution. This is due to the existence of conserved constants of the motion. To "protect" these constants, certain contributions to the spin-curvature force vanish when the precession and orbit frequencies are commensurate. We speculate that such resonances (or an analog of such resonances) might have more interesting effects if the analysis were taken to higher order in the small body's spin, and might be connected to the onset of chaos that has been seen in past analyses of the orbits of spinning bodies [66][67][68].
We conclude in Sec. VII with a sketch of our plans for future work on this problem. Chief among these plans is to marry this scheme for computing the motion of a small body to an osculating geodesic integrator [64,65] in order to compute the full worldline describing a small spinning body orbiting a Kerr black hole. This analysis is underway, and early results indicate that it is a robust and accurate way to describe these orbits. We then hope to use it to broaden the class of sources that can be incorporated into future EMRI models, and to model two-body dynamics in general relativity more broadly.

II. BACKGROUND: CONVENTIONS, NOTATION, AND GEODESIC MOTION
We begin by laying out the conventions and notation that we will use, and describe in some detail the Kerr geodesics that constitute the background motion we use for orbits of spinning bodies. This material has appeared in great depth in quite a few other sources, so we simply summarize what is relevant for us here and point the reader to relevant literature.

A. Generalities
Our goal is to understand the motion of a small body of mass µ orbiting a Kerr black hole with mass M and spin parameter a. We will use Boyer-Lindquist coordinates for our analysis, in which the metric is (We use geometric units with G = 1 = c throughout this paper.) From Eqs. (2.1) and (2.2), it is straightforward to compute the Kerr metric's connection coefficients and Riemann curvature components in these coordinates, both of which we need in explicit form. Semerák provides a catalog of these quantities (Ref. [74], Appendix A; see also footnote 2 of Ref. [21]). The Kerr spacetime admits two Killing vectors, ξ α t and ξ α φ , corresponding to timelike and axial diffeomorphisms. It also admits a Killing tensor K µν , given by is a Newman-Penrose null tetrad leg 4 , andm µ is the complex conjugate of m µ . This tensor satisfies a generalized form of Killing's equation, (2.5) 4 Here and below, we use " . =" to mean "the indexed object on the left-hand side is represented by the components on the right-hand side in Boyer-Lindquist coordinates." To our knowledge, K αβ does not have a simple geometrical interpretation as the Killing vectors do. Note that K αβ is symmetric under exchange of indices. For later discussion, it is useful to note that we can write where the anti-symmetric Killing-Yano tensor f αβ is given by [60] f µν = a cos θ e 1 µ e 0 ν − e 0 µ e 1 ν + r e 2 µ e 3 ν − e 3 µ e 2 ν , (2.7) and where The Killing-Yano tensor satisfies Associated with the Killing vectors and the Killing tensor are three constants which are conserved along a geodesic worldline. These are the orbit's energy, axial angular momentum, and "Carter constant": 14) The name "Carter constant" is commonly used for both K G and Q G . In all of these quantities, the superscript "G" indicates that these are quantities associated with motion along a geodesic. This is to be contrasted with analogs we will introduce in the next section which are associated with a more general class of trajectories. In much of our analysis, it is useful to normalize these quantities to the rest mass of the orbiting body µ. We indicate their values so normalized with a hat: (2.17) The existence of these constants of motion makes it possible to separate the equations governing Kerr geodesics. In first order form, they become [Ref. [75], Eqs. (33.32a)-(33.32d)] (2.21) These equations use proper time τ as the independent parameter along the geodesic. Another time parameter which is very useful for studying strong-field Kerr black hole orbits is λ, defined by dλ = dτ /Σ. The geodesic equations parameterized in this way are By using λ as our time parameter, the r and θ coordinate motions completely separate. The parameter λ is often called "Mino time," following Mino's use of it to untangle these coordinate motions [76].

B. Parameterizing geodesics
We have found it very useful to introduce the following reparameterization of r and θ: , cos θ = cos θ m cos(χ + χ 0 ) .
(2.24) We further find it helpful to remap θ m to θ inc , defined by This angle encodes whether the orbit is prograde (L z > 0) or retrograde (L z < 0). It smoothly varies from θ inc = 0 • for prograde equatorial (θ m = 90 • , L z > 0) to θ inc = 180 • for retrograde equatorial (θ m = 90 • , L z < 0). There is a fairly simple mapping between the constants of motion (Ê G ,L G z ,Q G ) and the orbital geometry parameters (p, e, θ m ). See Appendix B of Ref. [77] for explicit formulas relating these two parameterizations; Ref. [78] also provides many valuable results for studies of Kerr geodesic motion. Notice that r oscillates between r min and r max , given by (2.26) Likewise, θ oscillates from θ min = θ m to θ max = π − θ m . The transformations (2.24) replace the variables r and θ with secularly accumulating angles ψ and χ. As ψ and χ evolve from 0 to 2π, r and θ move through their full ranges of motion. By combining Eq. (2.24) with Eq. (2.22), it is straightforward to develop equations for dψ/dλ, and dχ/dλ. In this parameterization, a geodesic worldline is defined by a set of four functions [ψ(λ), χ(λ), φ(λ), t(λ)], three orbital parameters (p, e, θ m ), and four initial conditions. We define χ = ψ = 0 at λ = 0; the value of r and θ at λ = 0 is then set by the angles ψ 0 and χ 0 defined in Eq. (2.24). We likewise define t = t 0 and φ = φ 0 at λ = 0.

C. Geodesics and functions of geodesics in the frequency domain
Bound Kerr orbits are triperiodic, as was first explicitly shown by Schmidt [77]. It is simplest to compute the three orbital periods in Mino time, which (as discussed above) separates the radial and polar motions. Using this time variable, it is not too difficult to compute three periods Λ r , Λ θ , and Λ φ associated with the radial, polar, and axial motions, respectively, as well as associated frequencies Υ r,θ,φ = 2π/Λ r,θ,φ . One can also compute a quantity Γ which converts these periods and frequencies to their Boyer-Lindquist coordinate time analogues: are the periods in Boyer-Lindquist time, and the associated frequencies are Since Boyer-Lindquist time corresponds to time as measured by distant observers, this formulation is useful for describing the time-frequency behavior of observable 5 Due to periodicity, λ = λ r0 ± 2πjΛr will also satisfy this condition for any integer j (where Λr is the Mino-time radial period). A similar statement holds for λ θ0 .
quantities. Schmidt [77] first showed how to compute the frequencies Ω r,θ,φ using elegant but somewhat abstract techniques. Drasco and Hughes [69] showed how to understand these periodicities in simpler terms, with easyto-compute quadratures describing Λ r,θ,φ and Γ (from which computing Ω r,θ,φ is straightforward). Fujita and Hikida [78] then showed that these quadratures can be evaluated analytically, yielding easy-to-use formulas for Υ r,θ,φ and Ω r,θ,φ as functions of p, e, and θ m . These frequency-domain expansions are excellent tools for characterizing the behavior of functions associated with geodesic black hole orbits. The Mino-time periodic expansion is particularly good for describing strongfield orbital dynamics, and will be used extensively 6 in this paper. Let f (λ) ≡ f [r(λ), θ(λ)] be a function of r and θ that is computed along a geodesic worldline (2.30) A more complicated expansion must be used for functions which depend in addition on t and φ; see Ref. [69] for discussion. We will not need this more complicated form for this paper. Let us denote byf kn the amplitude (2.30) computed along the fiducial geodesic, χ 0 = ψ 0 = 0. Oncef kn is known, it is simple to compute the amplitude along an arbitrary geodesic [79]: This means that one only needs to computef kn ; the Fourier amplitudes along all other geodesics can be obtained from this quite simply.

III. SPIN-CURVATURE COUPLING
Strictly speaking, geodesic motion only applies to a zero-size point body moving through spacetime (and neglecting self force effects). If the body has any structure beyond the point particle description, that structure will couple to the spacetime through which it moves. This coupling will push the small body away from the geodesic, appearing as a force driving the body's motion.
Papapetrou [17] pioneered a moment-based technique for determining these couplings (see also Mathisson's discussion in Ref. [16]). Following the discussion in [80], the key idea is to choose some representative worldline that passes through a body moving through spacetime. One computes moments of the body's stress-energy tensor T αβ about that worldline. If the body is sufficiently compact that moments beyond the nth can be taken to vanish, then one can find equations of motion in terms of those non-vanishing moments by enforcing ∇ α T αβ = 0. If all moments beyond n = 0 vanish (meaning that the small body is a monopole point mass), then the geodesic equations for motion in the background spacetime emerge from this procedure. The first moment to couple beyond this describes the small body's spin angular momentum. Spin couples to the curvature of the background spacetime, producing a force which pushes the small body away from the geodesic. One can consider the coupling of higher order moments as well [81]; for example, the small body's quadrupole moment couples to the gradient of the curvature [82,83].
Here we discuss the spin-curvature coupling and the force which results from it. Further detailed discussion and derivations can be found in Refs. [15][16][17][18]80]. We first summarize the equations which govern the motion of a small spinning body through curved spacetime in full generality, making no assumptions about the magnitude of the small body's spin or the worldline that it follows. We then introduce a perturbative expansion, developing these equations at leading order in the small body's spin, and treating its four-velocity as "close to" geodesic motion (in a sense that is quantified precisely below). A similar approach was developed in Refs. [84,85], though without the goal of then connecting this perturbative force to an osculating geodesic integrator.

A. General form
We begin by defining the spin tensor, which is central to our analysis. Let δz α = x α − z α , where x α is a general point in spacetime, and z α is position along the small body's worldline. The spin tensor is defined as In these equations, T αβ is the stress-energy tensor of the orbiting body, and square brackets around indices denote antisymmetrization: In the covariant form (3.1), the spin is defined on an arbitrary spacelike hypersurface Σ, and x α are coordinates in that surface. In Eq. (3.2), we have taken Σ to be a "slice" of constant t, where t (not necessarily the Boyer-Lindquist time) is used to parameterize the orbiting body's worldline; g is the determinant of the metric. By enforcing ∇ α T αβ = 0, we find the following equations governing the motion of the small body [16][17][18]: The operator D/dτ denotes a covariant derivative along the worldline that the small body follows. The fourvelocity is defined, as usual, by u α = dx α /dτ . However, it is not the case that p α = µu α (where µ is the small body's rest mass). As discussed in Sec. I, Eq. (3.3) and (3.4) do not fully determine the small body's motion. One must also impose a spin-supplementary condition (SSC), which accounts for degrees of freedom which are implicit in its non-point-like structure. An outstanding review and discussion of this condition and its physical meaning is given by Costa and Natário [19]; further excellent review and comparison is given in Kyrian and Semerák [21], and additional comparison of spin-supplementary conditions is provided by Ref. [86].
A commonly used choice is the Tulczyjew SSC [20], We use Eq. (3.5) for our analysis. We emphasize that this is an arbitrary choice, and is adopted primarily because it is often used in literature that examines gravitationalwave generation from spin-enhanced orbits (e.g., Refs. [24,[56][57][58][59][60]). Equivalent motion can be shown to follow from a range of SSCs, though one must be careful to compare properly, as described in detail in Ref. [19]. Using Eq. (3.5), one can show that The parameters µ and M both have the dimensions of mass, and are related 7 to u α and p α by Note that u α u α = −1. Using Eq. (3.5), one can show that µ is constant along the orbiting body's worldline, but that M is not. From these definitions, it is simple to see that µ = M + O(S 2 ), and that u α = p α /µ + O(S 2 ). The mass parameters µ and M are identical for geodesic orbits (as they must be), and cannot be distinguished at linear order in the small body's spin. 7 We have reversed the definitions of µ and M relative to how they are defined in Ref. [21]. We do this so that −p µ pµ has the same name in both the spinning and spinless cases.
From the spin tensor, one finds the spin vector [21] Here, where the metric determinant √ −g = Σ sin θ for Kerr, and where [αβγδ] is the totally antisymmetric symbol: [αβγδ] = +1 if αβγδ is an even permutation of 0123 = −1 if αβγδ is an odd permutation of 0123 = 0 otherwise .
Notice that p µ S µ = 0. Less obviously, is a constant. Equations (3.3), (3.4), and (3.5) can be integrated to build the worldline of a spinning body. The worldline which one finds from such an integration admits a constant of the motion for each of the spacetime's Killing vectors: If ξ α is one of these Killing vectors, then is constant along the spinning body's worldline. For the Kerr spacetime, these constants are 14) These quantities define conserved energy and axial angular momentum for a spinning body orbiting a Kerr black hole; the superscript "S" distinguishes them from the conserved quantities associated with geodesics. In general, there is no conserved quantity analogous to the Carter constant for orbits of spinning bodies.

B. Leading order in small body's spin
We now imagine that the small body's spin can be taken to be "small" in some meaningful way. Let us define the dimensionless spin parameter s by where S, defined in Eq. (3.12), is the magnitude of the spin vector. If the small body is itself a Kerr black hole, we expect that s ≤ 1. Even for non-Kerr small bodies, this bound is likely to be a useful guide as long as the small body is compact. We thus expect that spin effects will have an impact on the motion at order µ 2 . It's worth noting other literature uses different conventions to normalize S. For example, Refs. [66] and [71] define In these works, the authors are not restricting their analysis to the large mass ratio limit, so µ is the reduced mass of the system, and M is its total mass. If the binary's secondary is in fact a Kerr black hole, then this condition implies that σ ≤ (µ/M ). It should be noted that these conditions are somewhat arbitrary, since it is the spin S that matters in the equations of motion, not the parameters s and σ. The definition of these dimensionless parameters is simply a convenience to guide our intuition.
The scaling of S with µ for the case of the small body being a Kerr black hole suggests that, especially for large mass ratio systems, it would be fruitful to neglect terms that are of O(S 2 ) and higher. With our convention for the scaling of the small body's spin, the terms we are neglecting would introduce effects at fourth order in the small body's mass. Truncating at O(S), we find It follows that so that the spin tensor is parallel transported along the small body's worldline at O(S).
We now rewrite the equation of motion Eq. (3.3) to leading order in S. We emphasize again that the derivative operator D/dτ is a derivative along the small body's worldline, and that our goal is to develop a force term that can be used with an osculating geodesic integrator. As such, our small body's motion will be regarded as tangent to some geodesic at every moment along its worldline. The geodesic to which the worldline is tangent will differ at each step, but there will be some "reference geodesic" defining D/dτ at each moment.
We begin by writing the small body's 4-velocity where u α G satisfies the geodesic equation, and where u α Equation (3.22) defines the force which, in the osculating geodesic picture, carries the small body's worldline from one reference geodesic to another. It is worth emphasizing that Eq. (3.22) only makes sense in this picture, since the derivative operator D/dτ is defined quite strictly with respect to this geodesic. The results that we present in Sec. V B for the spin-curvature force should thus be regarded as the force tending to push a spinning body away from a specified reference geodesic. The SSC (3.5) becomes when we linearize in S. The spin vector (3.9) becomes This equation can be inverted: It follows from Eq. (3.24) that so that the spin vector is parallel transported along the body's worldline at O(S). It also follows that at this order. This means that we only need to know three of the components of S α ; the fourth component is determined by the constraint. Equation (3.27) tells us that, in the orbiting body's rest frame, the spin vector is purely spatial. At linear order in S, the motion acquires a new conserved constant [60,87]: where f µν is the Killing-Yano tensor introduced in Eq.
(2.7). Although there is in general no analog of the Carter constant for spinning bodies orbiting Kerr black holes, there is such an analog at leading order in the smaller body's spin.
The leading-order-in-spin approach we have described is designed to be used for integrating orbits using the technique of osculating geodesics [64,65]. The small body's true worldline is, at every moment, tangent to some geodesic, the "osculating" geodesic. Using Eq. (3.22), we compute the force which pushes the small body away from that osculating geodesic. We use this to step the system by some amount ∆τ along its worldline, updating the spin via Eq. (3.26). At the end of this step, we update the osculating geodesic, and then repeat.
Integrating orbits with this scheme, one will find that the geodesic values for the energy, axial angular momentum, and Carter constant will oscillate as the system moves between different osculating geodesics. However, the spin-enhanced energy, axial angular momentum, and Carter constant, given by Eqs. (3.14), (3.15), and (3.28), will all be constant.

IV. FREQUENCY DOMAIN TREATMENT OF SPIN PRECESSION
We now expand on the details of our approach to analyzing spin-curvature forces. In addition to treating the small body's spin as a perturbative parameter, we take advantage of the fact that Kerr geodesics are triperiodic, so any function arising from these geodesics can be computed in a Fourier expansion. In the next section, we will use this to show how the components of the spin curvature force behave in the frequency domain.
We start by analyzing the precession of the small body's spin. Precession introduces new frequencies into our analysis, which can be represented as eigenvalues of a matrix representation of the master precession equation. Once the precession frequency is known, it is straightforward to include it in a frequency-domain expansion of the spin-curvature force.
Begin by expanding the covariant derivative in Eq. (3.26), and changing the independent parameter from proper time τ to Mino time λ. The "master equation" governing precession becomes This can be regarded as an eigenvector equation. To expedite solving it, write where S . = (S r , S θ , S φ ), and where P is a 3 × 3 matrix which depends only on r and θ. Here and throughout the paper, boldface symbols denote vectors and matrices whose indices are associated with the spatial coordinates r, θ, and φ. It should be clear from context which quantities are two-index matrices, and which are one-index vectors.
As discussed in Sec. III, the timelike component S t is determined by the constraint u G α S α = 0. Taking this into account, the elements of P are given by Each element of this matrix varies with time as the small body moves through its orbit. Because the matrix is itself a function of Kerr geodesic motion, it can be expanded in the frequency domain, so we write Harmonics of Υ φ do not enter due to the spacetime's axisymmetry. Each Fourier component P kn is itself a 3 × 3 matrix. The indices k and n formally run from −∞ to ∞, although in practice the sums converge at finite value. We now discuss how to solve this. There are three cases which we discuss separately: circular equatorial orbits, for which we need only the k = n = 0 terms; eccentric equatorial, circular inclined, or "resonant" orbits, for which we need one index; and generic orbits, for which we need both.

A. Circular and equatorial orbits
When the osculating geodesic is circular and equatorial, the small body's precession has a particularly simple analytic solution. We develop this solution first in a straightforward way, and then re-analyze this solution as an eigenvector problem. This allows us to begin developing the tools we will need to study more complicated orbit geometries.

Analytic solution for the spin
Consider an orbit of constant radius r, with constant polar angle θ = π/2. The precession matrix P is quite simple in this limit: For an equatorial orbit, S θ is the component of the small body's spin normal to the orbital plane. Equation (4.10) tells us that this spin component is constant for circular, equatorial orbits. As we'll see in Sec. V B, the product rS θ is constant for all equatorial orbits, consistent with this result. For the other two components, we combine Eqs. (4.9) and (4.11) yielding The values ofÊ G andL G z appearing in P r φ and P φ r are given by [88] The upper sign is for prograde orbits, the lower is for retrograde. Combining Eqs. (4.6), (4.7), (4.13), and (4.14), we find the somewhat remarkable simplification In terms of Mino-time λ, this means that for circular equatorial orbits, S r,φ undergo simple harmonic oscillations at the frequency The "CE" superscript is a reminder that this quantity only applies to circular and equatorial orbits. By matching to initial conditions, it is a straightforward exercise to construct S r,φ (λ).
Using dλ = dτ /Σ, with Σ = r 2 for an equatorial orbit, we convert to frequency conjugate to proper time along the orbit: Finally, using Eqs. (2.21), (4.13), and (4.14), we can convert this to a frequency conjugate to observer time: (4.18) As before, upper sign labels prograde orbits, and lower labels retrograde.
The rate at which the spin vector is seen to precess is given by the difference between this frequency and the orbital frequency. Using we have The term at O(r −5/2 ) reproduces well known results for the geodetic precession of an orbiting gyroscope; the term at O(a/r 3 ) reproduces the "gravitomagnetic" or Lense-Thirring correction to this precession arising from the larger body's spin. See, for example, Eqs. (3.3) and (3.4) of Ref. [25] (noting that, for an equatorial Kerr orbit, Barker and O'Connell's n and n (2) are parallel to one another, and point along the black hole's spin axis).

The spin as an eigenvector
Let us now reorganize the above analysis, introducing notation that will generalize to more complicated orbit geometries. We begin by writing S = S a e −iΥ a s λ . In other words, S a is an eigenvector of P with eigenvalue −iΥ a s . Using Eqs. (4.6), (4.7), and (4.8), we find the following three eigenvectors and eigenvalues: where Υ CE S = √ M r as before. A general spin vector can then be written with c −1,0,1 determined by initial conditions. As we saw above, the θ component of S is constant, while the components S r and S φ undergo simple harmonic oscillations [see Eqs. In the circular equatorial case, the matrix P does not vary with time, so the frequency-domain expansion (4.5) is trivial. We now consider cases in which the matrix varies, but is characterized by one dynamically important frequency.
If the orbit is circular but inclined, then the angle θ varies over the orbit. Only the polar frequency matters in this case, and (4.5) simplifies to If the orbit is equatorial but eccentric, then the radius r varies over the orbit, and Eq. (4.5) becomes In the general case, we expect both radial and polar frequencies to be important, and we must use Eq. (4.5) with no simplifications. However, for the special case of resonant orbits [72,73], the two frequencies are simply related to one another: with β r and β θ both integers. A harmonic kΥ θ + nΥ r = N Υ, with N = kβ θ + nβ r . In this case, the expansion takes the form This sum runs from −∞ to ∞, but the spacing between frequencies will depend in detail on the nature of the resonance.
In these three cases, the expansion can be written 8 for an appropriate choice of Υ x . To solve Eq. (4.3) with this form of P, we assume that there exist solutions S a of the form The index a labels different eigenvectors and eigenvalues.
The common factor of e −iΥ a s λ on both sides cancels. Multiply both sides by e iqΥxλ , and integrate from λ = 0 to 2π/Υ x : Performing the sums over j ′ , this simplifies further: Let us truncate the sum over j at some finite value ±j max . Expand this equation and rearrange: P −jmax · S a q+jmax + P −jmax+1 · S a q+jmax−1 + · · · + (P 0 + iqΥ x I) · S a q + · · · + P jmax−1 · S a q−jmax+1 + P jmax · S a q−jmax = −iΥ a s S a q . (4.35) The eigenvalue Υ a s must satisfy this for each of the 2j max + 1 Fourier components. To solve Eq. (4.35), we rewrite it as a bigger matrix equation. We first assemble a "vector of vectors" by combining the different Fourier components S a j as follows: There are 2j max + 1 elements in S a , each of which is itself a 3 element spin vector. We likewise define the (2j max + 1) × (2j max + 1) matrix of matrices P whose elements are each 3 × 3 precession matrices: (4.37) The indices g, h ∈ [−j max , j max ], and I is the 3 × 3 identity matrix. As a concrete example, for j max = 2 this matrix is To find the solution for the precessional motion of the small body, we then solve for the eigenvectors and eigenvalues of the system P · S a = −iΥ a s S a . (4.39) When we do this we find find 3 × (2j max + 1) eigenvalues and eigenvectors. On physical grounds 9 , we expect only 3 eigenvalues and eigenvectors, so we appear to have far more solutions than are needed.
On careful analysis, we find that these many solutions are simply related to one another: they fall into 3 groups, each of which has 2j max + 1 members. The eigenvalues in each group are simply shifted from one another by some multiple of Υ x . The surfeit of solutions originates in a relabeling ambiguity in the eigenvector expansion. Take Eq. (4.31) and shift the index j ′ by some integer ∆j ′ . The resulting eigenvector S a is unchanged by this shift if we take and likewise shift the eigenvector components: Strictly speaking, this shift leaves the system unchanged only in the limit j max → ∞. However, provided j max is large and we confine ourselves to solutions with ∆j ≪ j max , the different precessional solutions that we construct with these related eigensolutions will differ negligibly from one another. We emphasize here that this relabeling ambiguity means that in principle any eigenvalue/eigenvector set can be used to construct a valid solution, provided one correctly accounts for the index shift. In practice, we find that the set which has Υ 0 s = 0 and Υ −1 s = −Υ 1 s is computed most accurately, and is the one we have used in our analyses.
An example usefully illustrates how this relabeling ambiguity works. Consider an equatorial eccentric orbit around a black hole with a = 0.9M , p = 10M , and e = 0.1; for this orbit, Υ r = 2.824M . Note that e = 0.1 is fairly small, so we expect the expansion (4.27) to converge for a fairly small value of n max .
Solving Eq. (4.39) with n max = 3, we find 21 eigenvalues which we group as follows: . The eigenvalues corresponding to ∆n = ±n max fall slightly off of the trend given by Eq. (4.43). As stated above, we expect the shifted solutions to differ negligibly from one another when ∆n ≪ n max . The errors may be large when ∆n ≈ n max .
Next examine one of the eigenvectors. We focus on those corresponding to the eigenvalues labeled Υ 0,∆n s in Eq. (4.42). Because data at ∆n = ±n max appears to be less accurate than the other solutions, focus attention on the "innermost" 5 solutions (those corresponding to Υ 0 s = ±5.648M , ±2.824M , and 0). These five eigenvectors are Recall that S is a "vector of vectors," constructed by combining the Fourier components S j which are each associated with a given eigenvalue. These individual Fourier components are indicated in the eigenvectors shown in Eq. (4.44); five components are included. Notice the cyclic nature of the eigenvectors' components. This demonstrates how these components must be shifted according to Eq. (4.41) when we shift the eigenvalues. Notice also the slight numerical differences between the smallest components: Although they should be identical, we see both 7.806 × 10 −7 and 7.805 × 10 −7 . This difference goes away if we use a larger value of n max , although similar errors then appear in even smaller components. Empirically, we find that the eigenvectors corresponding to larger shifts tend to exhibit more numerical error. We thus use the set requiring no shift, for which the eigenvalues have the form Υ 1 s = −Υ −1 s ≡ Υ s , Υ 0 s = 0. Once these eigenvalues are identified, it is not difficult to isolate the associated eigenvectors, (S 1 , S 0 , S −1 ). We then extract the coordinate-domain spin eigenvector Fourier components S a j and assemble the solution for spin precession along an orbit: This sum runs from −j max to j max . For the example discussed above, the coordinate-domain eigenvector Fourier components are given by Notice that only the θ components of S 0 are non zero; if we had examined the S ±1 solutions, we would find by contrast that their θ components were all zero. Recall that for circular equatorial orbits, S θ is constant, and S r,φ undergo simple harmonic oscillation. For an eccentric equatorial orbit, we find a similar behavior: the components S r,φ are strongly coupled and oscillatory, whereas S θ varies independent of the other two components. We will revisit this behavior in Sec. V B. Figures 1 -4 show representative examples of how Υ s computed by this procedure varies with respect to orbital parameters for circular inclined and equatorial eccentric orbits. We omit the case of generic but resonant orbits since they are more complicated to compute, and do not add much to these results. Since the amount of data involved is significant, we do not show examples of the eigenvectors, though of course they are found by this procedure as well.

Figures 1 and 2 show how Υ s varies with inclination
θ inc for orbits of constant radius r for several values of black hole spin. The finite span of data used in Fig.  2 is simply because no stable orbits exist at r = 5M beyond some inclination θ max inc for the spins we include here. The gap near θ inc = 90 • in Fig. 1 is because we need a large value of k max in Eq. (4.26) to accurately model highly inclined orbits in the frequency domain. We have fixed k max = 12 for this initial analysis. Empirically, we find that the expansion converges to nine or ten digits of precision for k max = 12 on the range θ inc ≤ 70 • and θ inc ≥ 110 • . This holds for all values of r and a that we consider here; a larger value of k max would be needed if we examined inclinations θ inc outside this range. In other words, the value of k max needed for the expansion to converge depends quite a bit on orbital inclination, but only weakly on black hole spin and orbital radius.
One trend we see for the circular orbits is that Υ s does not vary by much with inclination angle at constant r, especially for large radius and small spin. This reflects the fact that the spacetime is nearly spherically symmetric for these orbits, and so Υ s ≃ √ rM at all inclinations. Figure 1 also shows a near symmetry between prograde and retrograde orbits: Υ s (θ inc ) ≃ Υ s (180 • − θ inc ). This also follows from the fact that the spacetime is nearly spherically symmetric at r = 10M . In all cases, we find that Υ s decreases from √ rM as θ inc increases toward 90 • (presumably reaching a minimum at θ inc = 90 • ), and increases back to √ rM as θ inc increases toward 180 • . Figures 3 and 4 show the variation of Υ s with eccentricity e for several examples of orbits of constant semi-latus rectum p, for several spins. The general trend we see is that Υ s monotonically increases from Υ s (e = 0) = √ pM when p is large. This is apparent in the left-hand (prograde) panel of Fig. 3, the a = 0.1M curve in the righthand (retrograde) panel of Fig. 3, and in the a = 0.9M and a = 0.8M curves of Fig. 4.
The behavior we find deviates from this tendency as orbits approach the last stable orbit (LSO), which marks the boundary between stable and unstable orbits. For prograde orbits at p = 5M , the LSO is at e LSO = 0.503 for a = 0.5M , and at e LSO = 0.796 for a = 0.6M . All prograde p = 5M orbits are stable for a = 0.7M , though the LSO is close to these orbits: p LSO → 4.79M as e → 1, not far beyond the large e portion of the a = 0.7M curve. The turnaround we see in the a = 0.9M and a = 0.5M curves in the right-hand panel of Fig. 3 appears to be the same phenomenon. For retrograde orbits at p = 10M , e LSO = 0.469 at a = 0.9M . All retrograde orbits at p = 10M are stable for a = 0.5M , though this curve lies close to the LSO: p LSO = 9.90M as e → 1.

C. Generic orbits
Finally, consider generic orbits, which include both radial and polar motions. We now need to solve Eq. (4.3) including both r and θ frequencies. To do so, we begin FIG. 1. Spin precession frequency Υs for circular inclined orbits of radius r = 10M as a function of θinc. The gap near θinc = 90 • is because we truncate the expansion (4.26) at kmax = 12. More terms than this must be kept in order for the precession matrix to be accurately represented in the frequency domain when the orbit is nearly polar. Two interesting features we see here are that Υs varies very little as a function of θinc, and that there is a near symmetry between the prograde (θinc < 90 • ) and retrograde (θinc > 90 • ) branches: Υs(θinc) ≃ Υs(180 • − θinc). These features are particularly pronounced for small a, and arise because the Kerr spacetime deviates from sphericity only slightly at r = 10M . by assuming that there exist solutions S a of the form (4.48) As in Sec. IV B, we find a common factor of e −iΥ a s λ which cancels. To go further, we use a trick adapted from classical FIG. 2. Spin precession frequency Υs for circular inclined orbits of radius r = 5M as a function of θinc. The maximum allowed inclination at this radius is θinc ≃ 59.9 • for a = 0.5M , and θinc ≃ 82.6 • for a = 0.9M ; no stable orbits exist at this radius for a = 0.1M . The variation of Υs with inclination is stronger at this radius than it was at r = 10M . This is not surprising, since this is a stronger field region of the spacetime, and the deviations from spherical symmetry are be much larger here. Aside from the stronger variation and the cutoffs associated with the lack of stable orbits at large θinc, the trends we see are qualitatively similar to those seen in Fig. 1.
mechanics [69] for analyzing biperiodic functions: on both sides, change Υ θ λ to Υ θ λ θ and Υ r λ to Υ r λ r ; multiply by e ipΥ θ λ θ e iqΥr λr ; integrate both λ θ and λ r over one full period. The result is Performing the sums over k ′ and n ′ , this simplifies to As in the previous analysis, truncate these sums at ±k max and ±n max , then expand the equation and rearrange. If written out directly, the rather cumbersome result would be difficult to analyze. We again clean things up substantially by nesting matrices, though we will need to take things one layer "deeper" since we have a second dynamical frequency in our problem. Let us first define (4.51) Next define the matrix P k whose elements are given by Recall that I is the 3 × 3 identity matrix.

FIG. 4. Spin precession frequency
Υs for prograde eccentric equatorial orbits with p = 5M as a function of eccentricity e. We see considerably more variation in the behavior of Υs here than we did for prograde orbits at p = 10M (left-hand panel of Fig. 3). The trend we find supports the idea that Υs decreases for orbits that come close to the last stable orbit: For prograde orbits at p = 5M eLSO ≃ 0.796 at a = 0.6M , and eLSO ≃ 0.503 at a = 0.5M . (For a = 0.7M , p = 5M is outside the separatrix between stable and unstable orbits for all e, though it is close to the separatrix as e → 1.) We then define and a matrix P whose elements are given by The matrix I is the (2n max + 1) × (2n max + 1) identity matrix.
As a concrete example, imagine that we truncate our sums at k max = 2, n max = 3. The vector S a k has 2n max +1 = 7 elements, and is given by (4.55) The matrix P k has (2n max + 1) 2 = 49 elements, and is given by (4.56) Next, the vector S a has 2k max + 1 = 5 elements: (4.57) The matrix P has (2k max + 1) 2 = 25 elements: Finally, to solve for the precessional motion of our spinning body, we must find the eigenvalues and eigenvectors of the matrix equation Doing so, we now find an even greater abundance of eigensolutions: Eq. (4.59) yields 3×(2k max +1)×(2n max + 1) eigenvalues and eigenvectors. This originates from essentially the same relabeling ambiguity that was responsible for the extra eigensolutions we found for Eq. (4.39), but with an additional order of "extraness" arising from the additional frequency associated with the underlying geodesic. Starting with Eq. (4.47), we can shift the sum over k ′ by ∆k ′ and the sum over n ′ by ∆n ′ , provided we shift the eigenvalue by and correspondingly shift the components of the eigenvectors.
Bearing this in mind, we find that the eigensolutions can be organized into three groups. As we discussed in Sec. IV B, the relabeling ambiguity means that the multiple solutions are in principle equivalent to one another, provided that the eigenvector indices are properly shifted. In practice, we find the most accurate solution corresponds to the eigenvector set Υ 1 s = −Υ −1 s ≡ Υ s , Υ 0 = 0. From the associated eigenvectors (S −1 , S 0 , S 1 ), we extract the Fourier components S kn , and finally assemble the solution for spin precession along a generic orbit: These sums run from −k max to k max , and from −n max to n max . Figure 5 shows an example of how Υ s varies with e and θ inc for a sample of generic orbits. We consider orbits with p = 10M about a black hole with spin a = 0.9M . We have fixed k max = 5 and n max = 5. The restriction on k means that we can only examine orbits with θ inc 60 • and θ inc 120 • . Although no issue of principle prevents us from increasing k max and examining more highly inclined orbits, we hold to k max = 5 in order to keep the precession matrix small and our analyses simple for this initial exploration.
The trends we see in Fig. 5 combine the trends we see for circular inclined and equatorial eccentric orbits. In particular, we see that Υ s tends to increase with e as other parameters are held fixed; it tends to decrease with θ inc with all other parameters fixed; and tends to decrease sharply as the LSO is approached. This last point is particularly visible in the right-hand panel of Fig.  5, which shows data for retrograde orbits (θ inc > 90 • ). Each trend line we plot is terminated near 10 the LSO. We see Υ s sharply decreasing in all of the cases we show here as we move to large e and approach the LSO.
One new feature of the general case is apparent in the left-hand panel of Fig. 5, which shows data for prograde orbits (θ inc < 90 • ): Υ s varies much more weakly with θ inc as e → 1. For e near 1, much of the orbit is spent at large radius, where the spacetime is nearly spherical, and Υ s does not depend on θ inc in the spherically symmetric limit. This behavior suggests that it may be useful to examine the limit e → 1 (for whichÊ G = 1). One might find a simple solution for Υ s in this case, much as we found Υ s = √ pM for circular and equatorial orbits. If such a solution exists, it may be a useful constraint for computing Υ s more generally.

V. THE SPIN-CURVATURE FORCE IN THE FREQUENCY DOMAIN
We now move to a discussion of the spin curvature force. We begin with some notation: we denote by f α forces defined with respect a trajectory's proper time τ . The force corresponding to the spin-curvature interaction is given by This is just a rewriting of Eq. (3.22). Since we find it very useful to use Mino time λ as our independent parameter, it is very useful to also define a force defined with respect to λ. Let us put Combining Eqs. (5.1) and (5.2), we have where U α G ≡ dx α /dλ = Σu α G . Using Eq. (3.25), which relates the spin tensor S αβ to the spin vector S α , plus the rule which follows from the constraint u G α S α = 0 and from Eq. (3.25), Eq. (5.3) can be written The components of the 4 × 3 curvature coupling matrix C are given by We now introduce the Fourier expansion. Each matrix element C α j can be expanded using Eqs. (2.29) and (2.30). We further know that the spin vector can be written as (4.61). Putting all of this together, the force components can be written (5.8) For presenting our results, it is useful to project F α S onto the Killing vectors and tensors: we define As the small body moves along its orbit, its geodesic energy, axial angular momentum, and Carter constant will each vary according to Eqs. (5.9)-(5.11) due to the spin-curvature force. These quantities can likewise be Fourier expanded: for C ∈ (E G , L G z , K G ), we have dλ jkn e −ijΥsλ e −i(kΥ θ +nΥr) . . We fixed kmax = 5 for this initial analysis. With this choice of kmax, we find empirically that we can accurately compute Υs for θinc 60 • and θinc 120 • . The behavior we see combines trends we see in both the equatorial eccentric and the inclined circular cases. For the prograde orbits, we see that Υs decreases with θinc at fixed e (similar to what we see in Fig. 1), but increases with e at fixed θinc (similar to what we see in the left-hand panel of Fig. 3). Interestingly, for prograde orbits Υs varies much less with θinc as e becomes large. For the retrograde cases, we see the strong influence of proximity to the last stable orbit (LSO), much as we saw in the right-hand panel of Fig. 3 and in Fig. 4.
With one exception, we use the expansion (5.12) of the components (5.9)-(5.11) in our detailed discussion of the spin force components that we study in Secs. V B and V C. The exception is for the radial component of the spin force, F r S . It is not difficult to see that F r S completely 11 decouples from Eqs. (5.9)-(5.11). Since it nonetheless includes important information, we examine it along with the components dE G /dλ, dL G z /dλ, and dK G /dλ. Before beginning our analysis, it is useful to examine how the quantities we will study scale with the small body's mass µ and the black hole's mass M . Begin with the spin components themselves. As already discussed, the spin magnitude S ≡ sµ 2 , with s ≤ 1 a dimensionless parameter. We also know that Consider now a point at some coordinate r = AM (where A is some positive number). Noting how the different 11 There are two ways to write the tensor K αβ . In one way, dK G /dλ is independent of F r S , but involves the other three components; the other way, it is independent of F θ S and involves the other three. In either case, one of the force components decouples from the projection of F α S onto the Kerr metric's Killing quantities.
metric components scale with black hole mass M at this coordinate, we infer that S r ∼ µ 2 , S θ,φ ∼ µ 2 /M . (5.14) Next examine the scaling of the various quantities which enter into the two forms of the spin-curvature force, Eqs. (5.1) and (5.3). From the behavior of the Riemann tensor, the components of the 4-velocity, the spin tensor, and the factor dτ /dλ, we see that Projecting onto the Killing vectors and Killing tensor to assemble the rates of change of geodesic energy, angular momentum, and Carter constant associated with the spin force, we find (5.17) In the figures that follow, we divide all of the quantities we plot by these scaling rules, so that one can easily assess the impact of our analysis for different masses and different small body spins.

A. Initial conditions for the spin vector
Our first step is to select initial conditions for the small body's spin vector. We begin by picking initial spatial components as measured in an orthonormal frame, and then convert the components to a Boyer-Lindquist coordinate frame. We emphasize that the orthonormal frame components are merely a convenient tool for visualizing the initial spin vector. An ensemble of spins that have the same magnitude S but differing orientations will have orthonormal frame components of similar magnitudes. By contrast, in the Boyer-Lindquist coordinate frame, the components' magnitudes will vary strongly as a function of the orbit's radius r and angle θ.
Let S be an orthonormal frame component, and let S j be a coordinate frame component. Then, (5.18) We must next fix the timelike spin component. To do so, we enforce S α u G α = 0, and find We used u G t = −Ê G . Finally, we have the rule that S ≡ √ S α S α = constant along the small body's worldline. With all this in mind, we use the following algorithm to set initial conditions on the spin vector: • Select initial components in the orthonormal frame.
• Increase or decrease all components by whatever factor is needed that so that S is some prescribed value.

B. Results
Here we show results for three representative cases. All are for orbits about a black hole with spin parameter a = 0.9M , but we consider different orbit geometries in order to explore how precession and the spin-curvature force behave in these different cases.
To present our results, we compare frequency-domain expansions for the evolution of the spin vectors [computed using Eq. (4.45)] and the spin-curvature force [computed with Eq. (5.12)] to time-domain solutions computed by directly integrating Eq. (4.1), and combining with Eq. (5.3). The time-domain solutions were found using the Mathematica function NDSolve, which adaptively chooses integration methods to optimize accuracy. As we'll discuss below and in Sec. V C, our timeand frequency-domain solutions agree with one another to within 10 −6 or 10 −7 over many hundred of orbital periods, an indication that the solutions are quite accurate. It is worth noting that in the frequency-domain approach, most of the computational effort is spent solving for the frequencies, the precession eigenvectors, and the Fourier-modes of the force components. Although we have not yet done a careful study of computational cost, our present studies indicate that the frequency-domain expansion is likely to be useful for large-scale studies of spin-enhanced binary motion, especially if it is practical to precompute and store such quantities in order for them to be read into a code and used as needed.

Equatorial orbits
Equatorial orbits have two particularly nice properties. First, the spin component Sθ = rS θ is constant along the background geodesic. To see this, examine the θ component of the precession equation (4.4): We used the equatorial condition (θ = π/2, u θ = 0) to simplify the general expression for P θ j . Use U r G = dr/dλ and rearrange: so Sθ = constant. The second useful property of these orbits is that K G = constant, at least at leading order in S. We initially saw this empirically, finding dK G /dλ = 0 in our numerics for all equatorial configurations. We then realized it is simple to prove this analytically, which we do in App. A.
The constancy of Sθ and K G means that only the components S r , S φ , dE G /dλ, and dL G z /dλ are interesting for equatorial orbits. Figure 6 shows the behavior of the spin components S r and S φ , as well as the force components dE G /dλ and dL G z /dλ, for the equatorial eccentric case we examine. Results for other equatorial cases are qualitatively similar. Notice that the frequency-domain solutions (solid lines) and time-domain solutions (dots) agree quite well. Detailed analysis shows that the two solutions differ by at most 10 −7 out to M λ = 500; see also discussion in Sec. V C. This level of agreement is typical for the equatorial eccentric cases we have examined, provided we include enough terms in the frequency-domain expansion. We go to n = 10 in the cases shown in Fig. 6.
Notice that the harmonic content of the solution is fairly complicated, especially for S φ . One noteworthy feature is a beat between the radial frequency and the precession frequency, with a period Λ beat = 2π/(Υ s − Υ r ) ≃ 12.85M −1 . We show enough data to see about two full cycles of this beat.

Circular and generic orbits
No simplification allows us to disregard components of the spin vector or the spin-curvature force for circular or generic orbits. Figure 7 shows the behavior of the three spin components S r,θ,φ (λ) and the three force components dE G /dλ, dL G z /dλ, and dK G /dλ for the inclined circular case. As in the equatorial case, we see excellent agreement between the frequency-and time-domain data in the plot. We find similar levels of agreement for other circular inclined cases we have examined, provided we include enough harmonics. For this plot, our frequency domain solution includes terms out to k max = 20. Detailed analysis for this case shows agreement within 10 −6 for small λ, drifting to a disagreement of ∼ 10 −5 at M λ = 500. This drift is due to errors in Υ s , and can be improved by including more terms in the Fourier expansion; see Sec. V C for further discussion of this point. We find that Υ s is computed more accurately for smaller values of θ inc .
As in the equatorial case, our solutions have rather ornate harmonic structure, with complicated beats between the polar frequency and the precession frequency. We show enough data to again capture about two full beat cycles, which have period Λ beat = 2π/(Υ θ − Υ s ) ≃ 8.49M −1 in this case.
Finally, Fig. 8 shows the spin and force components for a generic orbit. We include terms to k max = n max = 10 here. The time-and frequency-domain solutions differ by about 10 −6 for M λ 50, but the difference drifts to several × 10 −5 for M λ = 500. This can be improved by including more terms in the Fourier sums.
Because two orbital frequencies are available to beat against the precession frequency, the harmonic structure we find is particularly ornate for generic orbits. The span we show is sufficient to show beat structure at the radial-precession beat period [Λ beat−rp = 2π/(Υ s − Υ r ) ≃ 15.23M −1 ] and at the polar-precession beat pe-

The radial component
As discussed above, we examine the spin-curvature force component F r S separately since it decouples from dE G /dλ, dL G z /dλ, and dK G /dλ. Figure 9 shows this component for the three configurations that were used to produce Figs. 6, 7, and 8. The main new feature that we see here is that this force component has a non-zero average value: for these orbits, F r S < 0, indicating that there is an attractive force between the body and the black hole that it orbits. In retrospect, this is not surprising: on average, the spin vector in these three configurations has a component that is parallel to the black hole's spin, and it is well known that the spin-curvature interaction enhances the gravitational attaction of two bodies in a way that depends on the relative alignment of their spin and orbital angular momenta (see, e.g., discussion in Sec. VA6 of Ref. [89]). To clearly illustrate this, let us examine Eq. (5.3) in the limit of a circular, equatorial orbit, for an orbiting body whose spin is oriented normal to the equatorial plane. For this case, U r G = U θ G = 0, θ = π/2, and the only non-zero spin vector component is S θ . We find 23) where the upper sign is for prograde orbits and the lower is for retrograde. We have written this using the orthonormal form Sθ = √ g θθ S θ , which is particularly simple and convenient for a circular, equatorial orbit. Recall also that our F α is expressed using Mino time λ rather than proper time τ , and hence this differs by a factor Σ from other analyses of this quantity. Converting to d/dτ and considering r ≫ M , Eq. (5.23) agrees with the "spinorbit" and "spin-spin" forces given in Eqs. (5.69b) and (5.69c) of Ref. [89], which in turn follows the discussion and derivation given in Ref. [90].
Because F r S has an average non-zero value, it will play a role similar to the conservative self force, causing a shift to quantities like orbital frequencies, and leading to a secularly growing term in the orbital phase. Especially for orbits at small radius, this could have important implications for gravitational-wave source modeling.

C. Convergence
The Fourier expansions we have introduced, Eq. The agreement between the time-and frequency-domain solutions is outstanding. This is typical for equatorial orbits provided the frequency-domain expansion is truncated at a sufficiently large value of n. We used nmax = 10 here. (Note that S ≡ √ S α Sα is constant, although this is not apparent from the data shown here.) spin-curvature force, are formally correct only when an infinite number of terms are kept in the expansion. In this section, we examine how well these quantities converge when the sums are truncated at a finite number of terms. Figure 10 shows a typical example of how frequencydomain expansions converge as more terms are kept in the Fourier expansion. We show the convergence of the spin vector component S r as a function of Mino-time λ, for an equatorial eccentric orbit with parameters a = 0.9M , p = 5M , e = 0.7 (the orbit that was used to generate data for Fig. 6). We examine ∆S r as a function of λ, where with S r TD the result of direct time-domain integration of Eq. (4.1) along the geodesic, and S r FD, nmax from the frequency-domain expansion (4.45), truncating at n = n max .
To assess convergence, we examine this difference for n max = 5, 7, and 10. As expected, the difference decreases as n max is increased, but the manner in which it decreases is quite interesting. We generally see that, for small n max , ∆S r monotonically increases with λ, meaning that the two solutions drift away from each other. This drift is due to errors in our determination of the precession frequency, Υ s . Recall that Υ s is the eigenvalue of the precession equation (4.3). Determining this eigenvalue accurately requires us to describe the underlying geodesic motion accurately; in the frequency domain, this description becomes progressively more accurate as more terms are kept in our Fourier expansion. When we use n max = 5, our solution drifts by |∆S r | ≃ several × 10 −5 when we integrate to M λ = 500. Increasing this to n max = 7, the drift is reduced by about two orders of magnitude. Increasing still further to n max = 10 removes the drift altogether, and the two solutions differ by roughly 10 −7 at all times.
The convergence behavior shown for S r (λ) is typical: we see quite similar behavior for the other spin components and for all the components of the spin-curvature force. This is true for all the orbits we have examined, although the values of n max and k max needed to achieve convergence depends on an orbit's detailed parameters. We find it is not difficult to achieve a convergent frequency-domain representation of all physically important quantities relating to spin-curvature coupling for astrophysically interesting orbits.  Fig. 6, the solid (red) curves show quantities reconstructed from the frequency-domain expansion. Dots (green) show the same data computed by direct time-domain integration. We show the spin-vector components in the left-hand panels, and components of the spin-curvature force on the right. The lowest frequency component is a beat between the polar frequency Υ θ = 2.9230M and the precession frequency Υs = 2.1833M ; the λ range we show is wide enough that features in the corresponding period, 2πM/(Υ θ − Υs) ≃ 8.49, can be seen. We again find outstanding agreement between the timeand frequency-domain solutions provided that the frequency-domain expansion is truncated at sufficiently large k. We used kmax = 20 for the results shown here.

VI. ARE SPIN-PRECESSION RESONANCES INTERESTING?
Precession introduces a new frequency into the kinematics of an orbiting body. This suggests that interesting effects may arise when the spin frequency Υ s is commensurate with the orbital frequencies Υ θ,r , i.e. for orbits such that k and n can be found satisfying jΥ s + kΥ θ + nΥ r = 0 , (6.1) with j ∈ (−1, 0, 1). At such spin-orbit resonances, the Fourier mode corresponding to the harmonics which satisfy (6.1) will be constant with time, suggesting that the system's evolution could differ substantially from its behavior away from a resonance. Such behavior has been seen in studies of the dissipative self force [72,73] and of gravitational-wave recoil [91] when the Kerr orbital frequencies Ω r and Ω θ are commensurate. Although important for the self force problem and for radiation recoil, one can quickly convince oneself that spin-curvature coupling cannot produce dynamically important resonances, at least at linear order in S. Consider Eq. (5.12) for an orbit that has frequency harmonics which satisfy Eq. (6.1). For these modes, we would have dC/dλ = constant: these quantities would grow linearly in λ without bound. Such growth is inconsistent with the existence of the conserved integrals E S , L S z , and K S , defined by Eqs. (3.14), (3.15), and (3.28). The contributions to E S , L S z , and K S proportional to the spin tensor are oscillatory, so unbounded growth of the "geodesic" terms would, before long, violate their constant nature. To protect the system from these resonances, it must be the case that (dC/dλ) jkn = 0 for modes and frequencies which satisfy Eq. (6.1). Figure 11 shows an example of (dC/dλ) jkn vanishing on resonance. We examine a sequence of orbits, each with a = 0.9M , e = 0.2, and θ inc = 20 • , and varying the semi-latus rectum p. The sequence we examine includes one value of p such that Υ s −Υ θ +Υ r = 0. We indeed find that (dC/dλ) j=1,k=−1,n=1 is zero on resonance, though it is non-zero on nearby orbits. (We only show C = E G in this figure, although our data shows that this behavior also holds for L G z and K G .) Although there are no interesting effects due to these resonances at linear order, there might be interesting effects at higher order. We speculate that the behavior near resonance of nonlinear terms (which we have neglected) may play a role in pushing the dynamical evolution of spinning bodies from integrable to chaotic motion. The KAM theorem [92,93] teaches us that an integrable system will remain integrable under the influence of a weak, non-integrable perturbing force as long as that force is not resonant with the integrable motion. This suggests that circumstances may change when the perturbing force is in fact in resonance. If this speculation is correct, an analysis of spin-curvature forces near these resonances may help us to understand the onset of chaos that has been seen in previous analyses of spin-enhanced orbital motion [66][67][68].

VII. SUMMARY AND FUTURE WORK
We have shown how to compute the coupling between a small body's spin and the curvature of the black hole spacetime through which it moves in the frequency domain. The spin vector precesses along its orbit, which introduces a new frequency that must be accounted for in a frequency-domain description of quantities associated with this motion. Using the fact that Kerr geodesic motion is itself characterized by three frequencies, we have shown that all the quantities relevant to spin-curvature coupling can be computed very accurately.
Certain aspects of the large mass-ratio limit have proven to be useful and surprisingly accurate tools for helping to understand the two-body problem more generally. In particular, results that have come from the self force program have been shown to agree very well with re-sults computed using numerical relativity and have been used to refine effective one-body binary models [2][3][4]. We speculate that the results we present here may similarly find use by providing a limit that can be modeled precisely for understanding spin precessions and the influence of spin-curvature coupling in the evolution of binary systems. A useful starting point may be to check whether the solutions we find agree with the elegant results found by Gerosa et al. [35]. If so, our methods may prove to be a useful starting point for examining precession and spin-curvature coupling for more generic binaries than have been considered so far.
We plan to extend this work by using an osculating geodesic integrator [64,65] to develop spin-enhanced orbits. An osculating geodesic integrator models a nongeodesic worldline as a sequence of geodesic orbits. The spin-curvature force (5.1) then acts to move the small body from geodesic to geodesic in this sequence. As the orbit evolves, the geodesic energy E G will oscillate, but will do so in such a way that the corresponding spinenhanced energy E S remains fixed; similar statements hold for L G,S z and K G,S . The orbits which we find in this procedure will be very useful tools for allowing us to understand the importance of the spin on observable aspects of small body orbits very generally, allowing us to go beyond the special spin-orbit configurations that have been analyzed in detail in earlier work. In all three cases, this force component has an average value F r S < 0, corresponding to an attractive force between the orbiting body and the black hole, consistent with a weak-field analysis of orbits with aligned spins. This component of the spin-curvature force will play a role in the kinematics of these orbits akin to that played by the conservative self force.
Of particular interest will be to compare the spincurvature force to other non-geodesic effects that have an impact on a binary's evolution, such as the self force. We can get a very rough idea of how such forces compare by simply examining typical spin-curvature force components that we have computed (e.g., our Figs. 6 -9), along with typical self force components found by others (e.g., those shown in Figs. 6, 7, and 8 of Ref. [94]). It should be emphasized that such a comparison is extremely crude. Even after correcting for the factor of dτ /dλ = Σ between the two force definitions, and noting that we examine the forces on rather different orbits, we must be concerned about gauge. Best of all would be to develop a gaugeinvariant measure of the influence these forces have on binary orbits, as was done for simpler spin-orbit configurations in Refs. [22][23][24]. At present, the strongest defensible statement we can make is that the spin-curvature components appear large enough when crudely compared to self force components that it is very plausible that the small body's spin will leave an observationally important imprint on a binary's evolution, especially when we consider motion through the strong field of the binary's larger black hole.
Once spin-enhanced orbits are fully in hand, we can consider radiation from such configurations. As a first pass, it may not be too difficult to couple these orbits We show the difference between S r (λ) computed in the time domain using Eq. (4.1) and S r (λ) computed in the frequency domain using Eq. (4.45). We truncate the Fourier expansion at three different values of nmax in order to explore how well this expansion converges. The data shown here is for the same orbit we used in Fig. 6, i.e. an equatorial eccentric orbit with a = 0.9M , p = 5M , e = 0.7. A major source of error is due to inaccuracies in our determination of the precession frequency Υs This frequency is an eigenvalue of Eq. (4.3); we cannot solve the eigensystem accurately with insufficient Fourier modes. This accounts for the drift with time that we see in the top two panels: when nmax = 5, the two solutions drift by several × 10 −5 over an integration time M λ = 500. When nmax = 7, there is still a drift, but it is reduced by two orders of magnitude. When nmax = 10, the drift has been eliminated, and we find a difference |∆S r /µ 2 | 10 −7 at all times. The convergence of other spin components and the force components behaves in this way for all orbits we have examined, although the detailed values of nmax and kmax needed to converge depends on orbit parameters.
to some kind of "kludge" model for the evolution of orbit constants and wave emission. For example, a useful first approximation might be to develop a spin-enhanced worldline by combining the spin-curvature force with a self force, allowing us to compute an inspiral with both spin-precession and spin-coupling effects. We could then use that worldline as the source of a time-domain black hole perturbation theory solver (as has been done in, e.g., [95]), or even using a cruder approach based on some approximate set of radiative multipoles (as in e.g. Refs. [96,97]). Such a tool would likely be a useful first cut at building spin-enhanced waveforms to quantify the role that the small body's spin has the system's waves. The frequency-domain description may allow us to go beyond FIG. 11. Example of the vanishing of spin-curvature coupling on resonance. We examine a set of orbits about a black hole with a = 0.9M ; each orbit has eccentricity e = 0.3 and inclination θinc = 20 • . We vary the semi-latus rectum: p/M ∈ (2.78461, 2.83461, 2.88461, 2.98461). The case p/M = 2.83461 is resonant: Υs − Υ θ + Υr = 0 for that orbit. We find that the Fourier mode (dE G /dλ) j=1,k=−1,n=1 vanishes in the resonant case; the rates of change of L G z and K G likewise vanish on resonance for these mode indices. This behavior protects the conserved quantities E S , L S z , and K S , and guarantees that spin-orbit resonances aren't very interesting at least to linear order in the equations of motion and precession.
this and perhaps to extend black hole perturbation theory codes to include the influence of spin, much as is done in Ref. [71] but for general orbits and general spin orientations.