© The Authors 2026

1. Introduction

After years of monitoring S stars in the central parsec of the Milky Way, observational studies have demonstrated the presence of a supermassive compact object called Sgr A* at the center of the Galaxy (Eckart & Genzel 1996; Ghez et al. 1998, 2003, 2008; Schödel et al. 2002; Gillessen et al. 2009, 2017; GRAVITY Collaboration 2022), which is very likely a supermassive black hole (SMBH). According to the “no-hair” theorem, a black hole can be completely characterized by only three externally observable classical parameters: mass, angular momentum (hereafter, referred to as the spin), and electric charge. We assume that the BH charge can be neglected and, thus, we consider that the BH could be fully characterized by only its mass and spin and, furthermore, it can described within the framework of general relativity (GR) by the Kerr metric. This theorem implies that if we consider a higher moment of the gravity field, such as the mass quadrupole, it must be linked to the mass and spin, meaning that this theorem can be tested by independently measuring these three quantities and verifying the relation between them. Dozens of S star orbits are currently known (Gillessen et al. 2017), including the highly elliptical one of the star S2 with a 16-year period, reaching R ≈ 1400R_S ≈ 120 AU from Sgr A* at its pericenter, where R_S = 2GM_•/c², with G and M_• being the gravitational constant and the black hole’s mass respectively. The mass of Sgr A* had already been measured from AO-assisted imaging and spectroscopy of S stars, but by combining the infrared light collected by the four Unit Telescopes of the Very Large Telescope (VLT) at Paranal, the interferometric instrument GRAVITY has provided a much more precise determination, M_• = (4.2996 ± 0.0118)×10⁶ M_⊙ (GRAVITY Collaboration 2024). This increase in precision is a crucial step: it brings us to the level where future monitoring of S stars could, in principle, constrain the spin and quadrupole moment of the black hole (Will 2008; Angélil & Saha 2014), thereby enabling tests of the “no-hair” theorem.

As stated in Yu et al. (2016), understanding the spin direction of a supermassive black hole is also important because it can reveal clues about its growth history. For example, if the black hole spin is aligned with the young stellar disk in the Galactic Center, it could suggest that the black hole grew mainly from a gas disk that once matched the stellar disk. However, if the spin direction is very different from the stellar disk, it could mean that the black hole’s growth came from multiple chaotic accretion events, rather than one major episode. The different relativistic effects that can be observed in the vicinity of a strong gravitational field, such as the one surrounding Sgr A*, include:

• Schwarzschild precession;

• Shapiro time delay;

• Relativistic redshifts (transverse Doppler shift in special relativity as well as the gravitational redshift appearing in GR);

• Gravitational lensing effect;

• Relativistic aberration.

In the case of a rotating black hole, we need to add the Kerr effects that arise from the rotation of the black hole to this list, namely:

• Lense-Thirring (LT) effect, which results from the rotation of the body (e.g., a spinning black hole) distorts spacetime, causing nearby inertial frames to be “dragged” along with the rotation (also known as the frame-dragging effect). It impacts both the star orbit and photon trajectory and varies depending on the norm and the direction of the spin relative to the star’s orbit and affects both the astrometry and the spectroscopy;

• Quadrupole moment effect, which is due to the oblateness of the black hole arising from its rotation. Similarly to the LT effect, this one also impacts both the star orbit and photon trajectory, varying depending on the norm and the direction of the spin relative to the star’s orbit, affecting both the astrometry and the spectroscopy.

In Sect. 3.1, we discuss the post-Newtonian (PN) expression of the Schwarzschild and spin-related accelerations. From a precession point of view, the Schwarzschild precession is the most dominant relativistic effect, leading to a pericenter advance. If the black hole is rotating, then the LT and the quadrupole moment effects end up generating both apsidal (in-plane) precession and nodal (out-of-plane) precession; as opposed to the Schwarzschild case, which only has an in-plane effect. This absence of out-of-plane precession in the Schwarzschild case is a direct consequence of the spherical symmetry of the Schwarzschild solution. The detection of the Schwarzschild precession was made possible with S2 and other stars (GRAVITY Collaboration 2020, 2024); however, the higher order Kerr effects prove to be more challenging, especially with respect to the quadrupole moment. By measuring the evolution of the orbital orientations of multiple closer-in stars (Will 2008), it is possible to constrain the value and orientation of the spin, along with the magnitude of the quadrupole moment. However, considering how quickly these higher order effects fall with distance from Sgr A*, we would likely need to work with stars that have shorter periods and/or higher eccentricities if we want to test the no-hair theorem in a relatively short timescale (Grould et al. 2017; Waisberg et al. 2018). Therefore, we considered putative stars multiple times closer to Sgr A* and, thus, much more affected by Kerr effects. Such stars might exist if they are too faint to have been already detected by GRAVITY. Furthermore, GRAVITY+ (the upgrade of GRAVITY and of the VLTI infrastructure) is now close to completion and will very soon have the potential to track stars on orbits closer to the black hole than S2, which have thus far been too faint to observe (GRAVITY Collaboration 2022). Therefore, it might become possible to not only constrain the spin parameters, but also the quadrupole moment. We note that there are different effects that have the most potential to interfere with the detection the LT or quadrupole moment effects, but that we did not consider in this work, include (Alush & Stone 2022): mass-precession (Merritt et al. 2010), star spin (Dixon 1970), vector resonant relaxation (Kocsis & Tremaine 2015; Merritt et al. 2010), and tidal disruption (Psaltis et al. 2013; Fabrycky & Tremaine 2007).

Astrometric and spectroscopic data both have the power to detect spin and quadrupole relativistic effects. If current spectroscopic data do not yet have the required precision to detect these effects at the present time, the perspective of the exquisite spectroscopic precision of ELT/MICADO ESO (2025) will most likely allow a very significant contribution of spectroscopy in the constraint of Sgr A*’s spin. However, in this article, we restricted ourselves to studying astrometric effects only. This choice is motivated by two reasons. First, GRAVITY+ already has the power to strongly constrain the spin of Sgr A* based on putative future-detected closer-in stars, without necessarily needing to wait for the first light of MICADO. Second, even in the ELT era, astrometry will still play a crucial role in stellar fittings aimed at constraining Sgr A*’s spin. Indeed, the impact of spectroscopy for spin detection is particularly striking when data are available close to the pericenter passage. In this situation, spectroscopy might overcome the power of astrometry. On the contrary, astrometry is able to significantly constrain the spin effect along the full trajectory. As a consequence, astrometry and spectroscopy appear as two very complementary probes, with astrometry being particularly crucial if spectroscopic data closest to pericenter cannot be secured.

This paper is organized as follows. We present in Sect. 2 an overview of the framework used in our study. The derivation of the relativistic orbital precessions in the PN approximation is laid out in Sect. 3. In Sect. 4, we use the OOGRE code to study how the spin and quadrupole moment affect the evolution of orbits. Finally, we summarize and discussing the applicability of our results to upcoming observations in Sect. 5.

2. Reference frames

We begin this work by properly laying out the framework and present the four reference frames¹ that are needed. Following the convention of Poisson & Will (2014) and Merritt (2013), we introduced a so-called fundamental frame (X, Y, Z) and used it as an arbitrary reference for describing the problem. Considering that we eventually want to model the observable in the sky plane, it is natural to choose the observer’s frame as the fundamental one. For this purpose (as illustrated in Fig. 1), we took Z as pointing towards the observer and (X, Y) = (δ, α) = (Dec, RA) corresponding to the observer’s screen with α or RA and δ or Dec referring to the right ascension and declination, respectively. In other words, we defined the fundamental plane as the sky plane and center it on the apparent position of Sgr A*. This choice has a crucial impact on the definition of orbital parameters considering that the angular parameters are defined with respect to the fundamental frame. In the literature, we can find different definitions of the observer’s frame and fundamental frame. This leads to different results that cause apparent contradictions, as noted in the following sections and appendices.

The problem we study in this work refers to a star in a bound orbit around a rotating black hole. Thus, we need to describe two vectors in the fundamental frame: the angular momentum of the orbit and the angular momentum of the black hole, which naturally define two other frames, namely, that of the orbit and that of the black hole.

The orbit frame (x_orb, y_orb, z_orb) is illustrated by Fig. 1. The axis z_orb is along the angular momentum of the orbit, such that the plane (x_orb, y_orb) matches the instantaneous plane² of the orbit. Also, x_orb always points to the instantaneous pericenter of the orbit (we recall that the pericenter of a relativistic orbit is not fixed; see the notion of osculating orbital parameters noted in Appendix A). Further details on the line of nodes ℒ_orb//sky and the ascending node of the orbit 𝒜𝒩_orb//sky are provided in Appendix B. The orbital parameters that are included in Fig. 1 are discussed further in Appendix A.1. For the time being, it is useful to note that the orientation of an orbit is characterized by three angles: the two angles ι and Ω define the orientation of the orbital plane relative to the fundamental frame, while the third angle ω defines the position of the pericenter inside the orbital plane. In addition, it is useful to introduce the Gaussian frame³ (n_orb, m_orb, z_orb), which follows the movement of the star as illustrated in Fig. 1; n_orb is the radial vector and m_orb completes the orthonormal basis.

The black hole frame (x_bh, y_bh, z_bh), labeled in cartesian coordinates, is illustrated in Fig. 2. It is defined such that the angular momentum (or spin vector) of the black hole is along z_bh. The orientation of x_bh and y_bh is arbitrary. We introduce the angles (θ, β) that allow to orient z_bh in the orbit frame as defined in Fig. 2:

• θ ∈ [0; π] is the rotation angle from z_orb to z_bh;

• β ∈ [0; 2π] is the rotation angle about z_orb, from ℒ_orb//sky to z_bh//orb.

We also introduce an alternative definition of angles (i′,Ω′) in Appendix B that are useful for performing multi-star fits. These angles are not considered further in this article.

3. Post-Newtonian expressions of the secular shifts

3.1. Post-Newtonian acceleration

Here, we focus on the PN formalism, which allows us to separate the LT effect from the quadrupole moment effect. For the PN formalism, we used harmonic coordinates, which are commonly used in works on PN (Poisson & Will 2014). In the PN approximation, the equation of motion of a star of negligible mass, m_star, compared to the mass of the black hole, M_•, (Poisson & Will 2014) is expressed as

$$\begin{array}{c}\hfill \ddot{\mathbf{r}}=-\frac{\mathit{Gm}}{r^2}{\mathbf{n}}_{\mathrm{orb}}+{\mathbf{a}}_{\mathrm{PN}},\end{array}$$(1)

where m = M_• + m_star ≈ M_•. This approximation is referred to as ❶ and can be justified by the fact that the small PN parameter ϵ = GM/(c²p) is at least one or two orders of magnitude greater than η = m_star/M_• ∼ 10⁻⁶ − 10⁻⁵ for any hypothetical star relevant for our study (m_star ∼ 1 − 10 M_⊙ and with p ≤ p^S2). The separation vector from the black hole to the star is denoted by r, while n_orb = r/r is expressed as

$$\begin{array}{c}\hfill r=|\mathbf{r}|=\frac{p}{1+ecosf},\end{array}$$(2)

where p = a_sma(1 − e²) is the semi-latus rectum. Sometimes, p and e are alternatively defined as constants even in GR using the pericenter and apocenter radii r₋ and r₊ by: p = r₋(1 + e) = r₊(1 − e). However, this definition will yield a different definition of the pericenter advance compared to p = a_sma(1 − e²), (see Tucker & Will 2019). A dot in (1) denotes derivation with respect to the harmonic-coordinate time t and the PN acceleration a_PN accounts for deviations from Keplerian motion due to GR. Thanks to the use of the osculating orbital elements method, we are able to use the Newtonian relation of Eq. (2) at each date, while keeping in mind that the involved parameters are time-varying as well as the fact that we are using the weak-field and small-velocity approximations when working with the PN formalism.

If we want to understand, at leading order, how the Schwarzschild precession and the spin-related effects impact the orbit of S stars, it is sufficient to look at the dominant order of the analytical expressions of the secular shifts of orbital elements (Will 2008). In practice, however, since the PN formalism appears to be the most suitable for the process of fitting the data of S stars⁴, we have to choose one specific PN order for all effects to be consistent in the fitting process. As discussed in this section, the dominant order of the LT and quadrupole moment effects are at the 1.5PN (corresponding to (v/c)³) and 2PN orders (corresponding to (v/c)⁴), respectively. Therefore, we need to first obtain the analytical expressions up to the 2PN order.

Considering the second-order PN approximation (2PN) corresponding to (v/c)⁴, the PN acceleration can be decomposed into three parts:

$$\begin{array}{c}\hfill {\mathbf{a}}_{\mathrm{PN}}\approx {\mathbf{a}}_{2\mathrm{PN}}={\mathbf{a}}_{\mathrm{Sch}}^{2\mathrm{PN}}+{\mathbf{a}}_{\mathrm{LT}}^{1.5\mathrm{PN}}+{\mathbf{a}}_Q^{2\mathrm{PN}},\end{array}$$(3)

with ${\mathbf{a}}_{\mathrm{Sch}}^{2\mathrm{PN}}$, ${\mathbf{a}}_{\mathrm{LT}}^{1.5\mathrm{PN}}$, ${\mathbf{a}}_Q^{2\mathrm{PN}}$ corresponding to the 2PN acceleration of the star due to the Schwarzschild precession, the spin and the quadrupole moment, respectively. This approximation is referred to as ❷. The kPN notation indicates the PN order up to which the various effects enter the equation of motion. For instance ${\mathbf{a}}_{\mathrm{Sch}}^{2\mathrm{PN}}$ contains both the 1PN and 2PN Schwarzschild contributions. By considering approximation ❶ and by defining the small PN parameter as

$$\begin{array}{c}\hfill ϵ=\frac{\mathit{Gm}}{c^{2p}},\end{array}$$(4)

the accelerations can be expressed in terms of ϵ as (Merritt 2013; Will 2008; Blanchet & Iyer 2003):

$$\begin{array}{cc}& {\mathbf{a}}_{\mathrm{Sch}}^{2\mathrm{PN}}=\frac{ϵ}{p}{(1+ecosf)}^2((4\frac{\mathit{Gm}}{r}-v^2){\mathbf{n}}_{\mathrm{orb}}+4v_r\mathit{v})\hfill \\ \hfill & +\frac{{ϵ}^2}{p}{(1+ecosf)}^3((2v_r^2-9\frac{\mathit{Gm}}{r}){\mathbf{n}}_{\mathrm{orb}}-2v_r\mathit{v}),\hfill \end{array}$$(5)

$$\begin{array}{cc}& {\mathbf{a}}_{\mathrm{LT}}^{1.5\mathrm{PN}}=-2{ϵ}^{3/2}\sqrt{\frac{\mathit{Gm}}{p^3}}{(1+ecosf)}^3\chi \hfill \\ \hfill & \times [2\mathit{v}\times {\mathbf{z}}_{\mathrm{bh}}-3({\mathbf{n}}_{\mathrm{orb}}·\mathit{v}){\mathbf{n}}_{\mathrm{orb}}\times {\mathbf{z}}_{\mathrm{bh}}-3{\mathbf{n}}_{\mathrm{orb}}({\mathbf{n}}_{\mathrm{orb}}\times \mathit{v})·{\mathbf{z}}_{\mathrm{bh}}],\hfill \end{array}$$(6)

