Analytic Post-Newtonian Astrometric and Spectroscopic Models of Orbits around Black Holes

Observations of the S-stars, the cluster of young stars in the inner 0.1 pc of the Galactic Center, have been crucial in providing conclusive evidence for a supermassive black hole at the center of our galaxy. Since some of the stars have orbits less than that of a typical human lifetime, it is possible to observe multiple orbits and test the weak-field regime of general relativity. Current calculations of S-star orbits require slow and expensive computations in order to numerically solve geodesic equations for many small time steps. In this paper, we present a computationally efficient, first-order post-Newtonian model for the astrometric and spectroscopic data gathered for the S-stars. We find that future, 30-m class telescopes -- and potentially even current large telescopes with very high spectroscopic resolution -- may be able to detect the Shapiro effect for an S-star in the next decade or so.

One of the closest stars to Sgr A* is S0-2 (also known as S2; Schödel et al. 2002;Ghez et al. 2003;Gillessen et al. 2009), which has an orbital period of around 16 years and eccentricity of 0.88. Long-term studies of its orbit led to the first detection of gravitational effects during its 2018 periapsis -namely gravitational redshift (Gravity Collaboration et al. 2018Do et al. 2019;Gravity Collaboration et al. 2021, 2022a and Schwarzschild precession (Gravity Collaboration et al. 2020) -in an S-star orbit. As both photometric and spectroscopic sensitivities improve and shorter-period S-star candidates are identified (e.g., Peißker et al. 2020a,b), additional tools are needed to analyze and detect higher-order general relativistic effects, such as the Shapiro delay, and additional precession due to the frame dragging and quadrupole moment of the spacetime (e.g., Wex & Kopeikin 1999;Weinberg & Milosavljevic 2004;Will 2008;Merritt et al. 2010;Angélil et al. 2010;Angélil & Saha 2014;Psaltis et al. 2016;Grould et al. 2017;Waisberg et al. 2018).
Current modeling of S-star orbits involves integrating numerically the general relativistic (GR) equations of motion for each time step (e.g., Gillessen et al. 2017;Do et al. 2019). This method requires integrating across the span of observations using very small time steps to avoid error buildup and results in slow, expensive computations. Furthermore, the computational cost of such numerical calculations increases rapidly when using the orbits of multiple stars to jointly constrain the shared properties of the system (e.g., the black hole mass), since this involves simultaneously solving the geodesic equations for each time step for each star. This approach could become prohibitive when searching the multi-dimensional orbital parameter space with a statistical sampling algorithm, such as a Markov Chain Monte Carlo (MCMC) to obtain optimal solutions and quantify uncertainties in orbital parameters.
A simplification to these calculations can be introduced because of the fact that S-star orbits have pericenter distances that range from 1,400 Sgr A* Schwarzschild radii to values that are larger by orders of magnitude (Gillessen et al. 2017). At such distances, the orbits can be modeled as primarily Keplerian orbits with small corrections caused by GR effects. A framework for describing this behavior is through a semi-analytic post-Newtonian (PN) model, which uses traditional Keplerian orbital equations with additional terms up to some order in v/c, derived from GR equations. Damour & Deruelle (1985;1986, referred to hereafter as D&D I and D&D II) obtained an elegant analytic solution to the two-body problem in the first post-Newtonian order (1PN), which incorporates a variety of relativistic effects. The timing model developed by D&D II has been the workhorse for the pulsar community for many years in detecting relativistic effects (e.g., Edwards et al. 2006). It was further expanded to include second-order post-Newtonian (2PN) terms (Damour & Schafer 1988;Wex 1995).
While the D&D II analytic solution has been implemented in a timing model for fitting pulsar arrival times, the same model is not readily applicable for fitting astrometric and spectroscopic data of stars. This is because the latter relies on the Doppler shifts of emission lines in the stellar spectra as opposed to time intervals between pulses. The beauty of the analytic D&D II timing model, however, makes it possible to derive the line-of-sight velocities that correspond to a variety of time delays.
In this paper, we use the framework of the D&D I and D&D II 1PN two-body model (in harmonic coordinates) to derive a new analytic astrometric and spectroscopic model that incorporates the relativistic and astrophysical effects relevant to modeling S-star orbits.
In addition to computational efficiency, there are significant scientific advantages to our approach. In the numerical approach, all relativistic effects of the same post-Newtonian order that are embedded in the geodesic equations are reported as a single observable (i.e., the position in the sky or spectroscopic line shift). In principle, the magnitudes of the individual effects can be disentangled by exploring the differences in the numerical solutions with and without each effect (Grould et al. 2017). In contrast, a post-Newtonian analytic approach, such as our model, allows for calculating directly the characteristic "fingerprint" of each effect on the observables separately from the others (see Section 6), while providing a direct analytic handle of the dependence of each effect on the various parameters of the system.
The following sections present the two-body orbital equations ( §2), the projection of those equations to the plane of the sky ( §3), and the equations for spectroscopic effects ( §4). We discuss the implications of the model for S-star observations in Section 6. Unless otherwise indicated, we use geometrized units, i.e., G = c = 1, where G is the gravitational constant and c is the speed of light.

ORBITAL MODEL
As discussed in Section 1, our model has three main components: the two-body orbital model, the astrometric model, and the spectroscopic model. Since we want to use our model to be able to fit observations of the S stars, we first identify what parameters are the observables. With telescopes, we are able to observe the projected positions of the stars, i.e., right ascension (R.A., α) and declination (Decl., δ), and the radial velocities derived from their spectra, where λ 0 is the rest-frame wavelength of the stellar absorption or emission line used for measuring radial velocities and ∆λ is the shift at time t obs between observed and rest-frame wavelengths. Four free parameters -orbital period (P ), total mass of the system (M ), mass ratio of the two bodies (q), and radial eccentricity (e R ) govern the shape, period, and rate of precession of the two-body orbits.
The orientations of the orbits with respect to Earth determine the two-dimensional motions in the sky (i.e., the astrometric model) and the line-of-sight motions, which we can derive through spectroscopy. The three orientation angles -inclination (i), argument of ascending node (Ω), initial argument of periapsis (ω 0 ), and initial time of periapsis (or epoch of position, t 0 ) -are free parameters for our model. We must also fit the projected and line-of-sight proper motions of the Galactic center with respect to Earth (µ α , µ δ , µ ).
Nineteen additional quantities (introduced in the following sections) are derived parameters that we calculate from the observed or free parameters. Table 1 in Appendix B lists all the parameters introduced in this paper with appropriate references and units. Figure 1 illustrates how the binary system orbital parameters relate to each other.
In this section, we present the orbital equations and parameters in geometrized units.
− BB periapsis − Figure 1. The binary orbital parameters in the barycenter frame of the system. The large red, precessing ellipse shows the motion of the smaller mass m2 over roughly two periods. Similarly, the smaller blue ellipse shows the motion of the larger mass m1 over the same amount of time. Dotted magenta lines mark the semimajor and semiminor axes of the ellipses for both objects around the binary barycenter, which is indicated by a black dot. Arrows between the binary barycenter and the two masses denote the distances r1 and r2 between the binary barycenter and orbiting objects. A black x marks the closest approach that the smaller mass m2 makes to the larger mass m1, which precesses by an angle ∆θ for each orbit. The position angle θ is the angle formed between the periapsis point and the location of one of the masses, which are offset from each other by 180 • .

Coordinate Systems
We use four main coordinate systems/reference frames (Figure 2), all of which are described in Edwards et al. (2006). These are: the star reference frame (denoted by the subscript "star"), the binary barycenter (BB), the solar system barycenter (SSB), and the observer reference frame ("obs"). The different reference frames are necessary to define the various time delays discussed in Section 4. All celestial sky coordinates are given in the International Celestial Reference System (ICRS, Luzum & Petit 2015).
The times (or "clocks") we use in this paper are directly related to the four reference frames. We define the time of emission as measured at the stellar location as t star , the same time as measured by an observer at the binary barycenter as t BB , the time of arrival at the binary barycenter as t a,BB , the time of arrival at the solar system barycenter as t SSB , and the time of arrival recorded by the observer on Earth as t obs .
The relation between the observer light-arrival time t obs and the star light-emission time t star (used for the time-dependent orbital equations in Section 2.7) is the star-frame emission time plus the sum of all time delays due to binary system motion, solar system motion, and motion between the binary barycenter and solar system barycenter/observer, i.e., t obs = t star +∆ RB + ∆ EB + ∆ SB [binary effects] +∆ KB [parallax effects] +∆ R .
[solar system effects] (2) In our model, the total binary system time delay includes the binary Roemer delay (∆ RB ), binary Einstein delay (∆ EB ), and binary Shapiro delay (∆ SB ). The total solar-system-related time delay is the Earth Roemer delay (∆ R ). The interstellar time delay comprises time delays due to parallax and proper motion, which is simply the Kopeikin effect (∆ KB ). We provide explicit formulae for the time delays in Section 4.
Additional relations between the other time variables are as follows. The binary Einstein delay relates the star emission time to the time as measured at the binary barycenter by an inertial observer as t BB = t star + ∆ EB . The time of arrival at the binary barycenter relates to the star emission time via the binary effects, i.e., t a,BB = t star + ∆ RB + ∆ EB + ∆ SB . Similarly, the observer time of arrival relates to the solar system barycenter time of arrival via the solar system effects, i.e., t obs = t SSB +∆ R . Since the spectroscopic timing model is less sensitive than pulsar timing, we neglect interstellar delays (such as dispersion). As a result, the solar system barycenter time of arrival and binary barycenter time of arrival are related by a constant offset, which we set to zero without loss of generality.

Solar system barycenter (SSB)
Observer frame (obs) Star frame (star) Binary barycenter (BB) Distance ( ) between SSB and BB Figure 2. The relationship between the observer frame, the solar system barycenter, the binary barycenter, and the star frame, as defined in this paper. Distances and object sizes have been rescaled to show effect. The red axes for each reference frame are arbitrary to show how coordinate systems may vary from frame to frame depending on orientation.

Mass Measures
Three of the key values in describing a two-body orbital system are the total system mass and the individual masses of the two objects. We denote the component masses by m 1 and m 2 , such that m 2 ≤ m 1 . In the case of modeling stars orbiting a supermassive black hole, m 1 is the mass of the black hole and m 2 is the mass of the star. The total mass of the system is M = m 1 + m 2 . In classical, two-body orbits, the mass ratio of the two objects, is an important parameter that we can use to rewrite the reduced mass, and dimensionless reduced mass (ν) as Note that in Sgr A*-S-star systems, the mass ratios can be q ∼ 10 −6 − 10 −5 . In these cases, one can take the limit of q 1, but in this paper we leave the full expressions for generality.

Period, Semimajor Axis, and Mean Motion
The easiest, most direct property to measure is the orbital period, P . With the total system mass and orbital period, we use the 1PN version of Kepler's law (Blanchet et al. 1998, their Eq. 8.6) to calculate the center-of-mass semimajor axis, implicitly, via Another useful related quantity is the mean motion, the constant angular speed needed for an object to complete an equivalent circular orbit. It relates to the inverse period as n = 2π P . (4b)

Energy and Momentum
In addition to calculating the semimajor axis a R from P , we also fit for the eccentricity e R . In combination with the system mass M and the dimensionless reduced mass ν, these two parameters set the orbital behavior and we use them to calculate the energy and angular momentum of the system. Here, the total energy is and total angular momentum is (5b) We also define the quantity K, which is the general relativistic correction to the total angular momentum of the system, rewritten in natural units from D&D I Equation (4.14) and given by This is a particularly important quantity, as the value of K is what governs the orbital precession , where ∆θ is the angle the orbit precesses in each period (D&D I).

Other Eccentricities
While Keplerian orbits have only one effective eccentricity, e R , in GR, there are additional eccentricities that result from the curved spacetime (D&D I, their Eqs. 3.6b-c, 4.13). In the 1PN model, these are the time eccentricity (e t ; D&D I, their Eq. 3.6c), and the angular eccentricity (e θ ; D&D I, their Eq. 4.13), which in natural units is: In some cases, the differences between the eccentricities are negligible and all eccentricity expressions give comparable answers (see Section 4 for examples). The radial eccentricity can be determined readily from observational astrometric data.

Individual Objects Parameters
As noted in Section 2.3, the semimajor axis and radial eccentricity defined here are with respect to the center of mass of the system. Since observations of the S stars result in tracking the orbits of the individual stars, we need the derived parameters (i.e., semimajor axis and radial eccentricity) that give the orbital shapes of both objects in a two-body system, which we can derive from the corresponding effective one-body parameters (i.e., a R and e R ) and the mass ratio (q).
For the more massive of the two bodies, m 1 , the semimajor axis of its respective orbit around the binary barycenter is and the radial eccentricity is (D&D I, their Eqs. 6.3a-b, in natural units and mass ratio q). Similarly, the less massive of the two bodies, m 2 , follows an orbit around the binary barycenter with a semimajor axis of and a radial eccentricity of (D&D I, their Eqs. 6.3a-b, in natural units and mass ratio q).

Time-dependent Orbital Motion
In Sections 2.2 through 2.6, we presented the equations necessary for calculating many of the derived parameters in the model. In this section, we use those parameters to obtain the time-dependent orbital motion of the individual objects and the binary barycenter.
The heart of this time dependence comes from Kepler's equation, which relates time eccentricity (e t ), mean motion (n), mean anomaly (u), star time of emission (t star ), and epoch of position The mean anomaly, which is the angle between the periapsis of an orbit and another position in the orbit at some time, is crucial for calculating the other values in the polar orbital equations (i.e., radius and angle). The previous equation does not have an analytical solution for u but can be solved using a fast algorithm, such as the Newton-Raphson method. Note that, unlike the time-dependent equations in the astrometric model ( § 3), which use the observer time t obs , Kepler's equation is evaluated at the star time of emission t star . This difference, described by Equation 2, takes into account the vacuum retardation effect for the astrometric model.
The time-dependent distances of the two orbiting bodies from the barycenter, as given in D&D I (their Eqs. 7.1d-e), are and Calculating the position angle θ of the orbiting objects is somewhat more complicated. D&D I present the rather straightforward equation (their Eq. 4.11b) but using this in computations requires special care because the evaluation of the term tan (u/2) results in floating-point errors for the calculated value of θ around the asymptotes (i.e., u = xπ for odd x). Instead, we use a series expansion (D&D II, their Eq. 17d) of this equation (D&D I, their Eq. 4.11a), which avoids these computational issues: where We performed convergence tests of the series A e (u) at different values of angular eccentricity (e θ ) and determined that only a small number of terms is needed for necessary computational accuracy, e.g., 30 terms for fractional errors of less than 10 −8 .  Much like the position angle θ, the argument of periapsis also precesses due to GR effects (Equation 5c). Given the initial argument of periapsis, the time-dependent argument of periapsis (in terms of the mean anomaly, u) is In the previous section ( §2), the 1PN orbital equations are expressed with respect to the binary barycenter. In order to model observations, however, we must transform them to a frame with respect to the plane of the sky, i.e., in terms of right ascension (α) and declination (δ).
We transform the positions from the Cartesian binary barycentric frame to the plane of the sky using using three angles: inclination (i), longitude of the ascending node (Ω), and argument of periapsis (ω 0 ). The inclination describes the tilt of the orbit with respect to the observer, while the longitude of the ascending node is the rotation of the location of the ascending node (i.e., where the orbit intersects with the reference plane) with respect to the center of mass. Similarly, the argument of periapsis is the rotation of the object point of closest approach to the center of mass in the orbital plane. See Figure 4 for a graphic depiction.
We use the Thiele-Innes constants (A, B, C, F, G, H) defined by Equations (20)-(25) in O'Neil et al. (2019) to describe this rotation: and H = cos ω sin i .  . The angled blue ellipse shows the orbital plane, with the periapsis point marked. The grey plane is parallel to the vector pointing toward the observer and gives perspective for how the orbital plane is inclined with respect to Earth.
The relative right ascension and declination values (∆α and ∆δ) with respect to the celestial coordinates of the binary barycenter in the sky (denoted by α 0 and δ 0 ) use combinations of these constants (as described in O'Neil et al. 2019), such that and where d is the distance between the solar system barycenter and the binary barycenter. These sky projection values are in units of radians. While proper motion and parallax do fall under the category of astrometry, since they affect the observed positions of the binary system in the sky, they are calibrated out by subtracting the binary barycenter coordinates from the observed right ascension and declination values, as done in Equations (9k) and (9l). The binary barycenter coordinates, α 0 and δ 0 , change over time due to proper motion. Given some sky coordinate, a i and δ i , at initial time t 0 , after time t obs , the new right ascension becomes and the new declination becomes Parallax and proper motions still contribute to lineof-sight motions, however, and we discuss the effects on observed radial velocities in Section 4.5.
The various Newtonian and relativistic effects enter the astrometric model in two ways, via their contribution to: (i) the relative positions of the stars and the black hole in the frame of the binary barycenter and (ii) the light propagation time between the star and the observer. Indeed, Equation (8a) for the calculation of the relative positions is written in terms of t star , which is determined by propagation effects and is related to the time of observation t obs via Equation (2). We provide explicit equations for the various propagation effects in the following section, since we will use the same equations to derive the various spectroscopic effects.
In the astrometric model, so far, we have neglected the effects of gravitational lensing of the stellar positions in the black-hole spacetime. These effects have been explored, e.g., in (Nusser & Broadhurst 2004) and Bozza & Mancini (2012), and, albeit non-detectable with current data, might be relevant for modeling future observations (Grould et al. 2017).

SPECTROSCOPIC MODEL
As discussed in Section 1, pulsar timing models (e.g., D&D II; Edwards et al. 2006) use the time delays between emission and observation of pulses from binary pulsar systems to fit for the orbital parameters. With the S stars, we have both astrometric as well as spectroscopic data. In the previous section, we presented the equations we use to model the astrometric data. In this section, we present new line-of-sight velocity equations, which we derive from the time delay equations given in Edwards et al. (2006).
Our spectroscopic model incorporates five timing effects (presented in the following subsections): the binary Roemer effect ( §4.1), the periastron precession ( §4.2), the binary Einstein effect ( §4.3), the binary Shapiro effect ( §4.4), and the Earth Roemer effect ( §4.6). We also derive velocity equations for the Kopeikin effect ( §4.5) and the solar Shapiro effect (which we do not report here for brevity) to check their magnitudes and confirm that they are negligible. This allows us to omit other timing delays listed in Edwards et al. (2006) that are caused by Earth, the solar system, interstellar space, as well as higher-order GR effects, such as frame dragging, quadrupole moment, and additional 2PN terms described in Angélil et al. (2010).
All time derivatives listed in this section are with respect to the observer time (t obs ), or the arrival time of the light from the star as seen by the observer.

Binary Roemer Effect
The binary Roemer delay component (i.e., the component of the Roemer delay that pertains only to the movement of the star around the binary barycenter) is denoted by ∆ RB . It is simply the light travel time for the line-of-sight distance across the orbit divided by the speed of light (i.e., ∆ RB = −n · r/c). Adapting notation from both D&D II and Edwards et al. (2006), we calculate the binary Roemer delay as: Taking the time derivative of Equation (10a) gives the equation of radial velocity resulting from the binary Roemer effect where du/dt is derived from the Kepler Equation (Equation 8a) such that The additional term v ω results from the time dependence of the argument of periastron, ω. It is known as the periastron precession, which we introduce below.

Periastron Precession
Periastron precession is a first-order general relativistic effect that causes the orbit of an object to shift or rotate around its pericenter over time. As the location of periapsis advances around the binary barycenter, the slight change in direction from a regular elliptical orbit results in an observed, line-of-sight velocity boost for the orbiting object. The dω/dt term in the time derivative of the binary Roemer delay (Equation 10b) describes this spectroscopic effect, The time derivative of the argument of periastron (dω/dt) relates to the time derivative of mean anomaly du/dt (Equation 10c) by where K is again the dimensionless parameter related to the total momentum of the system and determines the amount of orbital precession (Equations 5c-5d and 8g).

Binary Einstein Effect
The binary Einstein delay (∆ EB ) accounts for the difference between the star time of emission (t star ) and the binary barycenter time (t BB ), which occurs due to both gravitational redshift (from GR) and time dilation (from special relativity). We write the binary Einstein time delay as where m 1 is the companion mass to the orbiting star (since we define m 1 ≥ m 2 ). The derivative of this time delay, then, is simply which, when we substitute in the expression for du/dt (Equation 10c), becomes In the case of a circular orbit, Equation (12c) results in a null effect, even though gravitational and relativistic Doppler effects are still present. The reason is that this equation only describes the change in the combined wavelength shift along an eccentric orbit introduced by the changing separation and velocity magnitude. The baseline effect is absorbed into the systemic radial velocity of the object, i.e., the line-of-sight proper motion (µ ), which is a free parameter, and can be trivially obtained in post-processing.

Binary Shapiro Effect
The binary Shapiro delay (∆ SB ) describes the general relativistic effect introduced when light from one of the bodies in the model passes through the gravitational well of the other. We rewrite Equation (73) from Edwards et al. (2006) in geometrized units as ∆ SB = −2m 1 ln 1 − e cos u − sin i sin ω(cos u − e) where the usage of e without a subscript indicates that any eccentricity (i.e., e R , e θ , e t ) may be used, as the difference is negligible. This results from the fact that the expressions for the three eccentricities differ at the ∼ 1/c 2 order, which we omit in the 1PN limit since they already multiply a term of the same order. The time derivative of the binary Shapiro delay (∆ SB ) gives the binary Shapiro effect for spectroscopic measurements: where we have introduced the symbol A, i.e., A ≡ 1 − e cos u − sin i sin ω (cos u − e) + cos ω sin u 1 − e 2 . (13c)

Kopeikin Effect
The Kopeikin effect combines spectroscopic contributions from both proper motion of the binary system and the parallax as seen from Earth. As discussed in Section 3, this is an astrometric effect, as well, but is calibrated out by subtracting the binary barycenter location from all observed positions of the two bodies for each time step or observation.
We base our definition of the Kopeikin effect (∆ KB ) on Sections 2.7.1 and 2.7.2 in Edwards et al. (2006), where it is broken into three components (Edwards et al. 2006, their Eq. 72): Here, ∆ SR is caused by changes in the viewing angle geometry due to the proper motion, ∆ AOP is due to the annual orbital parallax, and ∆ OP is due to the orbital parallax (the orbital equivalent of the Shklovskii effect). These three components, rewritten in geometrized units from Edwards et al. (2006, their Eqs. 73-75), are ∆ SR = a r sin i (t a,BB − t 0 ) × (µ α * sin Ω + µ δ cos Ω) C csc i + (µ α * cos Ω − µ δ sin Ω) S cot i , (14b) and where C and S are given by (65) and (66) in Edwards et al. (2006) and µ α * = µ α cos δ.
The annual orbital parallax and orbital parallax distances (to the binary barycenter) are approximately equal such that d ≈ d AOP ≈ d OP . As with the previous effects, the total Kopeikin effect is simply the time derivative of the above equations, where: and (14h) Here, the symbols B and C are defined as: The time derivatives of the expressions for C and S are and

Earth Roemer Effect
The Earth Roemer effect is analogous to the binary Roemer effect ( § 4.1), except that the changes in light travel time are due to the motion of the Earth around the solar system barycenter. We define the Earthmotion Roemer delay component using similar notation to Equations (13)-(16) in Edwards et al. (2006): where r ⊕ is the vector from the solar system barycenter to the Earth: andR BB is the unit vector between the binary barycenter and the observer: The binary barycenter-observer unit vector (R BB ) breaks down to the primary component (η), which is the initial unit vector between the observer and the binary barycenter (in celestial coordinates), i.e., and the shift in the binary barycenter position over time due to proper motion (µ α , µ δ , µ ) of the binary system. While the line-of-sight proper motion µ is one of the free parameters for the model, the transverse proper motion µ ⊥ depends on the right ascension and declination proper motions, µ α and µ δ . We define the transverse proper motion as where the projected right ascension and declination vectorsα andδ are given by Edwards et al. (2006, their Eqs. 17-18) The Earth Roemer effect is the time derivative of the Earth Roemer delay (Equation 15a), such that where the symbols D 1 , D 2 , D 3 , and D 4 are shorthand for and

SUMMARY OF MODEL
Having presented all of the equations and parameters we use in our analytic, 1PN model for S stars, we summarize in Figure 5 the steps one should take to implement it.
The Earth Roemer and Kopeikin effects both depend on the position of the Earth around the Sun, and so to calculate them, we must use Earth ephemerides. For the time derivatives of any values that use these data, we use three-point midpoint differentiation.

Magnitude of Spectroscopic Effects
Due to the exquisite accuracy of pulsar observations, pulsar timing models must incorporate numerous time delays in addition to the geometric and relativistic effects from the binary system. Such effects include atmospheric delays, dispersion from both the interplanetary medium and interstellar medium, frequency-dependent delays, and gravitational effects from many bodies in our solar system, namely Venus, Jupiter, Saturn, Uranus, and Neptune (in addition to the Sun).
Timing models for spectroscopic observations of S stars, on the other hand, do not require the same accuracy. The current sensitivity level for Doppler shifts on the instruments capable of observing the S stars (i.e., Keck and the Very Large Telescope) are around 10 km s -1 (Thatte et al. 1998;Martin et al. 2018) and upcoming ground-based 25-40 m extremely large telescopes (ELTs) are anticipated to have velocity sensitivities of around 1 km s -1 (e.g., Mawet et al. 2019;Marconi et al. 2021). As a result, we may neglect timing effects that are significantly smaller than predicted ELT-class sensitivities. Note that the ELT-class sensitivities have been estimated for stars down to ∼ 19 mid-infrared magnitude, but will depend on the actual properties of stars with smaller orbital separations that may be discovered in the future.
In considering the possibility of detecting higher-order general relativistic effects, it is useful to consider the order-of-magnitude strengths of the spectroscopic effects described in Section 4. We derived scaling equations for the six velocity components of the spectroscopic model and present them in Appendix A.
Using these scaling relations, we show in Figure 6 the relative strengths of the binary Roemer effect, binary Einstein effect, periastron precession, binary Shapiro effect, Earth Roemer effect, and Kopeikin effect for orbital periods that range from 0.5 years to 500 years with a fixed system mass and two-body mass ratio that is characteristic of Sgr A* and the S stars. The widths of the bands result from a range of eccentricities, e = 0.7−0.98, with the smaller effects corresponding to lower eccentricities.
Another notable aspect of Figure 6 is the difference in the slopes of the various spectroscopic effects. As a purely geometric phenomenon, wavelength shifts introduced by the binary Roemer effect scale as ∆λ/λ ∼ a −1/2 R ∼ P −1/3 , growing larger with smaller semimajor axes. The magnitudes of the periastron precession and binary Shapiro effects also increase with smaller semimajor axes and periods, although they do so much more rapidly (∆λ/λ ∼ a −3/2 R ∼ P −1 ). This is because of the steeper scaling of the PN corrections with orbital separation. Similarly, the binary Einstein effect lies in between these scalings (∆λ/λ ∼ a −1 R ∼ P −2/3 ). The Kopeikin effect, on the other hand, is not a GR effect. Rather, it describes the radial velocity contributions from the parallax and proper motions of the binary system. As a result, since there is more time for parallactic effects and proper motions to affect observations for any given orbit, the magnitude of the corresponding Doppler effect grows as the orbital period increases (∆λ/λ ∼ a R ∼ P 2/3 ). The magnitude of the Earth Roemer spectroscopic effect is independent of binary orbital period and, as such, has a constant magnitude of ∼ 30 km s -1 .
One last property to note is the effect of eccentricity on the magnitudes of the spectroscopic effects. For any fixed orbital period, the binary Roemer effect has a span of two orders of magnitude, the binary Einstein effect ranges around two orders, the periastron precession stretches almost four orders of magnitude, the binary Shapiro effect extends over three orders, and the Kopeikin effect spans one order. The three GR effects -binary Einstein effect, periastron precession, and binary Shapiro effect -are most dramatically influenced by changes in eccentricity. Higher eccentricities bring orbiting stars increasingly closer to the central black hole, resulting in higher angular momenta of the system (Equation 5b-5c) and larger amounts of orbital precession (∆θ), as well as greater Doppler contributions from the periastron precession effect. Similarly, the smaller Period (yr) Figure 6. Absolute magnitudes of characteristic contributions to the radial-velocity corrections in spectroscopic models of stellar orbits around Sgr A* introduced by the various Newtonian and post-Newtonian effects. Each shaded area corresponds to orbital eccentricities in the range 0.7 -0.98, with lower limits corresponding to e = 0.7 and upper limits corresponding to e = 0.98. The orientation angles of the orbits are assumed to be those of the S0-2 star. Horizontal dashed lines show the current measurement limits for VLT/Keck observations as well as the expected limits for an ELT-class telescope. We see that while the binary Roemer effect (far left, solid red band) is the dominant velocity contribution, the three GR effects -binary Einstein effect, periastron precession, and binary Shapiro effect -rapidly increase at smaller semimajor axes/orbital periods. The periastron precession effect (second from right, purple band) has the strongest dependence on eccentricity and spans four orders of magnitude, overlapping with the binary Einstein effect (second from left, yellow band) and Shapiro effect (far right, light blue band), which span one and three orders of magnitude, respectively. While the binary system spectroscopic effectsbinary Roemer, binary Einstein, periastron precession, and binary Shapiro -decrease for larger orbital separations (semimajor axes), the Kopeikin effect describes velocity contributions due to parallax and proper motions that are due to the movement of the observer and, as such, increases for larger orbital separations.
the periapsis distance is (due to higher eccentricities), the closer the light from orbiting stars must pass through the potential well of the central black hole, which also increases the binary Shapiro effect. We can reach several significant conclusions from Figure 6. First, based on the expected sensitivities of ELT-class telescopes of ∼1 km s -1 , we find that the binary Shapiro effect should be detectable by these telescopes through spectroscopic observations alone during the next S0-2 periapsis passage in 2034. Second, if observations detect and confirm S stars with shorter periods ( 10 yr) and high eccentricities ( 0.85 for periods of less than a year and 0.95 for periods of less than 10 yr), existing telescopes with velocity sensitivities of ∼10 km s -1 may already be capable of detecting the first-ever binary Shapiro effect for any S star.

Fingerprints of Spectroscopic Effects
One of the key advantages of the spectroscopic model we presented in Section 4 is that each effect has a unique signature (or "fingerprint") that we can search for in the data. Figure 7 shows the four binary system spectroscopic effects for the S0-2 star and a hypothetical star with period less than ten years, which we refer to as S0-X. For the sake of illustration, we picked a set of fiducial parameters similar to those listed in Peißker et al. (2020a), although such a star has yet to be confirmed (Gravity Collaboration et al. 2022b). Table 3 in Appendix C lists the parameter values used for the S0-2 and S0-X models.
Clearly, the patterns of the three effects vary greatly, both with respect to each other and between stars. One of the primary factors driving this is the orientation of the stellar orbit with respect to the observer. Both the binary Roemer effect and the periastron precession relate to the geometry of the orbit (Newtonian and GR): the binary Roemer effect is the change in light-travel time and the precession of the periapsis point due to GR results in an additional line-of-sight velocity boost as the orbit precesses. As a result, both introduce radial velocity contributions that have the same signs. For S0-2, the effects are primarily positive, while for S0-X, they are primarily negative.
The binary Shapiro effect, on the other hand, is stronger when emitted light from an orbiting star passes closer to the black hole en route to the observer. This is opposite of the effect that governs the sign of the binary Roemer and periastron precession effects, and so it primarily has the opposite sign of the other two. Depending on the orientation of the orbit with respect to the Earth, the binary Shapiro effect rapidly changes sign as the star moves behind the black hole.
While the three aforementioned spectroscopic effects all exhibit sign changes, the binary Einstein effect is always positive. This is because, by the nature of the gravitational redshift, the emitted light can only be redshifted.
The net effect of the different behaviors and signs of the spectroscopic effects is that they act like fingerprints on radial velocity observations over time. The spectroscopic signature of the periapsis precession (panels second from top) cannot be confused with that of the binary Shapiro effect (bottom panels), as they exhibit very different functional forms and signs. This is valuable in the analysis of S-star orbits, as it enables the identification of both binary system geometric and GR effects and solar system geometric and GR effects and also minimizes the chance of confusing them. Figure 8 zooms in on the spectroscopic effects shown in Figure 7 to compare the timings of the maxima and minima of each of the radial velocity contribution effects.
The velocities are scaled to similar orders of magnitude to elucidate these comparisons. We can see that all three effects are offset from each other. The binary Roemer effect peaks in magnitude first, the binary Einstein, second, the periastron precession, third, and the binary Shapiro, fourth. This is crucial, as it demonstrates that the velocity contributions will not cancel out and, therefore, fitting spectroscopic effects via the residual method is a viable way to detect them. Heißel et al. (2022) have come to a similar conclusion for different effects in the astrometric domain.

CONCLUSIONS
The results from our astrometric and spectroscopic model have several key implications. First, calculations using the orbital parameters of the S0-2 star (Gravity Collaboration et al. 2020) show that the binary Shapiro effect should be detectable spectroscopically by ELTclass telescopes during its next periapsis period. For any stars discovered with shorter periods and/or higher eccentricities, these effects could already potentially be detected with current instruments. Although stellar winds result in rotational broadening of lines, winds from B stars like S0-2 (and other stars in the Galactic center) are weak and not directly measurable (Martins et al. 2008;Fang & Chen 2021). Furthermore, since observational studies of the S stars use broad absorption lines to determine line-of-sight velocities, the systematic uncertainties of ∼ 10 km s -1 mean that detecting the described spectroscopic effects are still realistic. Overall, our model shows that while astrometric data are crucial for constraining the orbital parameters (i.e., period, orientation, eccentricity) of a particular star, they are not essential for detecting general relativistic effects, which can be done entirely via spectroscopy.
Second, the analytic nature of the model allows us to evaluate directly the various observables at any time, without having to integrate the differential equations of motion for each star. This results in a remarkable reduction of computational cost, as it can be easily demonstrated with fiducial values for the hypothetical S0-X star. For example, due to the S0-X orbital parameters and orientation, the velocity correction from the periapsis precession changes dramatically over the course of one-hundredth of its orbit by 1000 km s -1 , i.e., by 100 times the magnitude of the current measurement uncertainties. Integrating the geodesic equations with a method of order n (e.g., such that n = 4 for a fourth order Runge-Kutta method) using a timestep ∆t would introduce a fractional truncation error in each timestep that is ∼ (∆t/10 −2 P ) n+1 , where we have used the fact that there is significant evolution over a frac- Figure 7. The contributions of different Newtonian and post-Newtonian effects to the line-of-sight velocities of (left) the S0-2 star and (right) S0-X, a hypothetical, sub-10 yr period star, around Sgr A*. Unique functional forms for the different spectroscopic effects make it easy to disentangle geometric, Earth-related, and solar-system-related effects from true GR phenomena. The eccentricity and orientation of its orbit make post-Newtonian effects in a star like S0-X potentially detectable even with current instruments. S0-X Components v RB × 10 −4 v EB × 10 −3 v ω × 10 −3 v SB × 10 −2 Figure 8. Velocity contributions of the binary Roemer effect, binary Einstein effect, periastron precession, and binary Shapiro effect for the periapsis passages of S0-2 (left) and S0-X, a hypothetical, sub-10 year period star (right), around Sgr A*. Velocities are scaled to better show the unique shapes of the velocity components over time. Timescale bars in the lower left of both plots show a one-month period and one-day period for S0-2 and S0-X periapsis passages, respectively, as an indicator of useful observational cadences. The peaks and troughs of the four spectroscopic effects are offset from each other in time. This implies that they will not cancel out, which means the residual-fitting method is a viable way to attempt to detect GR phenomena with this model.
tion 10 −2 of the orbital period P . The total accumulated error after integrating for a single orbit will be ∼ (P/∆t)(∆t/10 −2 P ) n+1 . Requiring for this error to be smaller than the measurement error, i.e., to be ∼ 10 −3 , leads to the conclusion that we would need at least P/∆t ∼ 100 2+5/n timesteps per orbit, which is 10 4 for n ≥ 2. With our analytic model, evaluating the precise correction at any point in time will take only a handful of floating point operations. This factor of ∼ 10 4 reduction in computational cost of using an analytic model compared to a numerical one implies that a Bayesian MCMC statistical study that takes an hour for the analytical model will need about a year for the numerical one. This allows more efficient parameter space exploration as well as the possibility of fitting multiple stellar orbits with simultaneous constraints. Upcoming ELT-class telescopes will likely be pivotal in producing more data for stars in the Galactic center over shorter integration times.
While the post-Newtonian approach has many computational and scientific advantages, it does not allow for a straightforward incorporation of the perturbing ef-fects of an extended mass distribution (see, e.g., Jiang & Lin 1985;Rubilar & Eckart 2001). However, the current limits on the mass of an unseen perturber are of the order of ∼3000 M (Gravity Collaboration et al. 2020Collaboration et al. , 2022aHeißel et al. 2022) for extended mass within the S2 orbit and ∼1000 M for a compact mass within the inner arcsecond of Sgr A* and within or around the orbit of S2 Gravity Collaboration et al. 2020). Given the current accuracy of the astrometric and spectroscopic data, the presence of such a perturbing mass is not significant and does not compete with the effects that we have considered at the 1PN order.
We thank Norbert Wex and Gunther Witzel for carefully reading the manuscript and for their comments and suggestions. We also thank the members of the Extreme Astrophysics group at the University of Arizona for useful comments and suggestions throughout this project. This work was supported in part by NSF PIRE award OISE-1743747 and NSF award AST-1715061. To estimate the relative strengths of the various velocity effects, we simplified the equations in Section 4 to scaling relations, which are given below. These relations were used to calculate the values in Figure 6.