Rapidly rotating bodies moving in curved spacetime experience the socalled spincurvature force which becomes important for the motion of compact objects in gravitationalwave inspirals. As a first approximation this effect is captured in the motion of a spinning test particle. We solve the equations motion of a spinning particle to leading order in spin in arbitrary static and spherically symmetric spacetimes in terms of onedimensional closedform integrals. This solves the problem and proves its integrability in a wide range of modified gravities and near exotic compact objects. Then by specializing to the case of bound orbits in Schwarzschild spacetime we demonstrate how to express the solution in the form of Jacobi elliptic functions.
['General Relativity and Quantum Cosmology', 'High Energy Physics - Theory']
Analytic Solutions for the Motion of Spinning Particles near Spherically Symmetric Black Holes and Exotic Compact Objects
{'Analytic solutions for the motion of spinning particles near spherically symmetric black holes and exotic compact objects': "Vojtˇech Witzany 1, ∗ and Gabriel Andres Piovano 2, †  \n1 Institute of Theoretical Physics, Faculty of Mathematics and Physics, Charles University, CZ-180 00 Prague, Czech Republic 2 School of Mathematics and Statistics, University College Dublin, Belfield, Dublin 4, Ireland (Dated: September 10, 2024)  \nRapidly rotating bodies moving in curved space-time experience the so-called spin-curvature force, which becomes important for the motion of compact objects in gravitational-wave inspirals. As a first approximation, this effect is captured in the motion of a spinning test particle. We solve the equations motion of a spinning particle to leading order in spin in arbitrary static and spherically symmetric space-times in terms of one-dimensional closed-form integrals. This solves the problem and proves its integrability in a wide range of modified gravities and near exotic compact objects. Then, by specializing to the case of bound orbits in Schwarzschild space-time, we demonstrate how to express the solution in the form of Jacobi elliptic functions.  \nPACS numbers: 04.20.-q, 04.25.-g, 04.25.Nx, 04.30.Db, 04.70.-s  \nIntroduction. - After a track record of spectacular successes in their first three observing runs [1-3], in May 2023 the LIGO and Virgo instruments were joined for the first time by the Japanese detector KAGRA to start the fourth run of observing coalescing compact objects through gravitational waves [4]. The so-called LVK collaboration is expected to detect more events in this run than have been amassed to date. Nextgeneration detectors on Earth and in space promise to multiply the sensitivity and reach of the detectors even further [5, 6]. As the number of events grow, we will observe binaries with different mass ranges, mass ratios, and dynamical setups and, consequently, we will need a more faithful and complete picture of the two-body dynamics in order to detect and interpret these signals correctly.  \nIn particular, the space-based observatory Laser Interferometer Space Antenna (LISA) will detect inspirals of stellarmass compact objects (henceforth dubbed as secondary) onto supermassive black holes [7]. This will allow to accurately test whether the supermassive objects in the centres of galaxies are truly black holes or perhaps some exotic compact objects, or whether their gravitational field is described by Einstein gravity [8]. These so-called extreme-mass-ratio inspirals (EMRIs) are best modelled in a mass-ratio expansion [9]. At leading order, the two-body dynamics is approximated by the adiabatic inspiral of a test particle in the space-time of the massive compact object, where the inspiral is driven by fluxes sourced by the secondary [10, 11]. Subleading corrections to the motion are collectively called post-1-adiabatic, and include the effects of second-order in the mass-ratio fluxes, conservative self-force, and the spin of the smaller companion [12-14]. Post-1-adiabatic corrections are fundamental to model waveform suited for data analysis for LISA [7, 15], and the spin of the secondary plays a crucial role [16-18]. The effects of the latter are fully captured by the motion of a spinning test particle in the space-time of the massive compact object (computed to linear order to spin), and the outgoing gravitational-wave flux it sources [19-27]. Here we present an analytical solution of this motion near spherically symmet-  \nric compact objects in terms of quadratures. In particular, our analysis implies that the spinning particle motion is integrable at linear order in spin in any static, spherically symmetric space-time. We then specialize to Schwarzschild space-time and re-express the solution in terms of Jacobi elliptic integrals. Even though a number of works have treated this and similar topics previously (see, e.g., Refs [27-32]), this is the first time an analytic, closed-form solution of generic bound motion of a spinning test particle in Schwarzschild space-time is presented.  \nWe use the G = c = 1 geometrized units and the (-+++) signature of the metric. Greek indices label coordinate components. The Riemann tensor is defined by a ν ; κλ -a ν ; λκ ≡ R µ νκλ a µ , where ; denotes the covariant derivative and a µ is an arbitrary form.  \nSpinning particle. - The motion of a spinning particle is given by the equations of motion [33-35]  \nD P µ d τ = -1 2 R µ νκλ ˙ x ν S κλ , (1)  \nD S µν d τ = P µ ˙ x ν -P ν ˙ x µ , (2)  \nwhere P µ and ˙ x ν are the 4-momentum and tangent vector to the worldline of the spinning test-particle, while S µν = -S νµ is the spin tensor. We are only interested in the dynamics to linear order in spin S µν and, unless specified, all formulas are to be assumed at most O ( S ) accurate. This is justified by considering a comparably light compact object in the field of a heavy black hole, since the truncated terms can be shown to scale as higher order in the mass ratio, as discussed in the Supplemental material [36]. We fix the relation of the centroid x µ ( τ ) and the momentum P µ by the Tulczyjew-Dixon supplemental spin condition S µν P ν = 0 , to obtain P µ = m ˙ x µ + O ( S 2 ) where m is the particle mass [35, 37, 38]. The spin tensor can then be expressed as S µν = mε µνκλ ˙ x κ s λ / 2 , where s µ is the specific spin vector and ε µνκλ the Levi-Civita pseudo-tensor. The equations then  \nreduce to  \nD 2 x µ d τ 2 = -1 4 R µ νγδ ε γδ κλ ˙ x ν ˙ x κ s λ , (3)  \nD s λ d τ = 0 . (4)  \nMotion in static, spherically symmetric metric. - A general static spherically symmetric spacetime metric can be locally expressed in the form  \nd s 2 = -f ( r )d t 2 + h ( r )d r 2 + r 2 ( d θ 2 +sin θ 2 d ϕ 2 ) . (5)  \nThe Schwarzschild metric is contained within this class by setting f ( r ) , h ( r ) to f ( r ) = 1 /h ( r ) = 1 -2 M/r . The Killing vectors of this metric are generators of time translations and rotations around the x, y, z axes:  \nξ ( t ) = ∂ ∂t , (6)  \nξ ( x ) = -sin ϕ ∂ ∂θ -cos ϕ cot θ ∂ ∂ϕ , (7)  \nξ ( y ) = cos ϕ ∂ ∂θ -sin ϕ cot θ ∂ ∂ϕ , (8)  \nξ ( z ) = ∂ ∂ϕ . (9)  \nDixon [39] showed that for any Killing vector the spinning particle motion has a constant of the form C ( ξ ) = P µ ξ µ -ξ ρ ; σ S ρσ / 2 . In this case the Killing vectors correspond to a conserved specific spin-orbital energy and a formal angular momentum vector (we normalize each C ( ξ ) by m )  \nE = -f ˙ t + r 2 sin θf ' s ϕ ˙ θ -s θ ˙ ϕ 2 √ fh ,  \nJ x = -r 2 (sin ϕ θ +cos ϕ cos θ sin θ ϕ )  \n+ √ f h [ sin θ cos ϕh ( s t ˙ r -s r ˙ t ) + r sin ϕ sin θ ( s ϕ ˙ t -s t ˙ ϕ ) + r cos ϕ cos θ ( s t ˙ θ -s θ ˙ t ) ] ,  \nJ y = r 2 (cos ϕ ˙ θ -sin ϕ cos θ sin θ ˙ ϕ ) + √ f h [ sin θ sin ϕh ( s t ˙ r -s r ˙ t ) + r cos ϕ sin θ ( s t ˙ ϕ -s ϕ ˙ t ) + r sin ϕ cos θ ( s t ˙ θ -s θ ˙ t ) ] ,  \nJ z = r 2 sin 2 θ ϕ + √ f h [ cos θ h ( s t ˙ r -s r ˙ t ) + r sin θ ( s θ ˙ t -s t ˙ θ ) ] .  \n(10) ˙ ˙ (11) (12) ˙ (13)  \nThe linearisation in spin means that these integrals of motion are conserved up to O ( s 2 ) . Another integral of motion  \nis obtained by noticing that the following vector is paralleltransported along geodesics in general spherically symmetric static space-times:  \nl = r ˙ θ sin θ ∂ ∂ϕ -r sin θ ˙ ϕ ∂ ∂θ , (14)  \nand for spinning particles Dl µ / d τ = O ( s ) . As a result, the aligned component of the spin vector is an approximate constant of motion  \ns ∥ ≡ l µ s µ √ l ν l ν , d s ∥ d τ = 0 + O ( s 2 ) . (15)  \nAngular momentum aligned coordinates. - For general initial conditions, the formal angular-momentum vector ⃗ J = ( J x , J y , J z ) keeps pointing into a constant direction. Additionally, the orbital plane is always almost orthogonal to this vector up to O ( s ) corrections. Hence, without loss of generality, we can rotate into a new coordinate system θ, ϕ → ϑ, φ such that its axis of φ -rotation ( ϑ = 0 , π ) points in the direction of ⃗ J . When defined with respect to this new system, we obtain J ' x = J ' y = 0 , J ' z = √ J 2 x + J 2 y + J 2 z ≡ J , and the position of the particle will fulfil ϑ ( τ ) = π/ 2 + δϑ ( τ ) , where  \nδϑ = √ fh ( s r ˙ t -s t ˙ r ) r 2 ˙ φ + O ( s 2 ) . (16)  \nIn other words, the ϑ -motion of the particle is automatically expressed in terms of the other variables. This transformation was already implicitly used in the numerical studies of Refs [29, 40], even though the authors did not realise it also allows to analytically solve for the ϑ degree of freedom at O ( s ) . More details of the computation can be found in the Supplemental material.  \nThe equations of motion for the other orbital variables r ( τ ) , t ( τ ) , φ ( τ ) can now be expressed in first-order form up to O ( s 2 )  \n˙ r 2 = 1 h ( -1 + E 2 f -J 2 r 2 ) + s ∥ EJ (2 f -rf ' ) ( fh ) 3 / 2 r 2 , (17)  \n˙ t = E f + s ∥ J f ' 2 fr √ fh , (18)  \n˙ φ = J r 2 + s ∥ E r 2 √ fh . (19)  \nSolving parallel transport. - As for the spin degree of freedom, we need to solve the parallel transport equation at leading (geodesic) order. To do so, we construct a parallel-transported tetrad e µ ( A ) , A = 0 , ..., 3 inspired by Marck [41]. We start with the zeroth and third leg (component order t, r, φ, ϑ )  \ne µ (0) = ( E /f, ˙ r, J /r 2 , 0 ) , (20)  \ne µ (3) = ( 0 , 0 , 0 , -1 r ) , (21)  \nwhere ˙ r is given by eq. (17). Note that the full expression for ˙ r including O ( s ) corrections avoids singularities in the  \ndefinition of the tetrad at radial turning points. Now the first and second legs are given as e µ (1) = cos ψ ˜ e µ (1) +sin ψ ˜ e µ (2) and e µ (2) = -sin ψ ˜ e µ (1) +cos ψ ˜ e µ (2) , where  \n˜ e µ (1) = ( ˙ rr √ h √ f ( J 2 + r 2 ) , E r √ fh ( J 2 + r 2 ) , 0 , 0 ) , (22)  \n˜ e µ (2) = ( EJ f √ J 2 + r 2 , J ˙ r r , √ J 2 + r 2 r 2 , 0 ) , (23)  \nand the precession angle ψ ( τ ) is obtained by integrating  \n˙ ψ = EJ J 2 + r 2 . (24)  \nThe components of the spin with respect to the paralleltransported tetrad e µ ( A ) are constant at leading order. Additionally, since the spin magnitude s ≡ √ s µ s µ is conserved and s µ e µ (3) = s ∥ + O ( s 2 ) , the solution for the evolution of the spin vector can be expressed as  \ns t = √ s 2 -s 2 ∥ f ( J 2 + r 2 ) ( EJ cos ψ √ f + ˙ rr sin ψ ) ,  \n(25a) s r = √ s 2 -s 2 ∥ ( J ˙ r cos ψ r + E r sin ψ √ fh ( J 2 + r 2 ) ) , (25b) s φ = √ ( s 2 -s 2 ∥ )( J 2 + r 2 ) cos ψ r 2 , (25c) s ϑ = -s ∥ r . (25d)  \nOne can now also re-express δϑ from eq. (16) as  \nδϑ = √ ( s 2 -s 2 ∥ )( J 2 + r 2 ) sin ψ J r . (26)  \nSolution by quadrature. - The motion of the spinning test particle in any static, spherically symmetric metric is solved by quadrature as follows. First one integrates  \nτ ( r ) -τ ( r 0 ) = ± ∫ r r 0 d r ' √ R ( r ' ) , (27)  \nR ( r ) ≡ 1 h ( -1 + E 2 f -J 2 r 2 ) + s ∥ EJ (2 f -rf ' ) ( fh ) 3 / 2 r 2 , (28)  \nand inverts this relation to obtain r ( τ ) . This is then substituted into the r -parametrized solutions  \nt ( r ) -t ( r 0 ) = ± ∫ r r 0 d r ' √ R ( r ' ) [ E f + s ∥ J f ' 2 fr ' √ fh ] , (29)  \nφ ( r ) -φ ( r 0 ) = ± ∫ r r 0 d r ' √ R ( r ' ) [ J r ' 2 + s ∥ E r ' 2 √ fh ] , (30)  \nψ ( r ) -ψ ( r 0 ) = ± ∫ r r 0 d r ' √ R ( r ' ) [ EJ J 2 + r ' 2 ] . (31)  \nThe aligned coordinate position ϑ = π/ 2 + δϑ and the spin components s µ are then obtained by substituting the solutions for r ( τ ) and ψ ( τ ) into equations (25) and (26).  \nSchwarzschild space-time. - Let us demonstrate how to use this solution near a Schwarzschild black hole by setting f = 1 /h = 1 -2 M/r . We focus on bound motion between radial turning points r 1 > r 2 which are real roots of R ( r ) . Following Darwin's treatment of geodesics [42], we parametrize the motion by eccentricity e and semi-latus rectum p such that r 1 = p/ (1 -e ) , r 2 = p/ (1+ e ) . The relation between energy, angular momentum and the orbital elements e, p is obtained by inserting the geodesic relations from Darwin into the equation R ( r ) = 0 and computing the spin corrections  \nE 2 = ( p -2 M ) 2 -4 M 2 e 2 p [ p -M (3 + e 2 )] + s ∥ ( e 2 -1) 2 M √ Mp [ p 2 -4 Mp -4 M 2 ( e 2 -1)] p 2 [ p -M (3 + e 2 )] 2 , (32)  \nJ 2 = Mp p -M (3 + e 2 ) -s ∥ [ 2 p -3 M (3 + e 2 ) ] √ Mp [ p 2 -4 Mp -4 M 2 ( e 2 -1)] [ p -M (3 + e 2 )] 2  \n2 . (33)  \nInterestingly, these shifts exactly agree with those for a particle with spin fully aligned with orbital angular momentum when s ∥ = s (cf. Refs [43, 44]). The bound motion exists for e ∈ [0 , 1) and p ∈ ( p c ( e ) , ∞ ) , where the last stable orbits p c ( e ) are determined by the vanishing of the Jacobian of the transform in eqs (32) and (33), which yields (cf. the geodesic case in Ref. [45])  \np c = (6 + 2 e ) M +2 s ∥ √ 2(1 + e ) 3 + e . (34)  \nThis generalises the results for innermost stable circular orbits of particles with aligned spin from Refs [29, 30, 46, 47] to fully generic motion. The function R ( r ) can now be reexpressed as  \nR = 1 -E 2 r 4 ( r 1 -r )( r -r 2 )( r -r 3 )( r -r 4 ) , (35)  \nr 1 > r 2 > r 3 > r 4 , r 1 , 2 = p 1 ∓ e , r 4 = 0 , (36)  \nr 3 = 2 Mp p -4 M + 2 s ∥ √ Mp [ p 2 -4 Mp +4 M 2 (1 -e 2 )] ( p -4 M ) 2 . (37)  \nExpression as Jacobi elliptic integrals. - In Schwarzschild space-time, the most elegant parametrization of the motion is  \nthrough Carter-Mino time [48, 49] λ, d τ/ d λ = r 2 so that  \nλ ( r ) -λ ( r 2 ) = ∫ r r 2 d r ' √ R ( r ) = 2 F ( χ, k ) √ (1 -E 2 )( r 1 -r 3 ) r 2 , (38)  \nsin χ ≡ √ ( r 1 -r 3 )( r -r 2 ) ( r 1 -r 2 )( r -r 3 ) , (39)  \nk 2 = ( r 1 -r 2 ) r 3 ( r 1 -r 3 ) r 2 , (40)  \nR ( r ) ≡ (1 -E 2 )( r 1 -r )( r -r 2 )( r -r 3 ) r , (41)  \nwhere F ( χ, k ) is the elliptic integral of the first kind. All elliptic integrals F, K, E, Π and their inverses am , sn , cn will be defined in the angle-modulus convention of Byrd & Friedman [50] (see also Supplemental material). The λ ( r ) function can be inverted by the same arguments as given by van de Meent [51] for corresponding formulas for Kerr geodesics to yield  \nr ( λ ) = r 3 ( r 1 -r 2 )sn 2 ( K ( k ) π q r , k ) -r 2 ( r 1 -r 3 ) ( r 1 -r 2 )sn 2 ( K ( k ) π q r , k ) -( r 1 -r 3 ) , (42)  \nq r ≡ Υ r λ + q r 0 , (43)  \nπ  \n√  \n(1  \n-E  \n2  \nK  \n(  \nk  \n)  \n≡  \n,  \n(44)  \nwhere K ( k ) is the complete elliptic integral of the first kind and q r 0 is an integration constant determined by initial conditions. Again, following closely the notation and approach of Ref. [51] we express the other orbital variables as  \nt ( λ ) = q t + ˜ T r ( am ( q r π K ( k ) , k )) -˜ T r ( π ) 2 π q r , (45)  \nφ ( λ ) = q φ + ˜ Φ r ( am ( q r π K ( k ) , k )) -˜ Φ r ( π ) 2 π q r , (46)  \nψ ( λ ) = q ψ + ˜ Ψ r ( am ( q r π K ( k ) , k )) -˜ Ψ r ( π ) 2 π q r , (47)  \nq t ≡ Υ t λ + q t 0 , q φ ≡ Υ φ λ + q φ 0 , q ψ ≡ Υ ψ λ + q ψ 0 , (48)  \nwhere q t,φ,ψ 0 are again integration constants and  \n˜ T r ( χ ) = E √ (1 -E 2 ) r 2 ( r 1 -r 3 ) [ 2 M ( r 1 -r 3 )( r 2 -r 3 ) + r 3 ( r 3 ( r 2 + r 3 ) -r 1 ( r 2 -r 3 )) + 2 J E Ms ∥ r 3 -2 M F ( χ, k ) +( r 1 + r 2 + r 3 +4 M )( r 2 -r 3 )Π ( χ, r 1 -r 2 r 1 -r 3 , k ) + r 2 ( r 1 -r 3 ) E ( χ, k ) -2 M ( r 2 -r 3 ) ( 8 M 2 + J E s ∥ ) ( r 2 -2 M )( r 3 -2 M ) Π ( χ, ( r 3 -2 M )( r 1 -r 2 ) ( r 2 -2 M )( r 1 -r 3 ) , k ) -( r 1 -r 2 ) sin 2 χ √ r 2 ( r 1 -r 3 ) [ r 2 ( r 1 -r 3 ) -r 3 ( r 1 -r 2 ) sin 2 χ ] r 1 + r 2 -2 r 3 +( r 1 -r 2 ) cos 2 χ ] , (49)  \n˜ Φ r ( χ ) = 2( J + s ∥ E ) √ (1 -E 2 ) r 2 ( r 1 -r 3 ) F ( χ, k ) , (50)  \n˜ Ψ r ( χ ) = 2 EJ √ (1 -E 2 )( r 1 -r 3 ) r 2 ( J 2 + r 2 3 ) [ r 2 3 F ( χ, k ) + J 2 ( r 2 2 -r 2 3 ) J 2 + r 2 2 Re ( Π ( χ, r 3 -i J r 2 -i J r 1 -r 2 r 1 -r 3 , k )) -J ( r 2 -r 3 )( J 2 -r 2 r 3 ) J 2 + r 2 2 Im ( Π ( χ, r 3 -i J r 2 -i J r 1 -r 2 r 1 -r 3 , k )) ] , (51)  \n2  \n)(  \nr  \n-  \nr  \n)  \nr  \nr  \nΥ  \n1  \n3  \n2  \n[  \nΥ t = Υ r 2 π ˜ T r ( π ) = E 2 K ( k ) 2 M ( r 1 -r 3 )( r 2 -r 3 ) + r 3 ( r 3 ( r 2 + r 3 ) -r 1 ( r 2 -r 3 )) r 3 -2 M K ( k ) +( r 1 + r 2 + r 3 +4 M )( r 2 -r 3 )Π ( r 1 -r 2 r 1 -r 3 , k ) + r 2 ( r 1 -r 3 ) E ( k ) -( r 2 -r 3 )16 M 3 ( r 2 -2 M )( r 3 -2 M ) Π ( ( r 3 -2 M )( r 1 -r 2 ) ( r 2 -2 M )( r 1 -r 3 ) , k ) -J s ∥ E ( 1 -r 1 r 1 -2 M Π ( 2 M ( r 1 -r 2 ) r 2 ( r 1 -2 M ) , k )) ] , (52) Υ φ = Υ r 2 π ˜ Φ r ( π ) = J + s ∥ E , (53) Υ ψ = Υ r 2 π ˜ Ψ r ( π ) = EJ r 2 3 J 2 + r 2 3 + EJ 2 ( r 2 -r 3 ) K ( k )( J 2 + r 2 2 )( J 2 + r 2 3 ) [ J ( r 2 + r 3 ) Re ( Π ( r 3 -i J r 2 -i J r 1 -r 2 r 1 -r 3 , k )) -( J 2 -r 2 r 3 ) Im ( Π ( r 3 -i J r 2 -i J r 1 -r 2 r 1 -r 3 , k )) ] . (54)  \nFIG. 1. A generic inclined spinning particle orbit (left) and its description in angular-momentum-aligned coordinates (right) with p = 7 M,e = 0 . 3 , s ∥ = 0 , s = 0 . 3 M .  \n<!-- image -->  \nThis allows to plot fully explicit orbital solutions, as demonstrated in Fig. 1. It should be noted that Υ ψ , ˜ Ψ r ( χ ) can be evaluated at s ∥ = 0 at leading order (with no spin correction to r 3 , E , J ). On the other hand, Υ r,t,φ , ˜ T r , ˜ Φ r need to be explicitly expanded in s ∥ due to the spin corrections to r 3 , E , J at fixed p, e (see eq. (32),(33) and (37)). This procedure is straightforward, but the results are lengthy and thus relegated to the Supplemental material. From the ratios of these frequencies one can compute coordinate-time frequencies Ω r,φ,ψ ≡ Υ r,φ,ψ / Υ t or the nodal and periastron precession rates, which is also discussed in the Supplement.  \nDiscussion and outlooks. - The presented analytical solution for the motion of spinning test particles in Schwarzschild space-time can be taken 'as is' to source gravitational-wave fluxes or to compute the gravitational self-force on the particle [18-27]. Another useful output will be to extend the solution to compute action variables and scattering orbits in  \norder to compare to Effective one-body models [52-55] and other approaches to the relativistic two-body problem [56, 57]. The formulas provided in the notebooks along with this Letter can also be straightforwardly implemented into existing code-bases such as the BHPToolkit [58], thus providing an immediate impact in gravitational-wave modeling. In particular, this solution is a springboard to understand the radiationreaction on the fully generic motion of spinning particles in Kerr space-time [27, 32, 59]. Indeed, we have recently constructed semi-analytical orbits of spinning particles near Kerr black holes by using the Hamilton-Jacobi formalism presented in Ref. [59], and it is proving important to have independently derived analytical formulas for validation, even if just for the Schwarzschild case. An important application will also be to use the quadratures in equations (27)-(31) to study deviations of the particle motion in modified gravities such as in Einstein-Gauss-Bonnet gravity [60, 61], Horndeski and beyond-Horndeski theories [62, 63], or near boson stars [6467] and other exotic compact objects [68].  \nHowever, we can already use our results to draw general conclusions about the detectability of secondary spin in largemass-ratio inspirals onto spherically symmetric compact objects. Our solution demonstrates that time-averaged observables depend only on the aligned component of spin s ∥ , independently of the compact-object model, and the contributions of the orthogonal component s ⊥ , s 2 ⊥ ≡ s 2 -s 2 ∥ to any dynamical variable oscillate with the characteristic frequency Ω ψ generically different from the orbital frequencies. The spin contributions to the post-1-adiabatic phase of the waveform depends only the spin shifts to quantities such as average orbital frequencies and the average energy and angular momentum fluxes, which are independent of s ⊥ (see also Ref. [24]). Furthermore, the Ω ψ -oscillation due to s ⊥ ≲ m only appears in the O ( m/M ) sub-dominant corrections to the waveform amplitude [69, 70]. As a result, waveform templates that neglect s ⊥ contributions will have a mismatch with models that include s ⊥ of order O ( m 2 /M 2 ) . Hence, according to the Lindblom criterion, the waveforms that either include or ne-  \ne indistinguishable from each other in a matched filtering analysis unless the signal to noise ratio of the EMRI signal reaches the order of the large mass ratio M/m ∼ 10 4 -10 6 [71, 72]. This provides a key insight that will vastly simplify the treatment of spin in modified-gravity waveforms at large mass ratios.  \nAnother consequence of this work is that the motion of spinning particles in any static spherically symmetric spacetime is integrable to linear order in spin. In other words, if there is any chaos to be observed by the numerical integration of spinning particle motion [29, 31], it is an O ( s 2 ) effect that goes beyond the validity of the equations itself, which is in lines with the quantitative numerical scalings of resonances observed in Schwarzschild space-time [31].  \nAcknowledgements. - We would like to thank Viktor Skoup'y and Josh Mathews for pointing out typos in a draft of the paper. This work make use of the Wolfram Mathematica package 'KerrGeodesics' [73], which is a part of the BHPToolkit [58]. VW is grateful for the support from the Charles University Primus programme. GAP aknowledges support from an Irish Research Council Fellowship under grant number GOIPD/2022/496.  \n- ∗ vojtech.witzany@matfyz.cuni.cz † gabriel.piovano@ucd.ie\n- [1] B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. X 9 , 031040 (2019), arXiv:1811.12907 [astro-ph.HE].\n- [2] R. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. X 11 , 021053 (2021), arXiv:2010.14527 [gr-qc].\n- [3] R. Abbott et al. (KAGRA, VIRGO, LIGO Scientific), Phys. Rev. X 13 , 041039 (2023), arXiv:2111.03606 [gr-qc].\n- [4] B. P. Abbott et al. (KAGRA, LIGO Scientific, Virgo, VIRGO), Living Rev. Rel. 21 , 3 (2018), arXiv:1304.0670 [gr-qc].\n- [5] V. Kalogera et al. , (2021), arXiv:2111.06990 [gr-qc].\n- [6] LISA Consortium, (2017), arXiv:1702.00786 [astro-ph].\n- [7] S. Babak, J. Gair, A. Sesana, E. Barausse, C. F. Sopuerta, C. P. L. Berry, E. Berti, P. Amaro-Seoane, A. Petiteau, and A. Klein, Phys. Rev. D95 , 103012 (2017), arXiv:1703.09722 [gr-qc].\n- [8] K. G. Arun et al. (LISA), Living Rev. Rel. 25 , 4 (2022), arXiv:2205.01597 [gr-qc].\n- [9] L. Barack and A. Pound, Rept. Prog. Phys. 82 , 016904 (2019), arXiv:1805.10385 [gr-qc].\n- [10] S. A. Hughes, N. Warburton, G. Khanna, A. J. K. Chua, and M. L. Katz, Phys. Rev. D 103 , 104014 (2021), [Erratum: Phys.Rev.D 107, 089901 (2023)], arXiv:2102.02713 [gr-qc].\n- [11] S. A. Hughes, S. Drasco, E. E. Flanagan, and J. Franklin, Phys. Rev. Lett. 94 , 221101 (2005), arXiv:gr-qc/0504015.\n- [12] M. van de Meent, Phys. Rev. D 97 , 104033 (2018), arXiv:1711.09607 [gr-qc].\n- [13] N. Warburton, A. Pound, B. Wardell, J. Miller, and L. Durkan, Phys. Rev. Lett. 127 , 151102 (2021), arXiv:2107.01298 [gr-qc].\n- [14] B. Wardell, A. Pound, N. Warburton, J. Miller, L. Durkan, and A. Le Tiec, Phys. Rev. Lett. 130 , 241402 (2023), arXiv:2112.12265 [gr-qc].\n- [15] L. Barack et al. , Class. Quant. Grav. 36 , 143001 (2019), arXiv:1806.05195 [gr-qc].\n- [16] E. A. Huerta, J. R. Gair, and D. A. Brown, Phys. Rev. D 85 , 064023 (2012), arXiv:1111.3243 [gr-qc].\n- [17] E. A. Huerta and J. R. Gair, Phys. Rev. D 84 , 064023 (2011), arXiv:1105.3567 [gr-qc].\n- [18] G. A. Piovano, R. Brito, A. Maselli, and P. Pani, Phys. Rev. D 104 , 124019 (2021), arXiv:2105.07083 [gr-qc].\n- [19] N. Warburton, T. Osburn, and C. R. Evans, Phys. Rev. D 96 , 084057 (2017), arXiv:1708.03720 [gr-qc].\n- [20] S. Akcay, S. R. Dolan, C. Kavanagh, J. Moxon, N. Warburton, and B. Wardell, Phys. Rev. D 102 , 064013 (2020), arXiv:1912.09461 [gr-qc].\n- [21] J. Mathews, A. Pound, and B. Wardell, Phys. Rev. D 105 , 084031 (2022), arXiv:2112.13069 [gr-qc].\n- [22] V. Skoup'y and G. Lukes-Gerakopoulos, Phys. Rev. D 103 , 104045 (2021), arXiv:2102.04819 [gr-qc].\n- [23] V. Skoup'y and G. Lukes-Gerakopoulos, Phys. Rev. D 105 , 084033 (2022), arXiv:2201.07044 [gr-qc].\n- [24] V. Skoup'y, G. Lukes-Gerakopoulos, L. V. Drummond, and S. A. Hughes, Phys. Rev. D 108 , 044041 (2023), arXiv:2303.16798 [gr-qc].\n- [25] G. A. Piovano, A. Maselli, and P. Pani, Phys. Rev. D 102 , 024041 (2020), arXiv:2004.02654 [gr-qc].\n- [26] G. A. Piovano, A. Maselli, and P. Pani, Phys. Rev. D 107 , 024021 (2023), arXiv:2207.07452 [gr-qc].\n- [27] L. V. Drummond, P. Lynch, A. G. Hanselman, D. R. Becker, and S. A. Hughes, (2023), arXiv:2310.08438 [gr-qc].\n- [28] R. H. Rietdijk and J. W. van Holten, Class. Quant. Grav. 10 , 575 (1993).\n- [29] S. Suzuki and K.-i. Maeda, Phys. Rev. D 58 , 023005 (1998), arXiv:gr-qc/9712095.\n- [30] P. I. Jefremov, O. Y. Tsupko, and G. S. Bisnovatyi-Kogan, Phys. Rev. D 91 , 124030 (2015), arXiv:1503.07060 [gr-qc].\n- [31] O. Zelenka, G. Lukes-Gerakopoulos, V. Witzany, and O. Kop'aˇcek, Phys. Rev. D 101 , 024037 (2020), arXiv:1911.00414 [gr-qc].\n- [32] L. V. Drummond and S. A. Hughes, Phys. Rev. D 105 , 124041 (2022), arXiv:2201.13335 [gr-qc].\n- [33] M. Mathisson, Acta Phys. Polon. 6 , 163 (1937).\n- [34] A. Papapetrou, Proc. Roy. Soc. Lond. A209 , 248 (1951).\n- [35] W. Dixon, Il Nuovo Cimento (1955-1965) 34 , 317 (1964).\n- [36] See Supplemental Material at http://link.aps.org/ supplemental/ PLEASE-REPLACE , which includes additional Refs. [74-83], for further details of the derivations in this paper.\n- [37] W. Tulczyjew, Acta Phys. Pol. 18 , 393 (1959).  \n[38] O. Semer'ak, Mon. Not. R. Astron. Soc.  \n308  \n, 863 (1999).  \n- [39] W. G. Dixon, Proc. Roy. Soc. Lond. A314 , 499 (1970).\n- [40] O. Zelenka, G. Lukes-Gerakopoulos, V. Witzany, and O. Kop'aˇcek, (2019), arXiv:1911.00414 [gr-qc].\n- [41] J.-A. Marck, Proc. R. Soc. Lond. A 385 , 431 (1983).\n- [42] C. G. Darwin, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 249 , 180 (1959).\n- [43] L. V. Drummond and S. A. Hughes, Phys. Rev. D 105 , 124040 (2022), arXiv:2201.13334 [gr-qc].\n- [44] S. Mukherjee and S. Tripathy, Phys. Rev. D 101 , 124047 (2020), arXiv:1905.04061 [gr-qc].\n- [45] C. Cutler, D. Kennefick, and E. Poisson, Phys. Rev. D 50 , 3816 (1994).\n- [46] M. Favata, Phys. Rev. D 83 , 024028 (2011), arXiv:1010.2553 [gr-qc].\n- [47] O. Y. Tsupko, G. Bisnovatyi-Kogan, and P. Jefremov, Gravitation and Cosmology 22 , 138 (2016).\n- [48] B. Carter, Phys. Rev. 174 , 1559 (1968).\n- [49] Y. Mino, Phys. Rev. D 67 , 084027 (2003).  \n- [50] P. F. Byrd and M. D. Friedman, Handbook of elliptic integrals for engineers and physicists , Vol. 67 (Springer, 2013).\n- [51] M. van de Meent, Class. Quant. Grav. 37 , 145007 (2020), arXiv:1906.05090 [gr-qc].\n- [52] A. Buonanno and T. Damour, Phys. Rev. D 59 , 084006 (1999), arXiv:gr-qc/9811091.\n- [53] S. Akcay, R. Gamba, and S. Bernuzzi, Phys. Rev. D 103 , 024014 (2021), arXiv:2005.05338 [gr-qc].\n- [54] A. Albertini, R. Gamba, A. Nagar, and S. Bernuzzi, Phys. Rev. D 109 , 044022 (2024), arXiv:2310.13578 [gr-qc].\n- [55] M. Khalil, A. Buonanno, H. Estelles, D. P. Mihaylov, S. Ossokine, L. Pompili, and A. Ramos-Buades, (2023), arXiv:2303.18143 [gr-qc].\n- [56] L. Blanchet, A. Buonanno, and A. Le Tiec, Phys. Rev. D 87 , 024030 (2013), arXiv:1211.1060 [gr-qc].\n- [57] Z. Bern, A. Luna, R. Roiban, C.-H. Shen, and M. Zeng, Phys. Rev. D 104 , 065014 (2021), arXiv:2005.03071 [hep-th].\n- [58] 'Black Hole Perturbation Toolkit,' (bhptoolkit.org) ().\n- [59] V. Witzany, Phys. Rev. D 100 , 104030 (2019), arXiv:1903.03651 [gr-qc].\n- [60] D. Glavan and C. Lin, Phys. Rev. Lett. 124 , 081301 (2020), arXiv:1905.03601 [gr-qc].\n- [61] Y.-P. Zhang, S.-W. Wei, and Y.-X. Liu, Universe 6 , 103 (2020), arXiv:2003.10960 [gr-qc].  \n[62] E. Babichev, C. Charmousis,  \nand A. Leh'ebel, JCAP  \n04  \n, 027  \n(2017), arXiv:1702.01938 [gr-qc].  \n- [63] A. Bakopoulos, C. Charmousis, P. Kanti, and N. Lecoeur, JHEP 08 , 055 (2022), arXiv:2203.14595 [gr-qc].\n- [64] M. Kesden, J. Gair, and M. Kamionkowski, Phys. Rev. D 71 , 044015 (2005), arXiv:astro-ph/0411478.\n- [65] C. F. B. Macedo, P. Pani, V. Cardoso, and L. C. B. Crispino, Phys. Rev. D 88 , 064046 (2013), arXiv:1307.4812 [gr-qc].\n- [66] V. Diemer, K. Eilers, B. Hartmann, I. Schaffer, and C. Toma, Phys. Rev. D 88 , 044025 (2013), arXiv:1304.5646 [gr-qc].\n- [67] Y. Brihaye, V. Diemer, and B. Hartmann, Phys. Rev. D 89 , 084048 (2014), arXiv:1402.1055 [gr-qc].\n- [68] V. Cardoso and P. Pani, Living Rev. Rel. 22 , 4 (2019), arXiv:1904.05363 [gr-qc].\n- [69] E. E. Flanagan and S. A. Hughes, Phys. Rev. D 57 , 4535 (1998), arXiv:gr-qc/9701039.\n- [70] A. Pound and B. Wardell, (2021), arXiv:2101.04592 [gr-qc].\n- [71] L. Lindblom, B. J. Owen, and D. A. Brown, Phys. Rev. D 78 , 124020 (2008), arXiv:0809.3844 [gr-qc].\n- [72] B. F. Schutz, in The detection of gravitational waves (Cambridge University Press, 1991) pp. 406-451.\n- [73] 'KerrGeodesics,' (doi.org/10.5281/zenodo.8108265) ().\n- [74] I. Bailey and W. Israel, Commun. Math. Phys. 42 , 65 (1975).\n- [75] J. Ehlers and E. Rudolph, General Relativity and Gravitation 8 , 197 (1977).\n- [76] T. Binnington and E. Poisson, Phys. Rev. D 80 , 084018 (2009), arXiv:0906.1366 [gr-qc].\n- [77] J. Steinhoff and D. Puetzfeld, Phys. Rev. D86 , 044033 (2012), arXiv:1205.3926 [gr-qc].\n- [78] J. Steinhoff, T. Hinderer, A. Buonanno, and A. Taracchini, Physical Review D 94 , 104028 (2016).\n- [79] J. Steinhoff and G. Schafer, EPL 87 , 50004 (2009), arXiv:0907.1967 [gr-qc].\n- [80] S. Chakrabarti, T. Delsate, N. Gurlebeck, and J. Steinhoff, Phys. Rev. Lett. 112 , 201102 (2014), arXiv:1311.6509 [gr-qc].\n- [81] C. Møller, Communications of the Dublin Institute for Advanced Studies (1949).\n- [82] A. Morbidelli, Modern celestial mechanics: aspects of solar system dynamics (2002).\n- [83] R. Fujita and W. Hikida, Class. Quant. Grav. 26 , 135002 (2009).", 'Truncation of equations of motion at O ( s )': "The Mathisson-Papapetrou-Dixon equations of motion valid at quadrupolar order in a given background are [33-35, 39]  \nD P µ d τ = -1 2 R µ νκλ ˙ x ν S κλ + 1 6 R νκλγ ; µ J νκλγ + ... , (S1)  \nD S µν d τ = P µ ˙ x ν -P ν ˙ x µ -4 3 ( R µ κλγ J νκλγ -R ν κλγ J µκλγ ) + ... , (S2)  \nwhere J µνκλ is the Dixon's quadrupole, which encodes the deformations of an extended body due to external tidal fields or its rotation [see, e.g. 74-78]. The three dots denote octupole and higher order multipoles which are not explicitly discussed here for the sake of simplicity. The pole-dipole-quadrupole equations of motion are no longer universal because the quadrupole depends on the matter dynamics of the considered body. Adiabatic and dynamical tidal quadrupoles [76, 78] turn out to be suppressed by the fifth power of the mass ratio in the large-mass ratio expansion, so we will not discuss them here. However, J µνκλ also includes a contribution to the quadrupole sourced by the rotation of the body. This term turns out to be second-order-in-massratio as well as O ( S 2 ) [79]:  \nJ µνκλ = 3 m 2 M 3 c ES 2 P [ µ S ν ] γ S γ [ κ P λ ] , (S3)  \nwhere [] denotes antisymmetrization of indices, M 2 ≡ -P µ P µ , m ≡ -P µ ˙ x µ , and c ES 2 is a deformation coefficient that can be assumed to be approximately constant and ∝ m 3 almost up to the mass-shedding limit [see, e.g., 80]. The curvature is R µ νκλ ∼ 1 /r 2 c , where r c is the curvature length, and a covariant derivative adds an additional ∼ 1 /r c factor. The Christoffelsymbol (parallel transport/geodesic) terms in the equations of motion scale as ∼ 1 /r c . The dipole terms thus scale with a relative factor ∼ S/r c as compared to the leading order in the equations of motion, and the quadrupole terms as ∼ S 2 /r 2 c relative to leading order. Similarly, higher multipoles contain higher powers of S can be shown to scale with increasing powers of the ratio of spin and the curvature length.  \nThe linear momentum is related to the four-velocity in the following manner under the Tulczyjew-Dixon condition S µν P ν = 0 [77]:  \n˙ x µ = U µ + 2 S µν S κλ R κλνγ U γ 4 M 2 + S γδ S ωχ R γδωχ , (S4)  \nU µ ≡ m M 2 P µ -4 3 M 2 R [ µ νκλ J γ ] νκλ P γ + 1 6 M 2 S µν R κλχξ ; ν J κλχξ . (S5)  \nOne of the consequences of relation above is that M 2 ˙ x µ = mP µ + O ( S 2 ) so by projection into ˙ x µ we obtain M 2 = m 2 + O ( S 2 ) and both definitions of mass agree at leading order in spin. By taking a covariant time derivative of the relation above and of s µ ≡ ε µνκλ P ν S κλ / (2 M 2 ) we then obtain  \nD 2 x µ d τ 2 = -1 4 R µ νγδ ε γδ κλ ˙ x ν ˙ x κ s λ + O ( s 2 R ∇ R ) , (S6)  \nD s λ d τ = 0 + O ( s 2 R 2 ) . (S7)  \nIn particular, the terms like s 2 R ∇ R and s 2 R 2 have comparable contribution from the quadrupole ∝ c ES 2 coefficient as well from the non-trivial momentum-velocity relation in the Tulczyjew-Dixon supplementary spin condition. Nevertheless, all terms including the Riemann tensor and its derivatives scaled with the curvature length r c in the same way. The Riemann tensor scales as 1 /r 2 c while the magnitude of the specific spin s is set by the Møller radius r M , the minimal physical size of the rotating body [81]. Thus, for compact objects in the large mass ratio limit, the truncation of the equations of motion (3) and (4) are controlled by powers of r M /r c . A multipole expansion is valid only when r M /r c ≪ 1 , a condition that is always satisfied for extreme mass-ratio binaries. This is because for bound motion near a black hole r c ∼ √ r 3 /M ≳ M and thus also 1 /r c ≲ 1 /M , and the Møller radius r M ≲ m for compact objects. As such, the equations of motion (3) and (4) presented in the main text corresponds to the general form of the Mathisson-Papapetrou-Dixon equations of motion truncated at leading order in an expansion in m/M . Before adding the O ( S 2 ) terms we would thus have to add leading-order gravitational self-force for self-consistence; the crossinteractions between self-force and spin would then be of the same importance as the O ( S 2 ) terms in the equations of motion, which is beyond the scope of our Letter.", 'Aligning the coordinate frame': "Given the components of the angular-momentum vector ⃗ J computed from equations (11)-(13) in the main text and some initial conditions, the inclination angle ι ∈ (0 , π ) of ⃗ J with respect to the z direction is given by  \ncos ι = ⃗ J · ˆ z J = J z J , (S8)  \nJ ≡ √ J 2 x + J 2 y + J 2 z , (S9)  \nwhere ˆ z is the unit vector in the z direction. Now the aligned frame is obtained by rotating by ι around the axis ⃗n given as  \n⃗n = ⃗ J × ˆ z | ⃗ J × ˆ z | = ( -cos ξ, sin ξ, 0) , (S10)  \nξ ≡ arctan J x J y , (S11)  \nwhere ξ is related to the longitude of the ascending node known from classical celestial mechanics [82]. The angles ι, ξ can be used to parametrize the angular-momentum vector as  \n⃗ J = J (sin ξ sin ι, cos ξ sin ι, cos ι ) . (S12)  \nThe transformation of any vector into the aligned coordinates then reads  \n⃗v ' = ⃗n ( ⃗n · ⃗v ) + cos ι ( ⃗n × ⃗v ) × ⃗n -sin ι ( ⃗n × ⃗v ) . (S13)  \nIt is now easy to show that this map yields ⃗ J ' = (0 , 0 , J z ' = J ) . Note that the transformation also flips counter-rotating coordinate systems ( J z < 0 ) into co-rotating ones ( J z ' > 0 ), which corresponds to ι > π/ 2 . The inverse map is simply obtained by switching the sign on the angle of rotation ι  \n⃗v = ⃗n ( ⃗n · ⃗v ' ) + cos ι ( ⃗n × ⃗v ' ) × ⃗n +sin ι ( ⃗n × ⃗v ' ) . (S14)  \nFrom the transformation of the unit position vector ˆ r = (sin ϕ sin θ, cos ϕ sin θ, cos θ ) and the aligned vector ˆ r ' = (sin φ sin ϑ, cos φ sin ϑ, cos ϑ ) one can easily read off the aligned coordinate system and its inverse transformation. Specifically, one can take the t, r, φ, ϑ = π/ 2 + δϑ solution presented in the main text of the paper and map it back into the original coordinates ( t, r are unchanged)  \ncos θ = cos ϑ cos ι -sin ϑ sin ι cos( ξ -φ ) , (S15)  \ntan ϕ = cos ϑ sin ι sin ξ +sin ϑ [ cos 2 ι 2 sin φ -sin 2 ι 2 sin(2 ξ -φ ) ] cos ϑ sin ι cos ξ +sin ϑ [ cos 2 ι 2 cos φ -sin 2 ι 2 cos(2 ξ -φ ) ] . (S16)", 'Expressing δϑ': 'Now let us express ϑ ( τ ) = π/ 2 + δϑ ( τ ) from J x = J y = 0 . First, let us examine the J x = 0 expression  \n0 = cos φ [ r 2 δϑ ˙ φ + fh ( s t ˙ r -s r ˙ t ) √ fh ] +sin φ [ -r 2 δ ˙ ϑ + f ( s φ ˙ t -s t ˙ φ ) √ fh ] . (S17)  \nThe J y expression is identical up to sin φ ↔ cos φ . This means that the expressions in the square brackets in eq. (S17) must vanish independently. This may seem like an overconstrained system since both δϑ and δ ˙ ϑ fulfill apparently independent conditions. However, one can see that by taking the expression for δϑ from eq. (16), taking a time-derivative, and substituting equations of motion, the two equations yield the same information and are mutually consistent.  \nOnce we have expressed δϑ using the solution for r, ψ and the solution for φ , all presented in the main text (cf. eqs (26)-(31)), we obtain the explicit evolution of the generic coordinates θ, ϕ by linearizing the formulas (S15) and (S16) to obtain  \nθ = π -arccos [cos( ξ -φ ) sin ι ] + δϑ cos ι √ 1 -cos 2 ( ξ -φ ) sin 2 ι , (S18)  \nϕ = arctan [ cos 2 ι 2 sin φ -sin 2 ι 2 sin(2 ξ -φ ) cos 2 ι 2 cos φ -sin 2 ι 2 cos(2 ξ -φ ) ] -4 δϑ sin ι sin( ξ -φ ) 3 + cos 2 ι -2 cos(2 ξ -2 φ ) sin 2 ι . (S19)', 'Jacobi elliptic integrals': 'The convention of Byrd & Friedman [50] for incomplete Jacobi elliptic integrals of the first, second and third kinds F, E, Π is as follows  \nF ( χ, k ) ≡ ∫ χ 0 d ξ √ 1 -k 2 sin 2 ξ , (S20)  \nE ( χ, k ) ≡ ∫ χ 0 √ 1 -k 2 sin 2 ξ d ξ , (S21)  \nΠ( χ, n, k ) ≡ ∫ χ 0 d ξ (1 -n sin 2 ξ ) √ 1 -k 2 sin 2 ξ . (S22)  \n(S23)  \nThe complete elliptic integrals are then defined by K ( k ) ≡ F ( π/ 2 , k ) , E ( k ) ≡ E ( π/ 2 , k ) , Π( n, k ) ≡ Π( π/ 2 , n, k ) . The Jacobi amplitude is the inverse of the elliptic integral of the first kind F (am( u, k ) , k ) ≡ u and we also use sn( u, k ) ≡ sin(am( u, k )) , cn( u, k ) ≡ cos(am( u, k )) .  \nThis convention should be contrasted with the convention of Maplesoft Maple , which uses sin( ξ ) in the first argument of incomplete integrals, or Wolfram Mathematica , which uses k 2 in the last argument and switches the order of the first two arguments in the incomplete elliptic integral of the third kind Π .', 'Linear shift to the orbital frequencies': 'We write the linear correction to the Mino frequencies as  \nΥ r,t,φ = Υ r,t,φ (g) + s ∥ δ Υ r,t,φ , (S24)  \nwhere Υ r,t,φ (g) are the known expressions for geodesic Mino frequencies [51, 83], and we expand at fixed p, e . We first define the various arguments appearing in the elliptic integrals  \nk 2 0 ≡ 4 Me p +2 M ( e -3) , (S25)  \nα 0 ≡ 2 e ( p -4 M ) (1 + e ) [ p -2 M (3 -e )] , (S26)  \nβ 0 ≡ 16 M 2 e [ p -2 M (3 -e )] [ p -2 M (1 + e )] , (S27)  \nγ 0 ≡ 2 Me [ p -4 M +2 i √ p -3 M (1 + e 2 ) ] [ p -2 M (3 -e )] [ M (1 + e ) + i √ p -3 M (1 + e 2 ) ] . (S28)  \nThe frequency shifts are then  \nδ Υ r = πM √ ( p +2 M ( e -3))(( p -2 M ) 2 -4 M 2 e 2 ) [ ( p -2 M ( e +3)) K ( k 0 ) -( p -M ( e 2 +3)) E ( k 0 ) ] 4( p -2 M ( e +3))( p -M ( e 2 +3)) 3 / 2 K ( k 0 ) 2 , (S29)  \nδ Υ t = p √ M 2 √ p -M ( e 2 +3) ( p Π ( 4 e ( p -2 M (1 -e ) , k 0 ) ( p -2 M (1 -e )) K ( k 0 ) -1 ) -√ M A ( p, e ) E ( k 0 ) 2 + B ( p, e ) E ( k 0 ) + C ( p, e ) 8 M ( p -4 M ) 3 ( p -M ( e 2 +3)) 3 / 2 K ( k 0 ) 2 , (S30)  \nδ Υ φ = (3 + e 2 ) M √ p [( p -2 M ) 2 -4 M 2 e 2 ] 2 p ( p -M ( e 2 +3)) 3 / 2 , (S31)  \nwhere  \n( p, e ) = 2 p ( p -M ( e 2 +3))( p -2 M ( e +1))( p -4 M ) 2 ( p +2 M ( e -3))( p -2 M (1 -e )) (1 -e 2 )( p -2 M ( e +3)) ,  \nA (S32) B ( p, e ) = -2( p -4 M ) (1 -e )(1 + e ) 2 [ (1 + e ) p ( p +2 M ( e -3)) K ( k 0 ) ( e 4 M 2 ( p +4 M ) -2 e 2 M ( p -3 M )( p +4 M ) -28 M 3 + +33 M 2 p +2 p 2 ( p -7 M ) ) +2( p -M ( e 2 +3))( p -2 M (1 -e )) ( (1 -e )(1 + e ) 2 M ( p -4 M ) 2 Π( γ 0 , k 0 )+ +( p -2 M (1 + e ))(8 M 2 + Mp -p 2 + e 2 M (3 p -8 M ))Π( α 0 , k 0 ) )] ,  \n(S33)  \nC ( p, e ) = K ( k 0 ) (1 -e )(1 + e ) 2 ( p -2 M (1 + e )) [ ( p -2 M ( e +3)) ( (1 -e )(1 + e ) 2 ( p -4 M ) 4 ( -4 M 2 -( p -3 M ) p + e 2 M ( p +4 M ))Π( γ 0 , k 0 )+ +4( p -2 M (1 + e ))Π( α 0 , k 0 ) ( 3 M 3 e 6 p 2 + e 4 M 2 ( p ( -6 p 2 +3 Mp +80 M 2 ) -128 M 3 ) -128 M 5 + +144 M 4 p -p 4 ( p -6 M ) -2 M 2 p 3 -55 M 3 p 2 + Me 2 (6 p 4 -40 Mp 3 +113 M 2 p 2 -224 M 3 p +256 M 4 ) )) + +(1 + e )( p -2 M (1 + e )) pK ( k 0 ) ( 24 e 6 M 4 p + e 4 ( M ( p -10 M ) p 3 +184 M 4 p -32 M 3 p 2 -128 M 5 ) --e 2 (952 M 4 p -288 M 3 p 2 +76 M 2 p 3 +( p -16 M ) p 4 -1280 M 5 ) + 3( p -19 M ) p 4 +374 M 3 p 3 --1152 M 3 p 2 +1768 M 4 p -1152 M 5 )] . (S34)  \nThe radial shift δ Υ r agrees with the numerical integration of eq. (B26) in Ref. [43]. We checked the azimuthal shift δ Υ φ and the shift δ Υ t against the frequency shifts obtained from the Hamilton-Jacobi formalism [59]. The precession frequency written explicitly to leading order in terms of p, e reads  \nΥ ψ = √ Mp √ ( p -2 M ) 2 -4 M 2 e 2 K ( k 0 ) [ 4 MK ( k 0 ) + ( p -2 M ( e +3)) ( √ M (1 -e ) √ p -M ( e 2 +3) Im [Π( γ 0 , k 0 )] + Re [Π( γ 0 , k 0 )] ) ] , (S35)  \nwhich agrees with eq. (64) of Ref. [51].  \nUsing the previous expressions, it is easy to compute the nodal and periastron precession ν nodal and ν peri , respectively, which are defined as  \nν nodal = 2 π ( Υ φ Υ ψ -1 ) , ν peri = 2 π ( Υ φ Υ r -1 ) , (S36)  \nThe average polar fundamental frequency is Υ ψ , as can be seen from eq. (26). This leads to the explicit formulas  \nν nodal = 2 π ( 1 K ( p, e ) √ p ( p -2 M (1 -e ))( p -2 M (1 + e )) p -M ( e 2 +3) K ( k 0 ) -1 ) , (S37)  \nwhere  \nK ( p, e ) = 4 MK ( k 0 ) + ( p -2 M ( e +3)) ( √ M (1 -e ) √ p -M ( e 2 +3) Im [Π( γ 0 , k 0 )] + Re [Π( γ 0 , k 0 )] ) . (S38)  \nWe stress that ν nodal is evaluated at leading s ∥ = s = 0 order, since the precession itself is an O ( s ) effect. On the other hand, as can be seen from eq. (46), φ will evolve by exactly the same amount between two pericenter passages and no averaging needs to be invoked. The non-averaged periastron precession is then given by  \nν peri = ν g peri + s ∥ δν peri , (S39)  \nwhere  \nν g peri = 4 pK ( k 0 ) √ p ( p +2 M ( e -3)) -2 π , (S40)  \nδν peri = 2 √ ( p -2 M ) 2 -4 M 2 e 2 p ( p +2 M ( e -3)) ( E ( k 0 ) p -2 M ( e +3) -K ( k 0 ) p ) . (S41)', 'Linear shift to the orbital motion': "The expansion of the orbital motion in spin is non-unique. We choose to expand at fixed radial turning points (parametrized by p, e ) and fixed angle variables (homogeneous orbital phases) q t , q r , q φ , q ψ . This is because the angle variables smoothly parametrize the orbital motion and avoid singularities at turning points associated, e.g., with expansions at fixed velocities and coordinates u µ , x µ . Furthermore, this separates the 'oscillating' part of the spin corrections from the secular/averaged effect of the spin expressed by the frequency shifts δ Υ t,r,φ . By this method we obtain x = x (g) + s ∥ δx for x = t, r, φ . On the other hand, ψ is consistent to compute only to leading order, since it appears only in O ( s ) observables. The resulting expressions read  \nr (0) = p [ p -2 M (3 -e ) -4 Me sin 2 χ 0 ] p -2 M (3 -e 2 ) + e ( p -4 M ) cos 2 χ 0 ,  \nχ 0 ≡ am ( q r π K ( k 0 ) , k 0 ) .  \n(S42) t (0) = q t + √ ( p -2 M ) 2 -4 e 2 M [ p -2 M (3 -e )] [ p [ p 2 -14 Mp +2 Me 2 ( p -2 M ) + 36 M 2 ] ( F ( χ 0 , k 0 ) -q r π K ( k 0 ) ) ( e 2 -1)( p -4 M ) 2 + 2 [ p -2 M (3 -e )] [ p 2 -Mp (1 + 3 e 2 ) + 8 M 2 ( e 2 -1) ] ( Π( χ 0 , α 0 , k 0 ) -q r π Π( α 0 , k 0 ) ) (1 -e 2 )(1 + e )( p -4 M ) 2 -2 M [ p -2 M (3 -e )] ( Π( χ 0 , β 0 , k 0 ) -q r π Π( β 0 , k 0 ) ) p -2 M (1 -e ) + Mp [ p -2 M (3 -e )] ( E ( χ 0 , k 0 ) -q r π E ( k 0 ) ) (1 -e 2 )( p -4 M ) 2 -ep sin 2 χ 0 √ [ p -2 M (3 -e )] [ p -2 M (3 -e cos 2 χ 0 )] (1 -e 2 ) [ p -2 M (3 -e 2 ) + e ( p -4 M ) cos 2 χ 0 ] ] , (S43) φ (0) = q φ +2 √ p p -2 M (3 -e ) ( F ( χ 0 , k 0 ) -q r π K ( k 0 ) ) , (S44) ψ = q ψ + 2 √ [( p -2 M ) 2 -4 M 2 e 2 ] [ p -2 M (3 -e )] [ 4 M √ p -M (3 + e 2 ) ( F ( χ 0 , k 0 ) -q r π K ( k 0 ) ) + √ M (1 -e ) [ p -2 M (3 + e )] Im ( Π( χ 0 , γ 0 , k 0 ) -q r π Π( γ 0 , k 0 ) ) + √ p -M (3 + e 2 ) [ p -2 M (3 + e )] Re ( Π( χ 0 , γ 0 , k 0 ) -q r π Π( γ 0 , k 0 ) ) ] , (S45) δr = -e √ p [( p -2 M ) 2 -4 e 2 ] [ p -2 M (3 -e )] sin 2 χ 0 √ 1 -k 2 0 sin 2 χ 0 ( E ( χ 0 , k 0 ) -q r π E ( k 0 ) ) [ p -2 M (3 -e 2 ) + e ( p -4 M ) cos 2 χ 0 ] 2 , (S46) (S47)  \nThe δt, δφ functions are too long to be typed here explicitly, but they can be found in machine-readable form in the accompanying file SchwarzschildSpinAnalytic.m .  \nIt should be noted that if one is only interested in the approximate orbit, it is computationally more efficient to use the nonexpanded formulas (42)-(54) presented in the main text. On the other hand, in a two-timescale expansion of a gravitational-wave inspiral it is convenient to separate the various pieces contributing to the motion separately, which is the main reason to provide these formulas.  \nFinally, the expanded results can now be substituted to obtain δϑ ( q r , q ψ , p, e ) as  \nδϑ = √ M ( s 2 -s 2 ∥ ) [ Mp 2 +( p -M (3 + e 2 )) r 2 (0) ] sin ψ Mpr (0) . (S48)"}