$$\begin{array}{cc}& {\mathbf{a}}_Q^{2\mathrm{PN}}=\frac{3}{2}{ϵ}^2\frac{\mathit{Gm}}{p^2}{(1+ecosf)}^4{\chi }^2\hfill \\ \hfill & \times (5{\mathbf{n}}_{\mathrm{orb}}{({\mathbf{n}}_{\mathrm{orb}}·{\mathbf{z}}_{\mathrm{bh}})}^2-2({\mathbf{n}}_{\mathrm{orb}}·{\mathbf{z}}_{\mathrm{bh}}){\mathbf{z}}_{\mathrm{bh}}-{\mathbf{n}}_{\mathrm{orb}}),\hfill \end{array}$$(7)

with v² = v_r² + v_t² and χ being the spin parameter. The latter is positive and has a maximal value of 1 (for a Kerr black hole). The dimensionless spin vector, χ = χz_bh, is linked to the angular momentum, J = Jz_bh, of a black hole of mass, M_•, via

$$\begin{array}{c}\hfill \mathit{\chi }=\mathit{J}\frac{c}{G{M}_{\bullet }^2}·\end{array}$$(8)

As for the quadrupole moment parameter Q, we note that it enters Eq. (7) via

$$\begin{array}{c}\hfill Q=-\frac{1}{c^2}\frac{J^2}{M_{\bullet }}=-\frac{G^2{M}_{\bullet }^3}{c^4}{\chi }^2,\end{array}$$(9)

where the fact that Q depends only on M_• and χ is a consequence of the no-hair theorem (Poisson 1998; Krishnendu et al. 2019).

We note that the Schwarzschild acceleration is composed of two terms, namely, orders of 1PN and 2PN. The LT and quadrupole moment accelerations are composed of only one term of an order of 1.5PN and 2PN, respectively. The leading order of the various effects is thus 1PN, 1.5PN, 2PN, for the Schwarzschild, LT, and quadrupole moment effects, respectively. We need to consider the subleading order only for the Schwarzschild term. The LT and quadrupole moment subleading terms appear at order above 2PN. In this paper, we consider only the Kerr spacetime, and, thus, we plugged the expression of Q given by Eq. (9) into our equations. However, if we want to test the no-hair theorem, we should not assume Eq. (9) when generating the expression of the secular shifts⁵, since Q would need to be measured independently of χ.

3.2. Lagrange planetary equations

We can use the orbital perturbation theory (Merritt 2013; Will 2008) to express the derivative of each time-varying, osculating orbital parameter with respect to time, as a function of the various ℛ, 𝒮, and 𝒲 terms, the perturbative accelerations expressed in the Gaussian frame, introduced in Appendix C (also see Poisson & Will 2014, for details). They are called the Lagrange planetary equations⁶ and are expressed as

$$\begin{array}{cc}& \frac{\mathrm{d}p}{\mathrm{d}t}=2\sqrt{\frac{p^3}{\mathit{Gm}}}\frac{1}{1+ecosf}\mathcal{S};\hfill \end{array}$$(10)

$$\begin{array}{cc}& \frac{\mathrm{d}e}{\mathrm{d}t}=\sqrt{\frac{p}{\mathit{Gm}}}[sinf\mathcal{R}+\frac{2cosf+e(1+{cos}^{2}f)}{1+ecosf}\mathcal{S}];\hfill \end{array}$$(11)

$$\begin{array}{cc}& \frac{\mathrm{d}\iota }{\mathrm{d}t}=\sqrt{\frac{p}{\mathit{Gm}}}\frac{cos(\omega +f)}{1+ecosf}\mathcal{W};\hfill \end{array}$$(12)

$$\begin{array}{cc}& sin\iota \frac{\mathrm{d}\mathrm{\Omega }}{\mathrm{d}t}=\sqrt{\frac{p}{\mathit{Gm}}}\frac{sin(\omega +f)}{1+ecosf}\mathcal{W};\hfill \end{array}$$(13)

$$\begin{array}{cc}& \frac{\mathrm{d}\omega }{\mathrm{d}t}=\frac{1}{e}\sqrt{\frac{p}{\mathit{Gm}}}[-cosf\mathcal{R}+\frac{2+ecosf}{1+ecosf}sinf\mathcal{S}\hfill \\ \hfill & -ecot\iota \frac{sin(\omega +f)}{1+ecosf}\mathcal{W}].\hfill \end{array}$$(14)

No approximations are done to obtain Eqs. (10)–(14) from Eqs. (1).

As stated earlier, in the context of a small perturbing force, it is possible to get a good estimation of the orbital dynamics by only looking at the dominant contribution of each effect. To do so, we solve the equations within the framework of perturbation theory and find that, to express the leading order of such an effect for a given orbital parameter, it is appropriate to take the 0PN terms of all other orbital parameters in the right-hand side of Eqs. (10)–(15), before integrating these equations over t. Even more convenient would be to use f as an independent variable instead of t, considering it has the same fixed integration interval over one period for any star. It evolves as (see Poisson & Will 2014, for details):

$$\begin{array}{cc}\hfill \frac{\mathrm{d}f}{\mathrm{d}t}& ={\left(\frac{\mathrm{d}f}{\mathrm{d}t}\right)}_{\mathrm{Kepler}}-(\frac{\mathrm{d}\omega }{\mathrm{d}t}+cos\iota \frac{\mathrm{d}\mathrm{\Omega }}{\mathrm{d}t})\hfill \\ \hfill & =\sqrt{\frac{\mathit{Gm}}{p^3}}{(1+ecosf)}^2\hfill \\ \hfill & +\frac{1}{e}\sqrt{\frac{p}{\mathit{Gm}}}[cosf\mathcal{R}-\frac{2+ecosf}{1+ecosf}sinf\mathcal{S}],\hfill \end{array}$$(15)

where we recall the fact that f is the angle between the varying pericenter and the position vector of the star r (see Fig. 1). We can make use of the fact that the term ${\left(\frac{\mathrm{d}f}{\mathrm{d}t}\right)}_{\mathrm{Kepler}}=\sqrt{\frac{\mathit{Gm}}{p^3}}{(1+ecosf)}^2={\left(\frac{\mathrm{d}\phi }{\mathrm{d}t}\right)}_{\mathrm{Kepler}}$ is Keplerian (0PN) and that the non-Keplerian terms on the right-hand side of Eq. (15) start at 1PN, 1.5PN, and 2PN orders for the Schwarzschild, LT, and quadrupole moment effects, respectively. This means that we are allowed to express the temporal variation as

$$\begin{array}{c}\hfill \frac{\mathrm{d}t}{\mathrm{d}f}=\sqrt{\frac{p^3}{\mathit{Gm}}}\frac{1}{{(1+ecosf)}^2}[1+\mathcal{O}(ϵ)],\end{array}$$(16)

before multiplying Eqs. (10)–(14) by the zeroth-order term of Eq. (16) to get the dominant order derivative of each orbital parameter with respect to the true anomaly, f.

3.3. Analytical integration of the planetary equations

Now that we have the expression of the planetary equation as a function of f, we integrated them with f from f₀ to f₀ + 2π⁷ to find the dominant terms of the secular evolution of the orbital parameters. By considering the Keplerian dominant constant term of p, e, ω, and β on the right-hand side of Langrange planetary equations, we find that p and e (and, consequently, a_sma and P) do not show any secular variations after one orbit at the dominant orders. As for the other parameters, we find that the precessions per orbit of a star’s argument of periapsis, Δω, along with the longitude of the ascending node, ΔΩ, and inclination, Δι verify⁸, the following:

$$\begin{array}{cc}& \mathrm{\Delta }{\omega }_{\mathrm{Sch}}{|}^{1\mathrm{PN}}=6\pi ϵ,\hfill \end{array}$$(17)

$$\begin{array}{cc}& \mathrm{\Delta }{\omega }_{\mathrm{LT}}{|}^{1.5\mathrm{PN}}=-4\pi {ϵ}^{3/2}(2cos\theta +cot\iota sin\theta sin\beta )\chi ,\hfill \end{array}$$(18)

$$\begin{array}{cc}& \mathrm{\Delta }{\omega }_Q{|}^{2\mathrm{PN}}=-\frac{3\pi }{2}{ϵ}^2(1-3{cos}^2\theta -2cot\iota cos\theta sin\theta sin\beta ){\chi }^2,\hfill \end{array}$$(19)

$$\begin{array}{cc}& \mathrm{\Delta }{\iota }_{\mathrm{Sch}}{|}^{1\mathrm{PN}}=0,\hfill \end{array}$$(20)

$$\begin{array}{cc}& \mathrm{\Delta }{\iota }_{\mathrm{LT}}{|}^{1.5\mathrm{PN}}=4\pi {ϵ}^{3/2}sin\theta cos\beta \chi ,\hfill \end{array}$$(21)

$$\begin{array}{cc}& \mathrm{\Delta }{\iota }_Q{|}^{2\mathrm{PN}}=-3\pi {ϵ}^{2}cos\theta sin\theta cos\beta {\chi }^2,\hfill \end{array}$$(22)

$$\begin{array}{cc}& sin\iota \mathrm{\Delta }{\mathrm{\Omega }}_{\mathrm{Sch}}{|}^{1\mathrm{PN}}=0,\hfill \end{array}$$(23)

$$\begin{array}{cc}& sin\iota \mathrm{\Delta }{\mathrm{\Omega }}_{\mathrm{LT}}{|}^{1.5\mathrm{PN}}=4\pi {ϵ}^{3/2}sin\theta sin\beta \chi ,\hfill \end{array}$$(24)

$$\begin{array}{cc}& sin\iota \mathrm{\Delta }{\mathrm{\Omega }}_Q{|}^{2\mathrm{PN}}=-3\pi {ϵ}^{2}cos\theta sin\theta sin\beta {\chi }^2,\hfill \end{array}$$(25)

We note that Will (2008) adopted the same definition of fundamental frame as we have in our work, which is why we agree on the variation of the angular orbital parameters. However, our Eqs. (21) and (22) should not be compared with Eqs. (3.23a) of Will & Maitra (2017) and (5) of Tucker & Will (2019) since we have made a different choice of fundamental frame from theirs, leading to different definitions of ι and Ω among other parameters.

These secular trends are valid at 2PN order for the LT and quadrupole moment shifts, for which subleading contributions appear at orders higher than 2PN. However, for the Schwarzschild terms, we need to compute the subleading 2PN contribution, in addition to the leading 1PN contribution provided above.

To do so, the two-timescale method detailed in Appendix D can be used to account for the periodic evolution of the orbital parameters. This method requires introducing a small parameter similar to ϵ, but nonvarying, which we label with the orbit-averaged quantity, $\stackrel{\sim }{ϵ}$. Accordingly, the 2PN pericenter advance takes the following form:

$$\begin{array}{cc}& \mathrm{\Delta }{\omega }_{\mathrm{Sch}}{|}^{2\mathrm{PN}}=6\pi \stackrel{\sim }{ϵ}+\frac{3\pi }{2}{\stackrel{\sim }{ϵ}}^2(2+{\stackrel{\sim }{e}}^2),\hfill \end{array}$$(26)

$$\begin{array}{cc}& \mathrm{\Delta }{\iota }_{\mathrm{Sch}}{|}^{2\mathrm{PN}}=0,\hfill \end{array}$$(27)

$$\begin{array}{cc}& \mathrm{\Delta }{\mathrm{\Omega }}_{\mathrm{Sch}}{|}^{2\mathrm{PN}}=0,\hfill \end{array}$$(28)

which are more precise analytical expressions than those given in Alush & Stone (2022). When expressing the pericenter advance to 2PN accuracy, it is therefore essential to use the orbit-averaged parameters $\stackrel{\sim }{ϵ}$ and $\stackrel{\sim }{e}$, rather than their osculating or varying counterparts. Using the latter in the 1PN term would already introduce discrepancies at 2PN order. At leading order, however, we can safely identify $ϵ\simeq \stackrel{\sim }{ϵ}$, so this subtlety has no impact on Eqs. (17)–(25). A more detailed justification of this choice is provided in Appendix D.2. We note that if ι = 0, ℒ_orb//sky cannot be defined as it is degenerate, then ω and Ω, and thus Eqs. (18), (19), (24), and (25) cannot be defined either. This issue can be addressed by introducing new variables, as explained in Sect. 4.

3.4. Numerical integration of the planetary equations

To numerically integrate the planetary equations to simulate orbits, we used our Python-based code called Osculating Orbits in General Relativistic Environments (OOGRE), which was originally developed by Heißel et al. (2022). It is a perturbed Kepler-orbit model code that is based on the formalism presented in Sect. 3.1 and optimized for stars in the Galactic Center.

The goal of Sect. 4 is to investigate the impact of Kerr effects on the geometry of the orbit; the interest primarily lies in the dynamical interaction between the black hole and one star. Therefore, we do not consider the effect of mass-precession and vector resonant relaxation, while also neglecting the tidal disruption and effect of the star’s spin. Moreover, as the goal is not to simulate observations here, we turned off the following effects, as we find they would interfere with this goal:

• the Shapiro time delay;

• the Rømer effect;

• the gravitational lensing⁹.

Also, the relativistic aberration is not implemented in OOGRE and is not needed for the purposes of this discussion either. However, we extended this code, which initially only considers ${\mathbf{a}}_{\mathrm{Sch}}^{1\mathrm{PN}}$, to the Kerr metric by adding the spin and quadrupole moment contributions, ${\mathbf{a}}_{\mathrm{LT}}^{1.5\mathrm{PN}}$ and ${\mathbf{a}}_Q^{2\mathrm{PN}}$, to simulate these effects. We also pushed the Schwarzschild precession to the 2PN order with ${\mathbf{a}}_{\mathrm{Sch}}^{2\mathrm{PN}}$ for consistency.

We chose a set of initial conditions as in Heißel et al. (2022), and fix the osculating time as t_osc = t_peri − P_osc/2, an approximate¹⁰ value of the apocenter, granting the code more stability. The code integrates Eqs. (10)–(14) at the core of OOGRE while only using the approximations ❶ and ❷.

4. Spin impact on the orientation of the orbit

First, we want to understand analytically the effect of spin on the evolution of the orientation of the orbit. This leads us to introduce several precession movements, either instantaneous or secular. The secular evolution is well known (see e.g., Kocsis & Tremaine 2011) and consists of a precession movement of z_orb about z_bh. Instantaneous evolution is both less known and more useful in the perspective of orbital fits. Therefore, we focus on the latter in the discussion below. Second, we want to illustrate this analytical understanding by performing simulations for a hypothetical closer-in star, which we call “S2/10”.

4.1. Instantaneous orientation change

4.1.1. Orbital frame orientation change

We start by characterizing the evolution of the orientation of the orbit by expressing the temporal evolution of the orbit frame. Using the Lagrange planetary equations, it is possible to reexpress the temporal variations of these vectors as functions of the temporal variations of the orbital elements. We get (see Will et al. 2023):

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{\mathbf{x}}_{\mathrm{orb}}}{\mathrm{d}t}& =\frac{\mathrm{d}\varpi }{\mathrm{d}t}{\mathbf{y}}_{\mathrm{orb}}+\frac{\mathrm{d}\mathrm{\Theta }}{\mathrm{d}t}{\mathbf{z}}_{\mathrm{orb}}\hfill \\ \hfill & =\frac{\mathrm{d}\varpi }{\mathrm{d}t}{\mathbf{z}}_{\mathrm{orb}}\times {\mathbf{x}}_{\mathrm{orb}}-\frac{\mathrm{d}\mathrm{\Theta }}{\mathrm{d}t}{\mathbf{y}}_{\mathrm{orb}}\times {\mathbf{x}}_{\mathrm{orb}},\hfill \end{array}$$(29)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{\mathbf{y}}_{\mathrm{orb}}}{\mathrm{d}t}& =-\frac{\mathrm{d}\varpi }{\mathrm{d}t}{\mathbf{x}}_{\mathrm{orb}}+\frac{\mathrm{d}\mathrm{\Xi }}{\mathrm{d}t}{\mathbf{z}}_{\mathrm{orb}}\hfill \\ \hfill & =-\frac{\mathrm{d}\varpi }{\mathrm{d}t}{\mathbf{y}}_{\mathrm{orb}}\times {\mathbf{z}}_{\mathrm{orb}}+\frac{\mathrm{d}\mathrm{\Xi }}{\mathrm{d}t}{\mathbf{x}}_{\mathrm{orb}}\times {\mathbf{y}}_{\mathrm{orb}},\hfill \end{array}$$(30)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{\mathbf{z}}_{\mathrm{orb}}}{\mathrm{d}t}& =-s_{\mathcal{W}}\sqrt{{\left(\frac{\mathrm{d}\iota }{\mathrm{d}t}\right)}^2+{sin}^2\iota {\left(\frac{\mathrm{d}\mathrm{\Omega }}{\mathrm{d}t}\right)}^2}{\mathbf{m}}_{\mathrm{orb}}\hfill \\ \hfill & =-\frac{\mathrm{d}\mathrm{\Theta }}{\mathrm{d}t}{\mathbf{x}}_{\mathrm{orb}}+\frac{\mathrm{d}\mathrm{\Xi }}{\mathrm{d}t}{\mathbf{y}}_{\mathrm{orb}}\hfill \\ \hfill & =-\frac{\mathrm{d}\mathrm{\Theta }}{\mathrm{d}t}{\mathbf{y}}_{\mathrm{orb}}\times {\mathbf{z}}_{\mathrm{orb}}-\frac{\mathrm{d}\mathrm{\Xi }}{\mathrm{d}t}{\mathbf{x}}_{\mathrm{orb}}\times {\mathbf{z}}_{\mathrm{orb}},\hfill \end{array}$$(31)

where s_𝒲 is the sign of 𝒲 and

$$\begin{array}{cc}& \frac{\mathrm{d}\varpi }{\mathrm{d}t}\equiv \frac{\mathrm{d}\omega }{\mathrm{d}t}+cos\iota \frac{\mathrm{d}\mathrm{\Omega }}{\mathrm{d}t};\hfill \end{array}$$(32)

$$\begin{array}{cc}& \frac{\mathrm{d}\mathrm{\Theta }}{\mathrm{d}t}\equiv sin\omega \frac{\mathrm{d}\iota }{\mathrm{d}t}-sin\iota cos\omega \frac{\mathrm{d}\mathrm{\Omega }}{\mathrm{d}t};\hfill \end{array}$$(33)

$$\begin{array}{cc}& \frac{\mathrm{d}\mathrm{\Xi }}{\mathrm{d}t}\equiv -cos\omega \frac{\mathrm{d}\iota }{\mathrm{d}t}-sin\iota sin\omega \frac{\mathrm{d}\mathrm{\Omega }}{\mathrm{d}t}·\hfill \end{array}$$(34)

We can add the last lines of Eqs. (29)–(31) to reveal the precession aspect of the equations and their different components.

Thus, we obtain

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\mathbf{x}}_{\mathrm{orb}}}{\mathrm{d}t}·{\mathbf{y}}_{\mathrm{orb}}=-\frac{\mathrm{d}{\mathbf{y}}_{\mathrm{orb}}}{\mathrm{d}t}·{\mathbf{x}}_{\mathrm{orb}}=\frac{\mathrm{d}\varpi }{\mathrm{d}t},\end{array}$$(35)

encoding the change of direction of the major and minor axes in the orbit plane;

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\mathbf{x}}_{\mathrm{orb}}}{\mathrm{d}t}·{\mathbf{z}}_{\mathrm{orb}}=-\frac{\mathrm{d}{\mathbf{z}}_{\mathrm{orb}}}{\mathrm{d}t}·{\mathbf{x}}_{\mathrm{orb}}=\frac{\mathrm{d}\mathrm{\Theta }}{\mathrm{d}t},\end{array}$$(36)

encoding the change of the major axis orthogonal to the orbit plane, or equivalently the change of the orbital angular momentum direction along the major axis;

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\mathbf{z}}_{\mathrm{orb}}}{\mathrm{d}t}·{\mathbf{y}}_{\mathrm{orb}}=-\frac{\mathrm{d}{\mathbf{y}}_{\mathrm{orb}}}{\mathrm{d}t}·{\mathbf{z}}_{\mathrm{orb}}=\frac{\mathrm{d}\mathrm{\Xi }}{\mathrm{d}t},\end{array}$$(37)

encoding the change of the minor-axis direction orthogonal to the orbit plane, or equivalently the change of the orbital angular momentum along the minor axis. We note that ${\mathbf{x}}_{\mathrm{orb}}·\frac{\mathrm{d}{\mathbf{x}}_{\mathrm{orb}}}{\mathrm{d}t}={\mathbf{y}}_{\mathrm{orb}}·\frac{\mathrm{d}{\mathbf{y}}_{\mathrm{orb}}}{\mathrm{d}t}={\mathbf{z}}_{\mathrm{orb}}·\frac{\mathrm{d}{\mathbf{z}}_{\mathrm{orb}}}{\mathrm{d}t}=0$, as it should for unit vectors. Therefore, we can decompose the overall precession into two components: the precession of the argument of periapsis within the changing orbital plane, with associated characteristic frequency $\frac{\mathrm{d}\varpi }{\mathrm{d}t}$, hereafter referred to as “in-plane¹¹ precession”; and the precessions of the orbit plane, with characteristic associated frequencies, $\frac{\mathrm{d}\mathrm{\Theta }}{\mathrm{d}t}$ and $\frac{\mathrm{d}\mathrm{\Xi }}{\mathrm{d}t}$, hereafter referred to as “out-of-plane precession”.

These scalar quantities are also presented¹² in Will et al. (2023), and are independent from the observer’s location, which is obvious from their definition: they depend only on (x_orb, y_orb, z_orb), which depends only on the orbit. As such, ϖ, Θ and Ξ appear to be very useful for expressing observer-independent statements on the variation of the pericenter and angular momentum directions.

4.1.2. Secular orientation change

We introduce the secular shifts Δϖ, ΔΘ, and ΔΞ, as the integration over an orbit of the temporal variations of ϖ, Θ, and Ξ. The secular in-plane shift of the pericenter can be expressed at 2PN order as

$$\begin{array}{cc}& \mathrm{\Delta }{\varpi }_{\mathrm{Sch}}{{|}^{2\mathrm{PN}}=\mathrm{\Delta }{\omega }_{\mathrm{Sch}}|}^{2\mathrm{PN}}=6\pi \stackrel{\sim }{ϵ}+\frac{3\pi }{2}{\stackrel{\sim }{ϵ}}^2(2+{\stackrel{\sim }{e}}^2),\hfill \end{array}$$(38)

$$\begin{array}{cc}& \mathrm{\Delta }{\varpi }_{\mathrm{LT}}{|}^{2\mathrm{PN}}=-8\pi {ϵ}^{3/2}cos\theta \chi ,\hfill \end{array}$$(39)

$$\begin{array}{cc}& \mathrm{\Delta }{\varpi }_Q{|}^{2\mathrm{PN}}=-\frac{3\pi }{2}{ϵ}^2(1-3{cos}^2\theta ){\chi }^2.\hfill \end{array}$$(40)

We note that Eqs. (39) and (40) will still be valid if ι = 0, as opposed to Eqs. (18), (19), (24), and (25). Using ϖ allows us to quantify the precession inside the orbital plane without using the notion of ℒ_orb//sky which is degenerate in the case of a face-on¹³ orbit. From Eqs. (39) and (40) we see that the LT and quadrupole moment in-plane precession depend on the polar angle, θ, but not on the azimuthal angle, ψ, of the spin vector spherical coordinates defined in Fig. 2. We see that the in-plane precession of the LT effect is maximal for orbits in the equatorial plane of the black hole, namely, for θ = 0 (mod π), and it vanishes for polar¹⁴ orbits, namely, for θ = π/2 (mod π). However, the quadrupole moment in-plane precession is maximized for both equatorial and polar orbits and it vanishes for $\theta =arccos(\pm 1/\sqrt{3})$.

The secular out-of-plane shift of the major axis can be expressed at 2PN order from Eq. (33) as

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathrm{\Theta }=sin\omega \mathrm{\Delta }\iota -sin\iota cos\omega \mathrm{\Delta }\mathrm{\Omega },\end{array}$$(41)

to characterize the out-of-plane component of the pericenter precession. When using the secular shifts of the orbital parameters derived in Sect. 3.3 in Eq. (41), we obtain the following expressions:

$$\begin{array}{cc}& \mathrm{\Delta }{\mathrm{\Theta }}_{\mathrm{LT}}{|}^{2\mathrm{PN}}=-4\pi {ϵ}^{3/2}sin\theta sin\psi \chi ,\hfill \end{array}$$(42)

$$\begin{array}{cc}& \mathrm{\Delta }{\mathrm{\Theta }}_Q{|}^{2\mathrm{PN}}=3\pi {ϵ}^{2}cos\theta sin\theta sin\psi {\chi }^2,\hfill \end{array}$$(43)

corresponding to the out-of-plane movement of the major axis precession (also referred to as out-of-plane apocenter shift) due to the LT and quadrupole moment effects, as illustrated in Fig. 3. The Schwarzschild precession does not contribute to the out-of-plane shifts at any PN order.

Finally, the secular out-of-plane shift of the minor axis can be expressed at 2PN order from Eq. (34) as

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathrm{\Xi }=-cos\omega \mathrm{\Delta }\iota -sin\omega sin\iota \mathrm{\Delta }\mathrm{\Omega }.\end{array}$$(44)

We obtain the following expressions:

$$\begin{array}{cc}& \mathrm{\Delta }{\mathrm{\Xi }}_{\mathrm{LT}}{|}^{2\mathrm{PN}}=-4\pi {ϵ}^{3/2}sin\theta cos\psi \chi ,\hfill \end{array}$$(45)

$$\begin{array}{cc}& \mathrm{\Delta }{\mathrm{\Xi }}_Q{|}^{2\mathrm{PN}}=3\pi {ϵ}^{2}cos\theta sin\theta cos\psi {\chi }^2.\hfill \end{array}$$(46)

The expressions of ΔΘ and ΔΞ above show that these quantities depend on both θ and ψ. If we take ψ = π/2 (mod π), namely, when the projection of the spin axis of the black hole on the orbit is along the minor axis, then ΔΘ would be maximized, while ΔΞ would be null (i.e., the major axis precesses in a plane orthogonal to the minor axis). Conversely, with ψ = 0 (mod π), namely, when the projection of the spin axis on the orbit is along the major axis, ΔΘ would be null, while ΔΞ would be maximized (i.e., the minor axis precesses in a plane orthogonal to the major axis). We note that due to the in-plane precession, this can be true only instantaneously.

The LT effect generates a maximal out-of-plane precession for polar orbits and a vanishing one for equatorial orbits. We also see that the configuration that allows for the most apocenter displacement due to the LT effect (the maximum of $\sqrt{\mathrm{\Delta }{\varpi }_{\mathrm{LT}}^2+\mathrm{\Delta }{\mathrm{\Theta }}_{\mathrm{LT}}^2}$) would be an equatorial orbit.

As for the quadrupole moment effect, it does not induce any out-of-plane precessions for polar orbits. In addition, we observe that the out-of-plane precession is maximized for θ = π/4 (mod π). This behavior is similar to the Newtonian case of a particle orbiting an oblate mass; if the particle is mid-inclined with respect to the equatorial plane of the central mass¹⁵, the particle would experience the highest level of asymmetries along the orbit leading to the highest out-of-plane shifts among all possible relative orientations (Boain 2004). In terms of overall apocenter angular shift, we also see that the configuration that maximizes $\sqrt{\mathrm{\Delta }{\varpi }_Q^2+\mathrm{\Delta }{\mathrm{\Theta }}_Q^2}$ would be an equatorial orbit.

Even though ω and β are observer-dependent, ψ = β − ω, the angle between z_bh//orb and x_orb, rotating around z_orb in Fig. 3, is observer-independent. Therefore, ΔΘ and ΔΞ are observer-independent as well.

4.2. Secular shifts for S2 and “S2/10”

Here, we consider the star S2, as well as a hypothetical star that can potentially be detected with GRAVITY+ and that we call “S2/10”. It is a star with an orbit identical to that of S2 but with a semimajor axis 10 times smaller. In Table 1, we compare the various maximal secular shifts of the star S2. To compare these shifts using a single set of orbital parameters, we used osculating orbital elements at apocenter.

When computing the secular shifts for “S2/10”, we can see that

$$\begin{array}{cc}\hfill \mathrm{\Delta }{\varpi }_{\mathrm{Sch}}^{\mathrm{S}2/10}& =10\mathrm{\Delta }{\varpi }_{\mathrm{Sch}}{{|}_{1\mathrm{PN}}^{\mathrm{S}2}+{10}^2\mathrm{\Delta }{\varpi }_{\mathrm{Sch}}|}_{2\mathrm{PN}}^{\mathrm{S}2}\approx 10\mathrm{\Delta }{\varpi }_{\mathrm{Sch}}^{\mathrm{S}2},\hfill \\ \hfill \mathrm{\Delta }{X}_{\mathrm{LT}}^{\mathrm{S}2/10}& ={10}^{3/2}\mathrm{\Delta }{X}_{\mathrm{LT}}^{\mathrm{S}2},\hfill \\ \hfill \mathrm{\Delta }{X}_Q^{\mathrm{S}2/10}& ={10}^2\mathrm{\Delta }{X}_Q^{\mathrm{S}2},\hfill \end{array}$$

with X denoting either in-plane (ϖ) or out-of-plane (Θ, Ξ) precession components. This scaling between the Schwarzschild precession and the LT and the quadrupole moment effects arises from their respective radial dependencies, which are encapsulated by the parameter ϵ. While the intrinsic relativistic precessions increase steeply for smaller semimajor axes (as shown by the scaling relations above), the corresponding astrometric signatures on sky do not necessarily become easier to detect. Indeed, in astrometry, we measure projected positions rather than orbital angles, and the apparent size of the orbit scales linearly with a_sma/R₀. For such a star as “S2/10”, the orbit would appear ten times smaller on sky than that of S2, partly compensating for the intrinsic amplification of the relativistic effects. As a result, the gain in the astrometric measurability is only increased by a factor $\sqrt{10}$ and 10 for the LT effect and quadrupole term, while the Schwarzschild precession would remain essentially unchanged in relative detectability. Nevertheless, having smaller periods allows us to observe more pericenter and apocenter passages, thereby accumulating more secular shifts in a given amount of time. For instance, the period of the star S2 (≈16 yr) corresponds to 10^3/2 ≈ 32 orbits of the star “S2/10”. Therefore, when working with “S2/10”, instead of S2 over one period of S2, the astrometric measurability of the Schwarzschild precession is increased by a factor ≈10^3/2 ≈ 32, the LT precessions by a factor ${10}^{3/2}\times \sqrt{10}=100$, and the quadrupole moment shifts by a factor 10^3/2 × 10 ≈ 320. These scaling behaviors emphasize the scientific value of observing closer-in stars, providing a powerful means to probe both the spin and higher order mass moments of Sgr A*.

4.3. Time-averaged orientation change

After describing the instantaneous orientation change of the orbital frame, we go on to discuss its secular evolution over long timescales (as compared to the orbital period). We know from Eqs. (12), (13) and (31), along with the corresponding expression of the perturbative acceleration given in Eq. (C.8), that the LT effect acts on the orbital angular momentum such that

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\mathbf{z}}_{\mathrm{orb}}}{\mathrm{d}t}{|}_{\mathrm{LT}}=-2\frac{G^2{m}^2}{c^3{r}^3}\chi [\frac{esinf}{1+ecosf}({\mathbf{z}}_{\mathrm{bh}}·{\mathbf{m}}_{\mathrm{orb}})+2({\mathbf{z}}_{\mathrm{bh}}·{\mathbf{n}}_{\mathrm{orb}})]{\mathbf{m}}_{\mathrm{orb}}.\end{array}$$(47)

From this we obtain¹⁶:

$$\begin{array}{c}\hfill {\langle \frac{\mathrm{d}{\mathbf{z}}_{\mathrm{orb}}}{\mathrm{d}t}{|}_{\mathrm{LT}}\rangle }_t={\omega }_{\mathrm{LT}}({\mathbf{z}}_{\mathrm{bh}}\times {\mathbf{z}}_{\mathrm{orb}}),\end{array}$$(48)

where¹⁷:

$$\begin{array}{c}\hfill {\omega }_{\mathrm{LT}}=2n{ϵ}^{3/2}\chi ,\end{array}$$(49)

and $n=\sqrt{Gm/a_{\mathrm{sma}}^3}$ is the Newtonian mean motion. Moreover, we have the expression (see Fig. 2):

$$\begin{array}{c}\hfill {\mathbf{z}}_{\mathrm{bh}}=sin\theta cos\psi {\mathbf{x}}_{\mathrm{orb}}+sin\theta sin\psi {\mathbf{y}}_{\mathrm{orb}}+cos\theta {\mathbf{z}}_{\mathrm{orb}},\end{array}$$(50)

which we can then insert into Eq. (48). Next, by identifying to the time average of Eq. (31), we can define:

$$\begin{array}{cc}& \frac{\mathrm{d}\stackrel{\sim }{\mathrm{\Theta }}}{\mathrm{d}t}{|}_{\mathrm{LT}}\equiv -{\omega }_{\mathrm{LT}}sin\theta sin\psi ,\hfill \end{array}$$(51)

$$\begin{array}{cc}& \frac{\mathrm{d}\stackrel{\sim }{\mathrm{\Xi }}}{\mathrm{d}t}{|}_{\mathrm{LT}}\equiv -{\omega }_{\mathrm{LT}}sin\theta cos\psi ,\hfill \end{array}$$(52)

and thus get

$$\begin{array}{cc}& \mathrm{\Delta }{\stackrel{\sim }{\mathrm{\Theta }}}_{\mathrm{LT}}=P\frac{\mathrm{d}\stackrel{\sim }{\mathrm{\Theta }}}{\mathrm{d}t}{|}_{\mathrm{LT}}=\frac{2\pi }{n}\frac{\mathrm{d}\stackrel{\sim }{\mathrm{\Theta }}}{\mathrm{d}t}{|}_{\mathrm{LT}}=-4\pi {ϵ}^{3/2}\chi sin\theta sin\psi ,\hfill \end{array}$$(53)

$$\begin{array}{cc}& \mathrm{\Delta }{\stackrel{\sim }{\mathrm{\Xi }}}_{\mathrm{LT}}=P\frac{\mathrm{d}\stackrel{\sim }{\mathrm{\Xi }}}{\mathrm{d}t}{|}_{\mathrm{LT}}=\frac{2\pi }{n}\frac{\mathrm{d}\stackrel{\sim }{\mathrm{\Xi }}}{\mathrm{d}t}{|}_{\mathrm{LT}}=-4\pi {ϵ}^{3/2}\chi sin\theta cos\psi .\hfill \end{array}$$(54)

This shows that $\frac{\mathrm{d}\stackrel{\sim }{\mathrm{\Theta }}}{\mathrm{d}t}$ and $\frac{\mathrm{d}\stackrel{\sim }{\mathrm{\Xi }}}{\mathrm{d}t}$ are the time-average of $\frac{\mathrm{d}\mathrm{\Theta }}{\mathrm{d}t}$ and $\frac{\mathrm{d}\mathrm{\Xi }}{\mathrm{d}t}$ and that $\mathrm{\Delta }{\stackrel{\sim }{\mathrm{\Theta }}}_{\mathrm{LT}}=\mathrm{\Delta }{\mathrm{\Theta }}_{\mathrm{LT}}$ and $\mathrm{\Delta }{\stackrel{\sim }{\mathrm{\Xi }}}_{\mathrm{LT}}=\mathrm{\Delta }{\mathrm{\Xi }}_{\mathrm{LT}}$. In other words, we end up with the same results given by the instantaneous formalism used in Sect. 4.1, which links the two aspects of instantaneous and averaged precession shifts.

If we do the same exercise for the quadrupole moment using Eqs. (31) and (C.9), we get

$$\begin{array}{c}\hfill {\langle \frac{\mathrm{d}{\mathbf{z}}_{\mathrm{orb}}}{\mathrm{d}t}{|}_Q\rangle }_t=-{\omega }_Q({\mathbf{z}}_{\mathrm{bh}}·{\mathbf{z}}_{\mathrm{orb}})({\mathbf{z}}_{\mathrm{bh}}\times {\mathbf{z}}_{\mathrm{orb}}),\end{array}$$(55)

where

$$\begin{array}{c}\hfill {\omega }_Q=\frac{3n}{2}{ϵ}^2{\chi }^2,\end{array}$$(56)

Once again, we see that $\mathrm{\Delta }{\stackrel{\sim }{\mathrm{\Theta }}}_Q=\mathrm{\Delta }{\mathrm{\Theta }}_Q$ and $\mathrm{\Delta }{\stackrel{\sim }{\mathrm{\Xi }}}_Q=\mathrm{\Delta }{\mathrm{\Xi }}_Q$ as with the the LT effect discussed above.

We see from Eqs. (48) and (55) that z_orb has an average precession around z_bh, meaning that the angle θ between z_orb and z_bh has no secular evolution. In other words, the secular out-of-plane precession is done at fixed¹⁸θ, of z_orb about z_bh (as illustrated in Fig. 3), and Eqs. (49) and (56) are the frequencies at which the out-of-plane precession spans the cone of angle θ represented in Fig. 3, due to the LT (at the 1.5PN order) and quadrupole moment (at the 2PN order) effects, respectively. We also see that the direction of secular rotation of z_orb due to the LT effect is the same as that of the spin. The direction of secular rotation of z_orb due to the quadrupole moment effect depends on the sign of −z_bh ⋅ z_orb = −cos θ, as shown by Eq. (55); in the case of θ > π/2 (a retrograde orbit relative to the black hole), the direction of rotation is the same as that of the spin, for example. However, the opposite is true in the case of prograde orbits relative to the black hole.

In conclusion, we stress that both the in-plane and out-of-plane precessions can be considered differently on two timescales: time-averaged (on the cone defined by Eqs. (48) and (55) for the out-of-plane precession) and instantaneously, which corresponds to a rotation of the axes around each other, given by Eqs. (29) and (31). We note that it can also deviate from the cone defined by Eqs. (48) and (55) for the out-of-plane precession.

4.4. Numerical illustration

We illustrate the instantaneous variations in the orientation of the orbital frame using numerical integrations of the PN Lagrange planetary equations (see Sect. 3.2). These simulations were tailored to maximize the impact of each relativistic effect (see the analysis of Sect. 4.1.2). To do so, we aligned z_bh on each axis of the orbit frame, z_bh//z_orb (maximizing the in-plane precession of all effects), z_bh//x_orb and z_bh//y_orb (maximizing the out-of-plane precession of the LT effect); these three configurations are referred to as “equatorial orbit”, “polar orbit with polar major axis”, and “polar orbit with equatorial major axis”, respectively. Finally, we refer to Appendix E for a discussion of the study of “mid-inclined orbits”, which maximize the out-of-plane precession of the quadrupole moment, where ${\mathbf{z}}_{\mathrm{bh}}=-{\mathbf{y}}_{\mathrm{orb}}^{\mathrm{pro}}+{\mathbf{z}}_{\mathrm{orb}}^{\mathrm{pro}}={\mathbf{y}}_{\mathrm{orb}}^{\mathrm{ret}}-{\mathbf{z}}_{\mathrm{orb}}^{\mathrm{ret}}$. This provides an intuitive approach to the way the spin of the black hole can impact the orbit, independently of the observer.

The simulated star is taken to have orbital parameters identical to “S2/10” in terms of the semimajor axis and eccentricity. We performed a Euclidean projection of the orbit onto a plane, noting that this was for illustrative purposes only and does not aim at reproducing observables. Effects such as the Rømer delay, gravitational redshift or lensing were deliberately omitted, as they fall outside the scope of this analysis. In the following subsections, we show that the simulated variations in these configurations reproduce the expected PN trends, both in direction and in structure. We defined a prograde and retrograde orbit relative to the black hole, with angular momentum along ${\mathbf{z}}_{\mathrm{orb}}^{\mathrm{pro}}$ and ${\mathbf{z}}_{\mathrm{orb}}^{\mathrm{ret}}$, as orbits verifying ${\mathbf{z}}_{\mathrm{orb}}^{\mathrm{pro}}·{\mathbf{z}}_{\mathrm{bh}}>0$ and ${\mathbf{z}}_{\mathrm{orb}}^{\mathrm{ret}}·{\mathbf{z}}_{\mathrm{bh}}<0$, respectively.

4.4.1. In-plane precession: Focusing on equatorial orbits

In this subsection, we examine the impact of the in-plane precession due to Kerr effects, independently of the projection on the screen of a distant observer. For this purpose, we simulated orbits that are in the equatorial plane of the black hole, thereby maximizing the in-plane precession and eliminating the out-of-plane ones¹⁹.

We considered two orbits in the equatorial plane of the black hole: one prograde (θ = 0) with $\mathrm{\Delta }{\varpi }_{\mathrm{Kerr}}^{\mathrm{pro}}=\mathrm{\Delta }{\varpi }_{\mathrm{LT}}^{\mathrm{pro}}+\mathrm{\Delta }{\varpi }_Q^{\mathrm{pro}}$ and the other retrograde (θ = π) with $\mathrm{\Delta }{\varpi }_{\mathrm{Kerr}}^{\mathrm{ret}}=\mathrm{\Delta }{\varpi }_{\mathrm{LT}}^{\mathrm{ret}}+\mathrm{\Delta }{\varpi }_Q^{\mathrm{ret}}$, making ψ degenerate. For understanding the in-plane precession, in Figs. 4 and 5, we offer a Euclidean projection of these orbits onto a face-on plane.

Let us start with the well-known Schwarzschild precession. We know that ω is measured in the same direction as the movement of the star in its orbit. When looking at the analytical expressions, we see that Eq. (38) is independent from θ, it leads to the same positive Δϖ_sch for prograde and retrograde orbits, regardless of the orientation of the black hole with respect to the orbit. This means that the prograde and retrograde apocenter displacements will go in opposite ways. We see in Figs. 4b and 5b that the Schwarzschild precession shifts the apocenters of the prograde and retrograde orbits counterclockwise and clockwise, respectively. In simpler words, the simulations agree with Eq. (38) in terms of orientation on the projection plane.

Now, let us examine the LT effect. Figure 4b indicates an apocenter shift opposite to the direction of the spin and the same is observed in Fig. 5b. Indeed, when looking at Eq. (39), we see that for the prograde orbit, $\mathrm{\Delta }{\varpi }_{\mathrm{LT}}^{\mathrm{pro}}<0$ and $\mathrm{\Delta }{\varpi }_{\mathrm{LT}}^{\mathrm{ret}}>0$, meaning that both the prograde and retrograde apocenters will be shifted clockwise (i.e., opposite to the direction of the spin).

Finally, we see from Figs. 4b and 5b that, similarly to the Schwarzschild precession, the quadrupole moment shifts the apocenters in the same direction as the movement of the star in its orbit. Indeed, when looking at Eq. (40), we see that the sign of Δϖ_Q depends only on the value of cos²θ meaning that Δϖ_Q^pro = Δϖ_Q^ret > 0.

Finally, it is important to remember that all these effects do not share the same dominant PN order and, therefore, they do not act on the same scale. This is evident in Figs. 4b and 5b, where the Schwarzschild precession dominates the LT effect, which, in turn, dominates the quadrupole moment effect. Also, it is interesting to note that the total secular shift of orbits in the equatorial plane of a rotating black hole will differ in absolute value in the prograde and retrograde configurations; in the first case, the LT and quadrupole moment effects cancel each other out, whereas in the second case, they push the apocenter in the same direction.

If we summarize all the above, we get the following statements illustrated in Fig. 6:

1. The Schwarzschild precession pushes the apocenter in the same direction as the movement of the star in its orbit.

2. The secular shift due to the LT effect pushes the apocenter of an equatorial orbit in the direction opposite to the black hole’s rotation.

3. The secular shift due to the quadrupole moment pushes the apocenter of an equatorial orbit in the same direction as the movement of the star in its orbit.

4. The analytical total Kerr secular shift of equatorial orbits verifies $|\mathrm{\Delta }{\varpi }_{\mathrm{Kerr}}^{\mathrm{pro}}|<|\mathrm{\Delta }{\varpi }_{\mathrm{Kerr}}^{\mathrm{ret}}|$.

4.4.2. Lense-Thirring out-of-plane precession

In this section, we want to illustrate the configurations that maximize the dominant spin-induced effect brought on by the LT precession. As suggested in Sect. 4.1.2, to better understand the out-of-plane precession of the LT and quadrupole moment effects, it is useful to split it into two subtypes, precession around the major axis of the orbit and precession around the minor axis, encoded by Eqs. (33) and (34), respectively. Therefore, we focus on two types of polar orbits, for which the out-of-plane component of this effect is strongest (see Sect. 4.1):

1. Polar orbit with polar major axis²⁰: z_bh//x_orb i.e., (θ, ψ) = (π/2, 0);

2. Polar orbit with equatorial major axis: z_bh//y_orb i.e., (θ, ψ) = (π/2, π/2).

To visualize the out-of-plane components independently of the projection effects, we performed an edge-on Euclidean projection, considering that the osculating orbit is a vertical line with the apocenter at the bottom of the projection plane (see Figs. 7a and 8a).

We start with the first subtype, illustrated in Fig. 7, where the orbit lies in the polar plane of the BH and verifies z_bh//x_orb. As opposed to the equatorial orbits of Sect. 4.4.1, here the notion of prograde and retrograde orbits is more delicate. Polar orbits cannot be categorized as prograde or retrograde relative to the black hole and, furthermore, an orbit with an edge-on view cannot be categorized as prograde or retrograde relative to the observer. Therefore, we instead considered an orbit where the star is moving towards the observer at the apocenter (as illustrated in Fig. 7a).

We see from Sect. 4.1.2 that when assuming (θ, ψ) = (π/2, 0), the LT effect maximally contributes to the out-of-plane precession, while it does not contribute to the in-plane one; interestingly enough, the quadrupole moment has the opposite tendency. We recall that having Δϖ_LT = 0, ΔΘ_LT = 0 and ΔΞ_LT ≠ 0 means that, instantaneously, the apocenter will not experience any secular precession attributed to the LT effect. However, as soon as the apocenter is shifted out of the spin axis due to other effects, the LT effect will also start contributing to the apocenter displacement. Indeed, from Fig. 7b, we can see that the LT effect consists of an out-of-plane precession that shifts the minor axis in the same direction as the spin of the black hole, gradually turning the orbit plane from being edge-on to face-on, having a counterclockwise rotation in the projection plane. Conversely, the quadrupole moment effect is in-plane and shifts the apocenters in the direction opposite to the movement of the star in its orbit. All this is in agreement with Eqs. (42)–(46).

Now, we go on to consider the second subtype of out-of-plane precession: the one around the minor axis of the orbit. To study this type of precession we take a polar orbit but this time with its minor axis being along the spin axis of the black hole (z_bh//y_orb), as illustrated in Fig. 8. Since we still have θ = π/2 in this configuration, we observe the same tendencies for the evolution of Δϖ. However, now we have ψ = π/2, Eqs. (42) and (43) yield ΔΘ ≠ 0 meaning that, as opposed to the previous case, we observe a direct out-of-plane apocenter displacement. Conversely, Eqs. (45) and (46) yield ΔΞ = 0, which means that there is no precession about the major axis and this is naturally because the spin is along the minor axis. Indeed, in Fig. 8b, we see that the LT effect shifts the major axis in the same direction as the spin of the black hole, while keeping the orbit plane edge-on this time. Similarly to the polar orbit with a polar major axis configuration described above, the quadrupole moment effect is in-plane and shifts the apocenters in the direction opposite to the movement of the star in its orbit. Once more, this is all in agreement with Eqs. (42)–(46).

Overall, this investigation demonstrates that although, in principle, a single randomly oriented orbit contains signatures of all relativistic effects, the relative contributions are often entangled. Observing orbits with diverse orientations provides complementary information and helps to disentangle the in-plane and out-of-plane precession components.

5. Summary and conclusions

In this paper, we present the various kinds of precessions undergone by a stellar orbit due to the 1.5- and 2PN-order spin effects. We show that on an orbital timescale, the orbit frame is subject to an in-plane precession parametrized by the rate $\dot{\varpi }$ and to out-of-plane precessions parametrized by the rates $\dot{\mathrm{\Xi }}$ and $\dot{\mathrm{\Theta }}$. Integrating over an orbit, we provide the 2PN expressions of Δϖ, ΔΞ, and ΔΘ, given in Eqs. (38)–(46), which correspond to the tilt angle between two successive pericenters, and to the out-of-plane tilt of the minor and major axes, respectively (see Fig. 3). We also highlight the contributions due to Schwarzschild, LT, and quadrupole effects in these formulas. Moreover, we show how these orbital timescale precessions are linked to the well-known secular timescale precession of the orbital angular momentum around the black hole spin axis. By considering orbital orientations that maximize the various kinds of precessions of the orbit, we have been able to illustrate these findings numerically and show the consistency of our PN orbit-integration code OOGRE with our analytical predictions. We demonstrate that no single orbital configuration is optimal for detecting all relativistic effects simultaneously; instead, a diversity of orbital orientations is essential for disentangling the different types of precessions. We also provide estimates of the on-sky astrometric effects associated to these various precessions for the particular case of the S2 star, as well as for “S2/10”, a hypothetical star on an orbit that is ten times smaller than that of S2. Our results, highlighting the orbital-timescale reorientation of stellar orbits due to the black hole spin, will be very relevant when we are attempting to detect the spin parameter of Sgr A* on S stars at the Galactic center through their astrometric signatures by the GRAVITY instrument. When comparing the Schwarzschild, LT, and quadrupole moment contributions, their relative magnitudes are governed by their PN order:

• Schwarzschild precession scales as ϵ ∼ 𝒪(v²/c²) (1PN);

• LT precession as ϵ^3/2 ∼ 𝒪(v³/c³) (1.5PN);

• Quadrupole moment precession as ϵ² ∼ 𝒪(v⁴/c⁴) (2PN).

Here, ϵ = GM/c²p is the small PN parameter. To quantify how these effects change for stars on tighter orbits, we considered “S2/10”. In this case, the PN parameter increases by a factor of 10 (i.e., ϵ^S2/10 ≈ 10ϵ^S2). Thus, we saw that all relativistic effects grow significantly stronger at smaller orbital radii and the shorter orbital periods of such stars enable more frequent pericenter and apocenter passages. As a result, secular precessions accumulate more rapidly over a fixed observational time span. For instance, when working with “S2/10” over one period of S2, the angular shifts are increased by a factor ≈300 for the Schwarzschild precession, by a factor 10³ for the LT precessions, and by a factor ≈3000 for the quadrupole moment shifts, compared to the star S2. These scaling properties underscore the high scientific value of monitoring stars on tighter orbits, thus offering a powerful probe of both the spin and multipolar structure of Sgr A*.

While this study concentrates on astrometric effects, relativistic dynamics also imprint measurable signatures on radial velocity curves, especially near pericenter passages where Doppler shifts are strongest. Future high-precision instruments such as ERIS and MICADO will enhance this spectroscopic leverage. Nevertheless, astrometry remains essential, as it provides sensitivity throughout the entire orbit and is less dependent on observing stars precisely at pericenter. A combined astrometric–spectroscopic approach will ultimately deliver the most robust constraints on the black hole’s spin parameters.

Our study highlights the importance of identifying and tracking new S stars in the inner regions of the Galactic Center, where relativistic effects are strongest. The next-generation GRAVITY+ (GRAVITY Collaboration 2023) and ERIS (Kravchenko et al. 2023) instruments will significantly enhance the detection and monitoring capabilities, making it possible to refine the orientation of Sgr A* and further test the no-hair theorem. Even if closer-in stars cannot be detected, an alternative approach is to use multi-star fitting with the currently known S stars. Since these stars have different orientations and thus exhibit different types of secular shifts, they can provide additional constraints on both the black hole’s spin magnitude and its orientation. This method has the potential to reduce the required monitoring time needed to test the no-hair theorem. Additionally, a deeper understanding of systematic uncertainties, including the role of stellar perturbations and external mass distributions, will be necessary to interpret high-precision orbital data robustly and provide a more comprehensive picture of the astrophysical environment surrounding Sgr A*.

Acknowledgments

For fruitful discussions we thank our colleagues Laura Bernard, Éric Gourgoulhon, Aurélien Hees and Alexandre Le Tiec, as well as Clifford M. Will.

References

---

Appendix A: Orbit parametrization

Let us present the different ways to parameterize our orbits.

A.1. Newtonian orbits

A Newtonian orbit is defined by six Keplerian parameters that correspond to the 6 Newtonian degrees of freedom (three positions + three velocities). The typical set of Keplerian parameters is the following:

• a_sma the semimajor axis;

• e the eccentricity of the orbit;

• ι ∈ [0; π] is the inclination²¹ between the fundamental and orbital planes, and encodes the rotation about ℒ_orb//sky, from Z to z_orb;

• Ω ∈ [0; 2π] is the position angle of the line of nodes ℒ_orb//sky, intersection between the orbital and fundamental planes; it encodes the rotation about Z, from X to ℒ_orb//sky;

• ω ∈ [0; 2π] is the angular position of the pericenter within the orbital plane, counted from the line of nodes; it encodes the rotation about z_orb, from ℒ_orb//sky to x_orb;

• t_peri the time of pericenter passage.

The first two parameters are sufficient to define (i) geometrically the ellipse of the orbit (that is, the semimajor and semiminor axes), and (ii) the energetics of the orbit, that is, the total energy E and angular momentum L of the two-body system. There is a one-to-one correspondence between (a_sma, e) and (E, L). Note that these parameters are independent of the choice of fundamental frame. In addition, we defined²² the three angles ι, Ω, and ω (see Fig. 1) that encode the orientation of the ellipse; they are the three Euler angles, allowing to rotate the orbit frame into the observer’s frame. We highlight this dependency of the three angular orbital parameters on the choice of fundamental plane. For instance, Will & Maitra (2017) used the black hole equatorial plane as their fundamental plane, while we use the plane of the sky. As a consequence, the inclination angle of our work is not the same as in theirs, and we can see that this has obviously profound consequences on the evolution equations that we end up obtaining. Moreover, we need the initial condition along the orbit, so we introduce t_peri the time of pericenter passage. Finally, to characterize the position of the star along the orbit, we use the true anomaly f (see Fig. 1) or the azimuthal angle φ = ω + f.

A.2. Relativistic orbits

For orbits considered in the framework of GR or PN theories, the orbital elements become time-varying. The values of the orbital parameters at a particular date can thus be seen as the orbital parameters of a Newtonian orbit, osculating the relativistic orbit at this particular date. It is customary to represent a relativistic orbit by providing the set of orbital parameters at a particular time, t_osc, when it is osculated by the Newtonian orbit described by this same set of orbital parameters. A relativistic orbit is thus described by the six usual orbital parameters, plus the osculating time t_osc. The time-varying orbital parameters of a relativistic orbit can therefore be seen as the orbital parameters of a family of osculating Newtonian elliptical orbits. From this point onward, we consider (a_sma, e, ι, Ω, ω) as functions of time such that (a_sma(t_osc), e(t_osc), ι(t_osc), Ω(t_osc), ω(t_osc)) = (a_{sma,  osc}, e_osc, ι_osc, Ω_osc, ω_osc).

As for the osculating t_peri, we take the most recent time of pericenter passage before the first observation date of each star (see table D.1 of GRAVITY Collaboration 2022, for the first observation dates).

Moreover, it is worth mentioning that adding a certain time T to t_osc does not mean that the orbit would be temporally translated by T; taking t_osc = t₀ + T instead of t_osc = t₀ does not mean that we would observe the same orbit if we wait for T. This would only be true if we also add T to t_peri as well. Since it is common not to find in literature²³ any mention of the osculating time when giving the numerical values of the orbital parameters of S stars, we could possibly end up misinterpreting the true values of these parameters (see Heißel et al., in prep.).

Finally, it is very important to realize that the same set of initial conditions will yield different GR orbits when using different coordinate systems. This is well known in the context of relativity theory and stems from the fact that orbital elements are not covariant quantities (Brumberg 2017). Therefore, a relativistic orbit should be described by a set of 6 orbital parameters, plus the osculating time, plus the mention of the coordinate system used in the integration. This information is necessary and sufficient for uniquely pinpointing the orbit (see Heißel et al., in prep.).

Appendix B: Details on the frames of reference

In this appendix, we gather some details about the frames of reference of Sect. 2, that are important for the comprehension and reproducibility of our results. First, concerning the orbit frame represented in Fig. 1, we define the line of nodes as the intersection between the orbit and sky planes, and its director vector ℒ_orb//sky goes from the black hole to 𝒜𝒩_orb//sky, the ascending node of the orbit with respect to the observer. The latter is the intersection between the sky plane and trajectory where the star passes from behind the sky plane (Z < 0) to the front (Z > 0).

Second, concerning the black hole frame, we mentioned that the choice of the orientation of x_bh and y_bh inside the black hole equatorial plane has no impact on the physics of the problem.

We have defined in Fig. 2 the orientation of z_bh in the orbit frame by introducing the angles (θ, β). These angles are very useful for simplifying the computations of the quantities appearing in Lagrange planetary equations. Here we introduce an alternative definition with the angles²⁴ (i′,Ω′) defined in Fig. B.1. These angles are more adapted when dealing with multi-star fit because they allow us to characterize the orientation of the black hole with respect to the observer we use the angles (i′,Ω′). This is particularly useful in the case of studying the spin of SgrA*; multiple stars would be used when constraining the spin and its orientation, meaning that defining the latter with respect to the orbit of one single star instead of the sky is unpractical. On the other hand, for the calculations presented in Sect. 3.2, the use of (θ, β), which define the orientation of the spin vector relative to the star’s orbit and ℒ_orb//sky, would be more practical. This allows us to define the black hole’s angles with respect to the sky plane as follows:

• i′∈[0; π] is the rotation angle from Z to z_bh;

• Ω′∈[0; 2π] is the rotation angle about Z, from X to z_bh//sky.

Third, we introduced in Fig. 1 the Gaussian frame (n_orb, m_orb, z_orb)²⁵. This allows us to write the star’s velocity vector v of norm v as

$$\begin{array}{c}\hfill \mathit{v}=v_r{\mathbf{n}}_{\mathrm{orb}}+v_t{\mathbf{m}}_{\mathrm{orb}},\end{array}$$(B.1)

with

$$\begin{array}{c}\hfill v_r=\sqrt{\frac{\mathit{Gm}}{p}}esinf;v_t=\sqrt{\frac{\mathit{Gm}}{p}}(1+ecosf).\end{array}$$(B.2)

This can be done, even in the presence of an "out-of-plane precession" (see Sect. 3) since at any time, we can consider the osculating Keplerian orbit (see Sect. 3) which contains the velocity of the star at that time. In other words, the velocity is always within the momentarily orbital plane, and thus never has a z_orb-component. In addition to the true anomaly f introduced in Sect. 2, it is also useful to define the angle φ = ω + φ ∈ [0; 2π], defined as the rotation angle about z_orb, from ℒ_orb//sky to n_orb, and characterizing the position of the star along the orbit relative to the ascending note (see Fig. 1). We know from Eq. (4.55) of Merritt (2013) that by defining this quantity we can write

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathbf{n}}_{\mathrm{orb}}=& (cos\phi cos\mathrm{\Omega }-cos\iota sin\phi sin\mathrm{\Omega })\mathbf{X}\hfill \\ \hfill & +(cos\phi sin\mathrm{\Omega }+cos\iota sin\phi cos\mathrm{\Omega })\mathbf{Y}\hfill \\ \hfill & +sin\iota sin\phi \mathbf{Z};\hfill \\ \hfill {\mathbf{m}}_{\mathrm{orb}}=& (-sin\phi cos\mathrm{\Omega }-cos\iota cos\phi sin\mathrm{\Omega })\mathbf{X}\hfill \\ \hfill & +(-sin\phi sin\mathrm{\Omega }+cos\iota cos\phi cos\mathrm{\Omega })\mathbf{Y}\hfill \\ \hfill & +sin\iota cos\phi \mathbf{Z};\hfill \\ \hfill {\mathbf{z}}_{\mathrm{orb}}=& sin\iota sin\mathrm{\Omega }\mathbf{X}\hfill \\ \hfill & -sin\iota cos\mathrm{\Omega }\mathbf{Y}\hfill \\ \hfill & +cos\iota \mathbf{Z}.\hfill \end{array}\end{array}$$(B.3)

Similarly to the expression of z_orb, we can write z_bh in the fundamental frame of the observer as

$$\begin{array}{c}\hfill {\mathbf{z}}_{\mathrm{bh}}=sini\prime cos\mathrm{\Omega }\prime \mathbf{X}+sini\prime sin\mathrm{\Omega }\prime \mathbf{Y}+cosi\prime \mathbf{Z}.\end{array}$$(B.4)

Therefore, we obtain

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathbf{z}}_{\mathrm{bh}}·{\mathbf{n}}_{\mathrm{orb}}=& cosi\prime sin\phi sin\iota \hfill \\ \hfill & +sini\prime (cos\phi cos(\mathrm{\Omega }-\mathrm{\Omega }\prime )-sin\phi cos\iota sin(\mathrm{\Omega }-\mathrm{\Omega }\prime ))\hfill \\ \hfill =& sin\theta cos(\beta -\phi );\hfill \\ \hfill {\mathbf{z}}_{\mathrm{bh}}·{\mathbf{m}}_{\mathrm{orb}}=& cosi\prime cos\phi sin\iota \hfill \\ \hfill & -sini\prime (sin\phi cos(\mathrm{\Omega }-\mathrm{\Omega }\prime )+cos\phi cos\iota sin(\mathrm{\Omega }-\mathrm{\Omega }\prime ))\hfill \\ \hfill =& sin\theta sin(\beta -\phi );\hfill \\ \hfill {\mathbf{z}}_{\mathrm{bh}}·{\mathbf{z}}_{\mathrm{orb}}=& cosi\prime cos\iota +sini\prime sin\iota sin(\mathrm{\Omega }-\mathrm{\Omega }\prime )\hfill \\ \hfill =& cos\theta .\hfill \end{array}\end{array}$$(B.5)

These expressions are very handy for the derivation of Eqs. (C.1) to (C.9).

Appendix C: Perturbations in the Gaussian frame

Let us express²⁶ the components of the various PN accelerations in the Gaussian frame (n_orb, m_orb, z_orb):

$$\begin{array}{cc}\hfill {\mathcal{R}}_{\mathrm{Sch}}^{2\mathrm{PN}}=& {\mathbf{a}}_{\mathrm{Sch}}^{2\mathrm{PN}}·{\mathbf{n}}_{\mathrm{orb}}\hfill \\ \hfill =& ϵ\frac{\mathit{Gm}}{p^2}{(1+ecosf)}^2(3(e^2+1)+2ecosf-4e^2{cos}^{2}f)\hfill \\ \hfill & -9{ϵ}^2\frac{\mathit{Gm}}{p^2}{(1+ecosf)}^4;\hfill \end{array}$$(C.1)

$$\begin{array}{cc}\hfill {\mathcal{R}}_{\mathrm{LT}}^{1.5\mathrm{PN}}=& {\mathbf{a}}_{\mathrm{LT}}^{1.5\mathrm{PN}}·{\mathbf{n}}_{\mathrm{orb}}=2{ϵ}^{3/2}\frac{\mathit{Gm}}{p^2}\chi {(1+ecosf)}^{4}cos\theta ;\hfill \end{array}$$(C.2)

$$\begin{array}{cc}\hfill {\mathcal{R}}_Q^{2\mathrm{PN}}=& {\mathbf{a}}_Q^{2\mathrm{PN}}·{\mathbf{n}}_{\mathrm{orb}}\hfill \\ \hfill =& -\frac{3}{2}{ϵ}^2\frac{\mathit{Gm}}{p^2}{\chi }^2{(1+ecosf)}^4(1-3{sin}^2\theta {cos}^2(\beta -\phi ));\hfill \end{array}$$(C.3)

$$\begin{array}{cc}\hfill {\mathcal{S}}_{\mathrm{Sch}}^{2\mathrm{PN}}=& {\mathbf{a}}_{\mathrm{Sch}}^{2\mathrm{PN}}·{\mathbf{m}}_{\mathrm{orb}}\hfill \\ \hfill =& 4ϵ\frac{\mathit{Gm}}{p^2}{(1+ecosf)}^{3}esinf-2{ϵ}^2\frac{\mathit{Gm}}{p^2}e{(1+ecosf)}^{4}sinf;\hfill \end{array}$$(C.4)

$$\begin{array}{cc}\hfill {\mathcal{S}}_{\mathrm{LT}}^{1.5\mathrm{PN}}=& {\mathbf{a}}_{\mathrm{LT}}^{1.5\mathrm{PN}}·{\mathbf{m}}_{\mathrm{orb}}=-2{ϵ}^{3/2}\frac{\mathit{Gm}}{p^2}\chi esinf{(1+ecosf)}^{3}cos\theta ;\hfill \end{array}$$(C.5)

$$\begin{array}{cc}\hfill {\mathcal{S}}_Q^{2\mathrm{PN}}=& {\mathbf{a}}_Q^{2\mathrm{PN}}·{\mathbf{m}}_{\mathrm{orb}}\hfill \\ \hfill =& -3{ϵ}^2\frac{\mathit{Gm}}{p^2}{\chi }^2{(1+ecosf)}^4{sin}^2\theta cos(\beta -\phi )sin(\beta -\phi );\hfill \end{array}$$(C.6)

$$\begin{array}{cc}\hfill {\mathcal{W}}_{\mathrm{Sch}}^{2\mathrm{PN}}=& {\mathbf{a}}_{\mathrm{Sch}}^{2\mathrm{PN}}·{\mathbf{z}}_{\mathrm{orb}}=0;\hfill \end{array}$$(C.7)

$$\begin{array}{cc}\hfill {\mathcal{W}}_{\mathrm{LT}}^{1.5\mathrm{PN}}=& {\mathbf{a}}_{\mathrm{LT}}^{1.5\mathrm{PN}}·{\mathbf{z}}_{\mathrm{orb}}\hfill \\ \hfill =& 2{ϵ}^{3/2}\frac{\mathit{Gm}}{p^2}\chi {(1+ecosf)}^{3}sin\theta \hfill \\ \hfill & \times [2(1+ecosf)cos(\beta -\phi )+esinfsin(\beta -\phi )];\hfill \end{array}$$(C.8)

$$\begin{array}{cc}\hfill {\mathcal{W}}_Q^{2\mathrm{PN}}=& {\mathbf{a}}_Q^{2\mathrm{PN}}·{\mathbf{z}}_{\mathrm{orb}}\hfill \\ \hfill & =-3{ϵ}^2\frac{\mathit{Gm}}{p^2}{\chi }^2{(1+ecosf)}^{4}cos\theta sin\theta cos(\beta -\phi ).\hfill \end{array}$$(C.9)

Appendix D: Secular evolution in the PN formalism and two-timescale analysis

It is possible to derive the secular evolution of the orbital parameters at an arbitrary PN order using a perturbation method called two-timescale analysis. This section tries to pedagogically introduce the material presented in Will & Maitra (2017), Tucker & Will (2019), while taking the orbital elements as the ones defined with respect to observer’s frame and not the black hole’s frame²⁷; this allows us to obtain secular shifts of quantities that are directly observable.

D.1. Two-timescale method

Let us label generically all orbital elements by X_α, where α is an index, running typically from 1 to 5 to label the 5 standard orbital elements (p, e, ι, Ω, ω). Here, we use the orbital plane azimuthal angle φ = ω + f as the parameter, rather than time. So we consider functions X_α(φ). Equations (10) to (14) give the time evolution dX_α/dt. By using Eq. (15) we can easily get

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\phi }{\mathrm{d}t}& =\frac{\mathrm{d}\omega }{\mathrm{d}t}+\frac{\mathrm{d}f}{\mathrm{d}t}\hfill \\ \hfill & ={\left(\frac{\mathrm{d}f}{\mathrm{d}t}\right)}_{\mathrm{Kepler}}-cos\iota \frac{\mathrm{d}\mathrm{\Omega }}{\mathrm{d}t}\hfill \\ \hfill & =\sqrt{\frac{\mathit{Gm}}{p^3}}{(1+ecosf)}^2-\sqrt{\frac{p}{\mathit{Gm}}}cot\iota \frac{sin\phi }{1+ecosf}\mathcal{W},\hfill \end{array}$$(D.1)

and convert dX_α/dt to dX_α/dφ. The expressions of the various perturbing accelerations up to 2PN order are given by Eqs. (C.1) to (C.9).

Finally, using all this material, the evolution of the orbital elements is expressed as

$$\begin{array}{c}\hfill \frac{\mathrm{d}{X}_{\alpha }}{\mathrm{d}\phi }=\stackrel{\sim }{ϵ}{Q}_{\alpha }^{(1)}+{\stackrel{\sim }{ϵ}}^{3/2}{Q}_{\alpha }^{(3/2)}+{\stackrel{\sim }{ϵ}}^2{Q}_{\alpha }^{(2)},\end{array}$$(D.2)

where we express the PN small parameter²⁸ $\stackrel{\sim }{ϵ}$ as

$$\begin{array}{c}\hfill \stackrel{\sim }{ϵ}=\frac{\mathit{GM}}{\stackrel{\sim }{p}{c}^2}·\end{array}$$(D.3)

Here, we use an orbit-averaged value of the semilatus rectum, $\stackrel{\sim }{p}$, which is constant in the evolution. The functions $\stackrel{\sim }{ϵ}{Q}_{\alpha }^{(1)}$, ${\stackrel{\sim }{ϵ}}^{3/2}{Q}_{\alpha }^{(3/2)}$, and ${\stackrel{\sim }{ϵ}}^2{Q}_{\alpha }^{(2)}$, denote the 1PN, 1.5PN, and 2PN components of the evolution equation, respectively, and the Q_α⁽ⁱ⁾ are by definition 0PN quantities.

The evolution of the orbital parameters under these PN components takes the form of a long-term (compared to the orbital period) secular evolution, on top of which there are quickly-varying (at the orbital timescale) oscillations. Such a behavior can be well captured in the framework of a two-timescale analysis. Based on this understanding of what the solution should look like, we consider that X_α is no longer a function X_α(φ) of the orbital plane azimuthal angle, but rather a function of two variables, X_α(Φ, φ), where $\mathrm{\Phi }=\stackrel{\sim }{ϵ}\phi$ captures the long-term behavior of the element, while φ captures the quickly-varying oscillations. We use the trick of considering that φ and Φ are independent variables (which is only a trick, they are actually directly related by $\mathrm{\Phi }=\stackrel{\sim }{ϵ}\phi$), and derive evolution equations for X_α, now considered as a function of two independent variables Φ and φ.

Let us rewrite with more details the evolution of X_α, which now reads

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{X}_{\alpha }(\mathrm{\Phi },\phi )}{\mathrm{d}\phi }=& \stackrel{\sim }{ϵ}{Q}_{\alpha }^{(1)}(X_{\beta }(\mathrm{\Phi },\phi );\phi )+{\stackrel{\sim }{ϵ}}^{3/2}{Q}_{\alpha }^{(3/2)}(X_{\beta }(\mathrm{\Phi },\phi );\phi )\hfill \\ \hfill & +{\stackrel{\sim }{ϵ}}^2{Q}_{\alpha }^{(2)}(X_{\beta }(\mathrm{\Phi },\phi );\phi ),\hfill \end{array}$$(D.4)

where we recall that $\stackrel{\sim }{ϵ}$ does not depend on any variable quantity, and Q_α⁽¹⁾, Q_α^(3/2) and Q_α⁽²⁾ now depend on the various orbital parameters considered as functions of two variables, X_β(Φ, φ), and explicitly also on φ. Next, we look for a solution that is expressed as

$$\begin{array}{cc}\hfill X_{\alpha }(\mathrm{\Phi },\phi )=& {\stackrel{\sim }{X}}_{\alpha }(\mathrm{\Phi })+\stackrel{\sim }{ϵ}{Y}_{\alpha }^{(1)}({\stackrel{\sim }{X}}_{\beta }(\mathrm{\Phi });\phi )+{\stackrel{\sim }{ϵ}}^{3/2}{Y}_{\alpha }^{(3/2)}({\stackrel{\sim }{X}}_{\beta }(\mathrm{\Phi });\phi )\hfill \\ \hfill & +{\stackrel{\sim }{ϵ}}^2{Y}_{\alpha }^{(2)}({\stackrel{\sim }{X}}_{\beta }(\mathrm{\Phi });\phi ),\hfill \end{array}$$(D.5)

where ${\stackrel{\sim }{X}}_{\alpha }(\mathrm{\Phi })$ captures the long-term, secular evolution of the orbital element, while Y_α⁽ⁱ⁾ capture the quickly-varying evolution due to the different perturbing accelerations, the upper index i reminding the PN order of the small oscillating correction to the secular evolution. We assume that these Y_α⁽ⁱ⁾ are 2π-periodic in φ.

There is not unicity of this ansatz of solution but it is natural to consider

$$\begin{array}{c}\hfill {\stackrel{\sim }{X}}_{\alpha }(\mathrm{\Phi })=\langle X_{\alpha }(\mathrm{\Phi },\phi )\rangle ,\langle Y_{\alpha }^{(i)}({\stackrel{\sim }{X}}_{\beta }(\mathrm{\Phi }),\phi )\rangle =0,\end{array}$$(D.6)

where

$$\begin{array}{c}\hfill \langle F(\mathrm{\Phi },\phi )\rangle =\frac{1}{2\pi }{\displaystyle {\int }_0^{2\pi }}F(\mathrm{\Phi },\phi )d\phi ,\end{array}$$(D.7)

is the orbit average of some function F(Φ, φ), where we consider that Φ remains constant in the average, and that the result does not depend on the value of Φ. So we consider that ${\stackrel{\sim }{X}}_{\alpha }$ is simply the orbit average of X_α, and that the oscillating term Y_α⁽ⁱ⁾ averages to zero over one orbit. Now that we have this ansatz of the solution, we can rewrite the evolution of the orbital parameter, X_α. We express the derivative with respect to Φ as

$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}\phi }=\stackrel{\sim }{ϵ}\frac{\partial }{\partial \mathrm{\Phi }}+\frac{\partial }{\partial \phi },\end{array}$$(D.8)

and Eq. (D.4) becomes

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{X}_{\alpha }}{\mathrm{d}\phi }& =\frac{\mathrm{d}{\stackrel{\sim }{X}}_{\alpha }}{\mathrm{d}\mathrm{\Phi }}\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}\phi }\hfill \\ \hfill & +\stackrel{\sim }{ϵ}{\displaystyle \sum _{\beta }}\left(\frac{\partial Y_{\alpha }^{(1)}}{\partial {\stackrel{\sim }{X}}_{\beta }}\frac{\mathrm{d}{\stackrel{\sim }{X}}_{\beta }}{\mathrm{d}\mathrm{\Phi }}\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}\phi }\right)+\stackrel{\sim }{ϵ}\frac{\partial Y_{\alpha }^{(1)}}{\partial \phi }\hfill \\ \hfill & +{\stackrel{\sim }{ϵ}}^{3/2}{\displaystyle \sum _{\beta }}\left(\frac{\partial Y_{\alpha }^{(3/2)}}{\partial {\stackrel{\sim }{X}}_{\beta }}\frac{\mathrm{d}{\stackrel{\sim }{X}}_{\beta }}{\mathrm{d}\mathrm{\Phi }}\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}\phi }\right)+{\stackrel{\sim }{ϵ}}^{3/2}\frac{\partial Y_{\alpha }^{(3/2)}}{\partial \phi }\hfill \\ \hfill & +{\stackrel{\sim }{ϵ}}^2{\displaystyle \sum _{\beta }}\left(\frac{\partial Y_{\alpha }^{(2)}}{\partial {\stackrel{\sim }{X}}_{\beta }}\frac{\mathrm{d}{\stackrel{\sim }{X}}_{\beta }}{\mathrm{d}\mathrm{\Phi }}\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}\phi }\right)+{\stackrel{\sim }{ϵ}}^2\frac{\partial Y_{\alpha }^{(2)}}{\partial \phi }\hfill \\ \hfill & =\stackrel{\sim }{ϵ}{Q}_{\alpha }^{(1)}+{\stackrel{\sim }{ϵ}}^{3/2}{Q}_{\alpha }^{(3/2)}+{\stackrel{\sim }{ϵ}}^2{Q}_{\alpha }^{(2)}.\hfill \end{array}$$(D.9)

Using $\mathrm{d}\mathrm{\Phi }/\mathrm{d}\phi =\stackrel{\sim }{ϵ}$ and dividing by $\stackrel{\sim }{ϵ}$, we get

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{\stackrel{\sim }{X}}_{\alpha }}{\mathrm{d}\mathrm{\Phi }}& =Q_{\alpha }^{(1)}+{\stackrel{\sim }{ϵ}}^{1/2}{Q}_{\alpha }^{(3/2)}+\stackrel{\sim }{ϵ}{Q}_{\alpha }^{(2)}\hfill \\ \hfill & -\stackrel{\sim }{ϵ}{\displaystyle \sum _{\beta }}\left(\frac{\partial Y_{\alpha }^{(1)}}{\partial {\stackrel{\sim }{X}}_{\beta }}\frac{\mathrm{d}{\stackrel{\sim }{X}}_{\beta }}{\mathrm{d}\mathrm{\Phi }}\right)-{\stackrel{\sim }{ϵ}}^{3/2}{\displaystyle \sum _{\beta }}\left(\frac{\partial Y_{\alpha }^{(3/2)}}{\partial {\stackrel{\sim }{X}}_{\beta }}\frac{\mathrm{d}{\stackrel{\sim }{X}}_{\beta }}{\mathrm{d}\mathrm{\Phi }}\right)\hfill \\ \hfill & -{\stackrel{\sim }{ϵ}}^2{\displaystyle \sum _{\beta }}\left(\frac{\partial Y_{\alpha }^{(2)}}{\partial {\stackrel{\sim }{X}}_{\beta }}\frac{\mathrm{d}{\stackrel{\sim }{X}}_{\beta }}{\mathrm{d}\mathrm{\Phi }}\right)-\frac{\partial Y_{\alpha }^{(1)}}{\partial \phi }-{\stackrel{\sim }{ϵ}}^{1/2}\frac{\partial Y_{\alpha }^{(3/2)}}{\partial \phi }-\stackrel{\sim }{ϵ}\frac{\partial Y_{\alpha }^{(2)}}{\partial \phi }·\hfill \end{array}$$(D.10)

Here, we keep in mind that we have divided by $\stackrel{\sim }{ϵ}$, so a term of order kPN in dX_α/dφ is then on the order of (k − 1)PN in $\mathrm{d}{\stackrel{\sim }{X}}_{\alpha }/\mathrm{d}\mathrm{\Phi }$.

We know that

$$\begin{array}{c}\hfill \langle \frac{\partial Y_{\alpha }^{(i)}}{\partial \phi }\rangle =\frac{1}{2\pi }{\displaystyle {\int }_0^{2\pi }}\frac{\partial Y_{\alpha }^{(i)}}{\partial \phi }d\phi =\frac{1}{2\pi }{\left[Y_{\alpha }^{(i)}\right]}_0^{2\pi }=0,\end{array}$$(D.11)

considering that Y_α⁽ⁱ⁾ is 2π-periodic. Also, if we assume that Φ and φ are independent, we can write

$$\begin{array}{c}\hfill \langle \frac{\partial Y_{\alpha }^{(i)}}{\partial {\stackrel{\sim }{X}}_{\beta }}\rangle \equiv \langle Y_{\alpha ,\beta }^{(i)}\rangle =0,\end{array}$$(D.12)

where we use the standard coma notation for partial derivatives. Then,

$$\begin{array}{c}\hfill \langle \frac{\partial Y_{\alpha }^{(i)}}{\partial {\stackrel{\sim }{X}}_{\beta }}\frac{\mathrm{d}{\stackrel{\sim }{X}}_{\beta }}{\mathrm{d}\mathrm{\Phi }}\rangle =\langle \frac{\partial Y_{\alpha }^{(i)}}{\partial {\stackrel{\sim }{X}}_{\beta }}\rangle \frac{\mathrm{d}{\stackrel{\sim }{X}}_{\beta }}{\mathrm{d}\mathrm{\Phi }}=0.\end{array}$$(D.13)

and by taking the orbit average, we get

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{\stackrel{\sim }{X}}_{\alpha }}{\mathrm{d}\mathrm{\Phi }}=& \langle Q_{\alpha }^{(1)}(X_{\beta }(\mathrm{\Phi },\phi ),\phi )\rangle +{\stackrel{\sim }{ϵ}}^{1/2}\langle Q_{\alpha }^{(3/2)}(X_{\beta }(\mathrm{\Phi },\phi );\phi )\rangle \hfill \\ \hfill & +\stackrel{\sim }{ϵ}\langle Q_{\alpha }^{(2)}(X_{\beta }(\mathrm{\Phi },\phi ),\phi )\rangle .\hfill \end{array}$$(D.14)

Let us now consider the Taylor expansions:

$$\begin{array}{cc}\hfill Q_{\alpha }^{(i)}(X_{\beta };\phi )=& Q_{\alpha }^{(i)}({\stackrel{\sim }{X}}_{\beta }+\stackrel{\sim }{ϵ}{Y}_{\beta }^{(1)}+{\stackrel{\sim }{ϵ}}^{3/2}{Y}_{\beta }^{(3/2)}+{\stackrel{\sim }{ϵ}}^2{Y}_{\beta }^{(2)};\phi )\hfill \\ \hfill =& Q_{\alpha }^{(i)}({\stackrel{\sim }{X}}_{\beta };\phi )\hfill \\ \hfill & +{\displaystyle \sum _{\gamma }}(\stackrel{\sim }{ϵ}{Y}_{\gamma }^{(1)}+{\stackrel{\sim }{ϵ}}^{3/2}{Y}_{\gamma }^{(3/2)}+{\stackrel{\sim }{ϵ}}^2{Y}_{\gamma }^{(2)})Q_{\alpha ,\gamma }^{(i)}+\dots \hfill \\ \hfill =& Q_{\alpha }^{(i)}({\stackrel{\sim }{X}}_{\beta };\phi )+\stackrel{\sim }{ϵ}{\displaystyle \sum _{\gamma }}\left(Y_{\gamma }^{(1)}{Q}_{\alpha ,\gamma }^{(i)}({\stackrel{\sim }{X}}_{\beta };\phi )\right)\hfill \\ \hfill & +{\stackrel{\sim }{ϵ}}^{3/2}{\displaystyle \sum _{\gamma }}\left(Y_{\gamma }^{(3/2)}{Q}_{\alpha ,\gamma }^{(i)}({\stackrel{\sim }{X}}_{\beta };\phi )\right)+O({\stackrel{\sim }{ϵ}}^2),\hfill \end{array}$$(D.15)

where the same coma notation for partial derivative is used. The summation is done on all parameters ${\stackrel{\sim }{X}}_{\gamma }$ that appear in Q_α⁽ⁱ⁾. Note that $Q_{\alpha }^{(i)}({\stackrel{\sim }{X}}_{\beta },\phi )$ in this expression means that we consider the orbital parameters to be only secularly varying, and constant at the orbital timescale considered in the average process. So considering all terms up to 1PN order, Eq. (D.14) becomes

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{\stackrel{\sim }{X}}_{\alpha }(\mathrm{\Phi })}{\mathrm{d}\mathrm{\Phi }}=& ⟨Q_{\alpha }^{(1)}({\stackrel{\sim }{X}}_{\beta };\phi )⟩\hfill \\ \hfill & +{\stackrel{\sim }{ϵ}}^{1/2}⟨Q_{\alpha }^{(3/2)}({\stackrel{\sim }{X}}_{\beta };\phi )⟩\hfill \\ \hfill & +\stackrel{\sim }{ϵ}[{\displaystyle \sum _{\gamma }}(⟨Y_{\gamma }^{(1)}{Q}_{\alpha ,\gamma }^{(1)}({\stackrel{\sim }{X}}_{\beta };\phi )⟩)+⟨Q_{\alpha }^{(2)}({\stackrel{\sim }{X}}_{\beta };\phi )⟩].\hfill \end{array}$$(D.16)

This is the final expression of the secular evolution of the orbital parameters. This expression is in agreement with the formulas of Will & Maitra (2017), Tucker & Will (2019) provided that only Q_α⁽¹⁾ is considered²⁹. Will & Maitra (2017), Tucker & Will (2019) derive the following expression for Y_α⁽¹⁾:

$$\begin{array}{c}\hfill Y_{\alpha }^{(1)}={\displaystyle {\int }_0^{\phi }}{Q}_{\alpha }^{(1)}d\phi \prime -(\phi +\pi )\langle Q_{\alpha }^{(1)}\rangle +\langle \phi Q_{\alpha }^{(1)}\rangle .\end{array}$$(D.17)

The first term in the square brackets of Eq. (D.16) then becomes

$$\begin{array}{cc}\hfill \langle Y_{\gamma }^{(1)}{Q}_{\alpha ,\gamma }^{(1)}({\stackrel{\sim }{X}}_{\beta },\phi )\rangle =& \langle Q_{\alpha ,\gamma }^{(1)}{\displaystyle {\int }_0^{\phi }}{Q}_{\gamma }^{(1)}\rangle +\langle \phi Q_{\gamma }^{(1)}\rangle \langle Q_{\alpha ,\gamma }^{(1)}\rangle \hfill \\ \hfill & -\pi \langle Q_{\gamma }^{(1)}\rangle \langle Q_{\alpha ,\gamma }^{(1)}\rangle -\langle Q_{\gamma }^{(1)}\rangle \langle \phi Q_{\alpha ,\gamma }^{(1)}\rangle ,\hfill \end{array}$$(D.18)

so that provided that the Q_α⁽ⁱ⁾ are known for all elements, the secular evolution can be derived. Finally we can come back to the evolution equation as a function of φ by remembering that Φ and φ are actually the same quantity up to a factor $\stackrel{\sim }{ϵ}$. So we get

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{\stackrel{\sim }{X}}_{\alpha }(\phi )}{\mathrm{d}\phi }=& \stackrel{\sim }{ϵ}⟨Q_{\alpha }^{(1)}({\stackrel{\sim }{X}}_{\beta };\phi )⟩\hfill \\ \hfill & +{\stackrel{\sim }{ϵ}}^{3/2}⟨Q_{\alpha }^{(3/2)}({\stackrel{\sim }{X}}_{\beta };\phi )⟩\hfill \\ \hfill & +{\stackrel{\sim }{ϵ}}^2[{\displaystyle \sum _{\gamma }}(⟨Y_{\gamma }^{(1)}{Q}_{\alpha ,\gamma }^{(1)}({\stackrel{\sim }{X}}_{\beta };\phi )⟩)+⟨Q_{\alpha }^{(2)}({\stackrel{\sim }{X}}_{\beta };\phi )⟩],\hfill \end{array}$$(D.19)

which is the second order secular evolution of the orbital element, here considered as a function of φ only (the variable Φ was just a mathematical trick to derive the equation).

D.2. Application to the 2PN Schwarzschild shift

In this section we neglect the spin effects and consider a 2PN Schwarzschild spacetime. So in particular the Q_α^(3/2) is zero.

Let us start from the expression of dω/dt given in Eq. (14). Equation (D.1) leads to

$$\begin{array}{c}\hfill \frac{\mathrm{d}\phi }{\mathrm{d}t}=\sqrt{\frac{\mathit{GM}}{p^3}}{(1+ecosf)}^2,\end{array}$$(D.20)

as Ω is constant in Schwarzschild (the evolution of Ω depends on 𝒲, which is linked to the spin). The evolution of ω as a function of φ then reads

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\omega }{\mathrm{d}\phi }=& \frac{1}{e}\sqrt{\frac{p}{\mathit{GM}}}[-cosf\mathcal{R}+\frac{2+ecosf}{1+ecosf}sinf\mathcal{S}]\hfill \\ \hfill & \times \sqrt{\frac{p^3}{\mathit{GM}}}{(1+ecosf)}^{-2},\hfill \end{array}$$(D.21)

by combining Eqs. (10) to (14) and (D.20). Using the 2PN expressions of ℛ and 𝒮 given in Eqs. (C.1) and (C.4), we get

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\omega }{\mathrm{d}\phi }& =ϵ\frac{1}{e}[-cosf(3+2ecosf+e^2(4{sin}^{2}f-1))\hfill \\ \hfill & +4\frac{2+ecosf}{1+ecosf}sinf(1+ecosf)esinf]\hfill \\ \hfill & +{ϵ}^2\frac{1}{e}[-cosf(-9{(1+ecosf)}^2)\hfill \\ \hfill & -2\frac{2+ecosf}{1+ecosf}sinf{(1+ecosf)}^{2}esinf].\hfill \end{array}$$(D.22)

We want to replace ϵ by $\stackrel{\sim }{ϵ}=GM/\stackrel{\sim }{p}{c}^2$, for a constant value of the semilatus rectum $\stackrel{\sim }{p}$ to obtain

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\omega }{\mathrm{d}\phi }& =\stackrel{\sim }{ϵ}\frac{1}{e}\frac{\stackrel{\sim }{p}}{p}[-cosf(3+2ecosf+e^2(4{sin}^{2}f-1))\hfill \\ \hfill & +4\frac{2+ecosf}{1+ecosf}sinf(1+ecosf)esinf]\hfill \\ \hfill & +{\stackrel{\sim }{ϵ}}^2\frac{1}{e}{\left(\frac{\stackrel{\sim }{p}}{p}\right)}^2[-cosf(-9{(1+ecosf)}^2)\hfill \\ \hfill & -2\frac{2+ecosf}{1+ecosf}sinf{(1+ecosf)}^{2}esinf],\hfill \end{array}$$(D.23)

which takes the exact same form as the general expression of Eq. (D.4). The quantity in factor of $\stackrel{\sim }{ϵ}$ in the first line corresponds to Q_ω⁽¹⁾, while that in factor of ${\stackrel{\sim }{ϵ}}^2$ in the second line corresponds to Q_ω⁽²⁾. They both depend on three orbital parameters, ω itself (through f = φ − ω), p and e, as well as explicitly on φ (again through f). They are 2π-periodic in φ. Let us simplify a bit to get

$$\begin{array}{cc}\hfill Q_{\omega }^{(1)}& =\frac{\stackrel{\sim }{p}}{p}(8+(e-\frac{3}{e})cosf-10{cos}^{2}f);\hfill \\ \hfill Q_{\omega }^{(2)}& ={\left(\frac{\stackrel{\sim }{p}}{p}\right)}^2(-4+3(\frac{3}{e}-2e)cosf+2(11-e^2){cos}^{2}f\hfill \\ \hfill & +15e{cos}^{3}f+2e^2{cos}^{4}f).\hfill \end{array}$$(D.24)

We now want to compute Eq. (D.19) for ω, which reads

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\stackrel{\sim }{\omega }}{\mathrm{d}\phi }& =\stackrel{\sim }{ϵ}\langle Q_{\omega }^{(1)}(p,e,\omega ,\phi )\rangle \hfill \\ \hfill & +{\stackrel{\sim }{ϵ}}^2[\langle Y_p^{(1)}{Q}_{\omega ,p}^{(1)}(p,e,\omega ,\phi )\rangle +\langle Y_e^{(1)}{Q}_{\omega ,e}^{(1)}(p,e,\omega ,\phi )\rangle \hfill \\ \hfill & +\langle Y_{\omega }^{(1)}{Q}_{\omega ,\omega }^{(1)}(p,e,\omega ,\phi )\rangle +\langle Q_{\omega }^{(2)}(p,e,\omega ,\phi )\rangle ],\hfill \end{array}$$(D.25)

where we remind that the orbital elements appearing in the various functions Q and Y are considered to vary only secularly (see the comment following Eq. (D.15)): they are functions of Φ alone, they do not depend on φ, so they are not affected by the orbit averaging. For instance, we have ⟨e cos f⟩ = e⟨cos f⟩, which would not be true if the full dependence of e on both φ and Φ was considered.

The first and last terms of Eq. (D.25) are trivial to compute

$$\begin{array}{cc}\hfill \langle Q_{\omega }^{(1)}\rangle & =3,\hfill \\ \hfill \langle Q_{\omega }^{(2)}\rangle & =7-\frac{e^2}{4},\hfill \end{array}$$(D.26)

where we use $p=\stackrel{\sim }{p}$, which is true here given that the orbital-timescale variability of p is discarded, only its secular variation is considered.

The value of ⟨Q_ω⁽¹⁾⟩ leads to the well-known textbook result: $\mathrm{d}\stackrel{\sim }{\omega }/\mathrm{d}\phi =3\stackrel{\sim }{ϵ}$ at 1PN order. However here, we need to go to 2PN, so a we also need to consider the three ⟨YQ⟩ terms in the bracket of Eq. (D.25). These are expressed as

$$\begin{array}{cc}& \langle Y_p^{(1)}{Q}_{\omega ,p}^{(1)}\rangle +\langle Y_e^{(1)}{Q}_{\omega ,e}^{(1)}\rangle +\langle Y_{\omega }^{(1)}{Q}_{\omega ,\omega }^{(1)}\rangle =\hfill \\ \hfill & \langle Q_{\omega ,p}^{(1)}{\displaystyle {\int }_0^{\phi }}{Q}_p^{(1)}\rangle +\langle Q_{\omega ,e}^{(1)}{\displaystyle {\int }_0^{\phi }}{Q}_e^{(1)}\rangle +\langle Q_{\omega ,\omega }^{(1)}{\displaystyle {\int }_0^{\phi }}{Q}_{\omega }^{(1)}\rangle \hfill \\ \hfill & +\langle \phi Q_p^{(1)}\rangle \langle Q_{\omega ,p}^{(1)}\rangle +\langle \phi Q_e^{(1)}\rangle \langle Q_{\omega ,e}^{(1)}\rangle +\langle \phi Q_{\omega }^{(1)}\rangle \langle Q_{\omega ,\omega }^{(1)}\rangle \hfill \\ \hfill & -\pi \langle Q_p^{(1)}\rangle \langle Q_{\omega ,p}^{(1)}\rangle -\pi \langle Q_e^{(1)}\rangle \langle Q_{\omega ,e}^{(1)}\rangle -\pi \langle Q_{\omega }^{(1)}\rangle \langle Q_{\omega ,\omega }^{(1)}\rangle \hfill \\ \hfill & -\langle Q_p^{(1)}\rangle \langle \phi Q_{\omega ,p}^{(1)}\rangle -\langle Q_e^{(1)}\rangle \langle \phi Q_{\omega ,e}^{(1)}\rangle -\langle Q_{\omega }^{(1)}\rangle \langle \phi Q_{\omega ,\omega }^{(1)}\rangle ,\hfill \end{array}$$(D.27)

where we use Eq. (D.18) to express the Y_α⁽¹⁾. We thus need to compute quite a few more terms. Let us start by:

$$\begin{array}{cc}\hfill Q_{\omega ,p}^{(1)}& =-\frac{1}{p}{Q}_{\omega }^{(1)},\hfill \\ \hfill Q_{\omega ,e}^{(1)}& =cosf(1+\frac{3}{e^2}),\hfill \\ \hfill Q_{\omega ,\omega }^{(1)}& =sinf(e-\frac{3}{e}-20cosf).\hfill \end{array}$$(D.28)

The last two have zero average, because ⟨cos f⟩ = ⟨sin f⟩ = ⟨cos f sin f⟩ = 0. So, this cancels a few terms in the second and third lines of the RHS of Eq. (D.27). Following the exact same procedure as for ω, we can write dp/dφ and de/dφ, and extract their 1PN component to get

$$\begin{array}{cc}\hfill Q_p^{(1)}& =8\stackrel{\sim }{p}esinf,\hfill \\ \hfill Q_e^{(1)}& =(3+7e^2)sinf+10ecosfsinf,\hfill \end{array}$$(D.29)

both of which have zero average (no secular evolution of p and e at 1PN order). This cancels the first two terms of the last two lines in Eq. (D.27). So we are left with:

$$\begin{array}{cc}& \langle Y_p^{(1)}{Q}_{\omega ,p}^{(1)}\rangle +\langle Y_e^{(1)}{Q}_{\omega ,e}^{(1)}\rangle +\langle Y_{\omega }^{(1)}{Q}_{\omega ,\omega }^{(1)}\rangle =\hfill \\ \hfill & \langle Q_{\omega ,p}^{(1)}{\displaystyle {\int }_0^{\phi }}{Q}_p^{(1)}\rangle +\langle Q_{\omega ,e}^{(1)}{\displaystyle {\int }_0^{\phi }}{Q}_e^{(1)}\rangle +\langle Q_{\omega ,\omega }^{(1)}{\displaystyle {\int }_0^{\phi }}{Q}_{\omega }^{(1)}\rangle \hfill \\ \hfill & +\langle \phi Q_p^{(1)}\rangle \langle Q_{\omega ,p}^{(1)}\rangle -3\langle \phi Q_{\omega ,\omega }^{(1)}\rangle .\hfill \end{array}$$(D.30)

It is easy to get

$$\begin{array}{cc}\hfill \langle \phi Q_p^{(1)}\rangle & =-8\stackrel{\sim }{p}ecos\omega ,\hfill \\ \hfill \langle \phi Q_{\omega ,\omega }^{(1)}\rangle & =-cos\omega (e-\frac{3}{e})+5cos2\omega ,\hfill \end{array}$$(D.31)

using

$$\begin{array}{c}\hfill \langle \phi sinf\rangle =-cos\omega ,\langle \phi cosfsinf\rangle =-\frac{1}{4}cos2\omega .\end{array}$$(D.32)

The terms with integrals of Q_β⁽¹⁾ are a bit more cumbersome, but it is still straightforward to obtain

$$\begin{array}{cc}\hfill \langle Q_{\omega ,p}^{(1)}{\displaystyle {\int }_0^{\phi }}{Q}_p^{(1)}\rangle & =-24ecos\omega +4e(e-\frac{3}{e}),\hfill \\ \hfill \langle Q_{\omega ,e}^{(1)}{\displaystyle {\int }_0^{\phi }}{Q}_e^{(1)}\rangle & =-12-\frac{7}{2}{e}^2-\frac{9}{2e^2},\hfill \\ \hfill \langle Q_{\omega ,\omega }^{(1)}{\displaystyle {\int }_0^{\phi }}{Q}_{\omega }^{(1)}\rangle & =\frac{25}{2}+(e-\frac{3}{e})(-3cos\omega +\frac{1}{2}(e-\frac{3}{e}))\hfill \\ \hfill & +15cos2\omega .\hfill \end{array}$$(D.33)

Putting everything together, we finally arrive at

$$\begin{array}{c}\hfill \frac{\mathrm{d}\stackrel{\sim }{\omega }}{\mathrm{d}\phi }=3\stackrel{\sim }{ϵ}-\frac{3}{4}{\stackrel{\sim }{ϵ}}^2(10-{\stackrel{\sim }{e}}^2).\end{array}$$(D.34)

Equation (D.34) reproduces the 2PN result of Will & Maitra (2017), see their Eq. (3.23a). To obtain the secular pericenter shift, we must integrate over a full radial cycle (f : 0 → 2π) rather than a full azimuthal cycle (φ : 0 → 2π), since in the PN regime the pericenter itself precesses. We expand the argument of pericenter as

$$\begin{array}{c}\hfill \omega (f)={\omega }_0(f)+{ϵ}_0{\omega }_1(f)+{ϵ}_0^2{\omega }_2(f)+\mathcal{O}({ϵ}_0^3),\end{array}$$(D.35)

with ϵ₀ = GM/(p₀c²) the PN small parameter. At first post-Newtonian order, we find

$$\begin{array}{c}\hfill \phi (f)=f+\omega (f)=f+{\omega }_0+{ϵ}_0{\omega }_1(f),\end{array}$$(D.36)

so that over a full radial period,

$$\begin{array}{c}\hfill \mathrm{\Delta }f=2\pi \Rightarrow \mathrm{\Delta }\phi =2\pi +6\pi {ϵ}_0.\end{array}$$(D.37)

Combining Eq. (D.37) with Eq. (D.34), the 2PN pericenter advance in Schwarzschild spacetime reads

$$\begin{array}{c}\hfill \mathrm{\Delta }{\stackrel{\sim }{\omega }}_{\mathrm{Sch}}=6\pi \stackrel{\sim }{ϵ}+{\stackrel{\sim }{ϵ}}^2(3\pi +\frac{3}{2}\pi {\stackrel{\sim }{e}}^2),\end{array}$$(D.38)

where the substitution $ϵ\to \stackrel{\sim }{ϵ}$ introduces only 3PN-level corrections. Thus, the final expression for the 2PN secular shift of the pericenter/apocenter is

$$\begin{array}{c}\hfill \boxed{\mathrm{\Delta }{\omega }_{\mathrm{Sch}}=6\pi \stackrel{\sim }{ϵ}+\frac{3\pi }{2}{\stackrel{\sim }{ϵ}}^2(2+{\stackrel{\sim }{\mathrm{e}}}^2)}.\end{array}$$(D.39)

Here, the quantities $\stackrel{\sim }{ϵ}$ and $\stackrel{\sim }{e}$ refer to orbit-averaged parameters. If, in the 2PN term of Eq. (D.39), we replace $\stackrel{\sim }{ϵ}$ and $\stackrel{\sim }{e}$ with their osculating counterparts ϵ_osc = ϵ(p_osc) and e_osc, this substitution would only modify the result at 3PN order–meaning the 2PN expression remains unchanged. However, replacing $\stackrel{\sim }{ϵ}$ with ϵ_osc in the 1PN term would already introduce a discrepancy at 2PN order.

This is why, when expressing the pericenter advance to 2PN accuracy, it is essential to keep the orbit-averaged quantity $\stackrel{\sim }{ϵ}$ in the 1PN term ($6\pi \stackrel{\sim }{ϵ}$ and not 6πϵ_osc). Alternatively, if the preference is to work with osculating elements throughout, the entire perturbative approach–including the ansatz for the solution (e.g., Eq. (D.6))–must be adjusted accordingly. In that case, the expression for the 2PN pericenter advance would take a different form, reflecting the difference in parameterization.

In contrast, for the 2PN-order secular shifts due to the LT effect and the black hole’s quadrupole moment, the expressions involve only leading-order contributions. Therefore, we could safely use either averaged or osculating values of p, θ, and ψ without affecting the accuracy at the given order.

Appendix E: Mid-inclined orbits

We investigate the case of a mid-inclined orbit, for example, the case where ${\mathbf{z}}_{\mathrm{bh}}=-{\mathbf{y}}_{\mathrm{orb}}^{\mathrm{pro}}+{\mathbf{z}}_{\mathrm{orb}}^{\mathrm{pro}}={\mathbf{y}}_{\mathrm{orb}}^{\mathrm{ret}}-{\mathbf{z}}_{\mathrm{orb}}^{\mathrm{ret}}$, namely, (θ, ψ) = (π/4, 3π/2) or (θ, ψ) = (3π/4, π/2). we see from Eqs. (39), (42), and (45) that with these values of θ and ψ, we end up observing both in-plane and out-of-plane precessions due to the LT effect, at moderate values. However, when looking at Eqs. (43) and (46) we see that the orbits will experience the most out-of-plane shift due to the quadrupole moment in this configuration. Indeed, when simulating a prograde and retrograde orbit relative to the black hole, we see in Figs. E.1 and E.2 that the LT and quadrupole moment effects each make the apocenters experience both an in-plane and out-of-plane precession. The in-plane precessions are similar the the equatorial orbit case. As for the out-of-plane precessions, they mainly act on the major axis but also on the minor axis on much longer timescales (see Eq. (E.1) below). We also see that the quadrupole moment appears to oppose the LT effect for the out-of-plane shift when the orbit is prograde, but not for the retrograde orbit.

In practice, unless we have an equatorial (no out-of-plane precession) or polar (surface of the out-of-plane precession cone represented in Fig. 3 being a plane for θ = π/2 (mod π)) orbit with the major axis being equatorial or polar, both the major and minor axis will precess. However, the out-of-plane shift remains mainly along the tangent to the surface of the out-of-plane precession cone as long as the period of the orbit is much less than the precession period P^cone of the angular momentum of the orbit around the out-of-plane precession cone. To check that, we derive P^cone/P using Eqs. (49) and (56) at the 2PN order, expressed as

$$\begin{array}{c}\hfill \frac{P_{\mathrm{LT}}^{\mathrm{cone}}}{P}=\frac{2\pi }{{\mathrm{\omega }}_{\mathrm{LT}}P};\frac{P_Q^{\mathrm{cone}}}{P}=\frac{2\pi }{{\mathrm{\omega }}_{Q}P}\end{array}$$(E.1)

If we apply this estimation for "S2/10" for example, we get $\frac{P_{\mathrm{LT}}^{\mathrm{cone}}}{P}(\chi =0.99)\sim {10}^3$, $\frac{P_{\mathrm{LT}}^{\mathrm{cone}}}{P}(\chi =0.1)\sim {10}^4$, $\frac{P_Q^{\mathrm{cone}}}{P}(\chi =0.99)\sim {10}^5$ and $\frac{P_Q^{\mathrm{cone}}}{P}(\chi =0.1)\sim {10}^7$.

In Fig. E.3, we summarize in a schematic way the results of Figs. E.1 and E.2 and consider how the quadrupole moment out-of-plane precession of the major axis, is studied through a mid-inclined orbits with z_bh//orb//y_orb. Alternatively, in order to maximize the quadrupole moment out-of-plane precession of the minor axis instead of the major axis, we would need to study a mid-inclined orbits with z_bh//orb//x_orb.

In conclusion, we now have insight into the different types of precession that can manifest around a rotating black hole and comprehend how a random relative orientation between the orbit and black hole can be decomposed into these precession categories. With this new understanding, we can have a better understanding of the possible black hole orientations that can generate the secular precessions that observed for the orbits of currently, and potential closer-in stars.

All Tables

All Figures