Relativistic effects in orbital motion of the S-stars at the Galactic Center

The Galactic Center star cluster, known as S-stars, is a perfect source of relativistic phenomena observations. The stars are located in the strong field of relativistic compact object Sgr A* and are moving with very high velocities at pericenters of their orbits. In this work we consider motion of several S-stars by using the Parameterized Post-Newtonian (PPN) formalism of General Relativity (GR) and Post-Newtonian (PN) equations of motion of the Feynman's quantum-field gravity theory, where the positive energy density of the gravity field can be measured via the relativistic pericenter shift. The PPN parameters $\beta$ and $\gamma$ are constrained using the S-stars data. The positive value of the $T_g^{00}$ component of the gravity energy-momentum tensor is confirmed for condition of S-stars motion.


Introduction
The S-stars cluster in the Galactic Center [1][2][3][4] is a unique observable object to investigate. Once for the orbital period the stars pass their pericenters, where they reach high velocities (∼ 0.01 or even ∼0.1 of the speed of light for the most recent discovered stars [5]) and get close to the relativistic compact object Sgr A*. Such conditions lead to manifestation of relativistic effects in proximity to a supermassive compact object (Super-Massive Black Hole). So we are able to explore them in the frame of Post-Newtonian formalism. For recent overviews of GR see, e.g., Debono and Smoot [6] and references therein.
The Post-Newtonian formalism is a tool which provides an opportunity to express the relativistic equations of gravity as the small-order deviations from the Newtonian theory. The Parametrised Post-Newtonian formalism can be implemented to various theories of gravity [7]. It is a version of the PN formalism that has a parameters that quantitatively express the difference of a theory of gravity and the Newtonian theory. In this work, the parameters of β and γ, which appear in the expansion of the Schwarzschild solution [8], are taken into consideration. Even in the Solar system, physical effects can be detectable via the PPN formalism. The most precise constrain of the γ parameter (with the accuracy of 2 × 10 −5 ) was obtained using the data from the Cassini experiment [9]. The detected effect was the the time delay, which is also known as the classical fourth test of general relativity by Shapiro [10]. The PPN parameters γ and β were also measured with the EPM2017 ephemerides [11,12]. It turns out that β is at the 10 −5 level as well, and that also γ measurement is reaching the Cassini experiment level. In particular, [13] delivers err β = 3.9 × 10 −5 and [14] yields err β = 1.8 × 10 −5 . Iorio [11], relying uponrecent EPM data, suggests a 10 −6 level. In the modern Solar System experiments, even the second-order PN effects are being investigated [15,16]. The PPN formalism is a promising approach in testing cosmological models [17].

Observational data
The S-stars cluster is being observed for 25 years now. This period of time is comparable with the orbital periods of the stars. Some of them have made a revolution of their orbit, and in particular cases it is even not a single time. The accumulated arrays of data differ in their meaningfullness for different S-stars. The star named S2 is the most popular choice among the researchers, and it is not unreasonable. The S2 is the most thoroughly observed, so it has the biggest time series of the observational data. We will consider this star, along with the stars S38 and S55.
There are two types of S-stars observational data: the astrometric and the spectroscopic one. The astrometric data is given as the visual positions of the S-stars on the sky. It is important that one should take into account not only the pair of coordinates, but also the time of the observations. At first sight it may seem to the reader that to obtain the best-fit trajectory, one should implement the fitting procedures that use the points of visual positions of the stars to fit the ellipses into them as the visual trajectories. But what do these procedures actually fit are not 2D ellipses, but spirals in the 3D space with the axes of not only right ascension and declination, but also the time of the observed positions. The spectroscopic data consists of radial velocities, corrected with respect to the Local Standard of the Rest (the so-called VLSR-correction) of the S-stars. But we cannot use the modelled velocities directly to fit the observed radial velocities. We shall interpret the latter as the redshift values. The modelled radial velocities must be transformed to observable values with the equation (9) so the relativistic Doppler effect and the gravitational redshift are taken into account.
The dataset that we use is described in previous work [25]. It is a dataset combined from the data of 4 works Gillessen et al. [4], Boehle et al. [28], Chu et al. [29], and Do et al. [30]. It contains the astrometric and the spectroscopic data of S2, S38 and S55 obtained on VLT, Keck observatory and Subaru telescope. The data presented in different works must be corrected for difference in the reference system [4,31].

Equations of motion
For the purpose of this work it is important to use the equations of motion, which will be integrated to simulate the S-stars trajectories. In the frame of geometrical General Relativity theory (GR), according to classical textbook by Brumberg (1991) , the Post-Newtonian equations of motion of a test particle in the static gravitational field of central massive body are given bÿ where α denotes the coordinate system choice. For example, α = 1 for the Standard (Schwarzschild) coordinates and for either isotropic or harmonic coordinates the parameter α = 0. Note that, though the secular effects do not depend on the parameter α, the shape of the orbit is different for different α. For the purposes of the relativistic celestial mechanics, the isotropic coordinate system is more commonly used, as it is stated in Misner et al. [8]. Here it is important to note that modern S-stars observations allow, at first time, to measure directly the shape of the S-star orbit, which means that to test the value of the the parameter α itself. This rises a question about the role of the coordinate system choice in GR. For isotropic or harmonic coordinates the parameter α = 0, so the corresponding PN equations of motion in the frame of GR arë Intriguingly, these equations of motion coincide with the PN equations of motion in the Feynman's field-theoretical approach to gravitation (FGT) ( [21,33]), where they do not depend on the choice of coordinate system. Importantly, the term 4 ϕ N c 2 includes the non-linear part: factor (4 = 3 + 1), where 3 is the linear part and 1 is the non-linear part. This non-linear part corresponds to the contribution of the effective radial mass-energy density distribution of the gravitational field itself around the central body This additional mass distribution gives correction to the 00-component of the tensor gravitational potential ψ µν in the form: which leads to the non-linear part of the relativistic pericenter shift. Derivations of all needed equations and possible applications of the PN equations of motion are presented in review [19].
The PPN generalization of the Schwarzschild metric in the frame of isotropic coordinates leads to the PPN equations of motion [25], which are given bÿ Substituting β = γ = 1 leads to GR PN equations in the frame of the isotropic or harmonic coordinate systems. We can see that in the limit c → ∞ these equations reduce to Newtonian equations of motion In the work of Brumberg [32], the author considered the problem in a more general way. He did also obtain PPN equations of motion, but the coordinate system that he considered was not fixed. This fact implies that there appears the corresponding third degree of freedom, which is related to the coordinate system choice. So Brumberg used to have 3 PPN parameters (A, B, K), instead of the 2 usual PPN parameters (β, γ).

Pericenter shift
The Post-Newtonian equations of motion lead us into the effect of the pericenter shift. The value of the pericenter shift per one turn in GR is well known [32], [34] where a is a semi-major axis, e is an eccentricity and r g = GM/c 2 is a gravitational radius. From coincidence of the PN equations of a test particle motion in GR and FGT it follows that the expression for the pericenter shift in FGT and GR is exactly the same. The difference from GR is that the FGT formula has two different terms: The first term with a factor of 7π corresponds to the linear approximation, when in the field equations one does not take into account the non-linear contribution (energy of gravitational field itself). The second term with a factor of π occurs after taking into account the non-linearity due to the positive energy density of the gravitational field. This term corresponds to the measurement of the field energy of the gravitation via pericenter shift observations. The pericenter shift values and parameter of The observable effect is similar to one with Mercury anomalous pericenter shift, that Einstein explained. The effect of GR Schwarzschild perihelion precession was also measured with artificial Earth satellites [35,36], so this effect is well studied in the Solar System. Now we can observe it with S-star cluster in the center of our Galaxy. For Mercury, the pericentral ϕ N c 2 is −3.21 · 10 −8 , and for binary pulsar PSR 1913+16 it is −2.7 · 10 −6 . We can see that the field at the pericenter of the S-stars orbits is stronger by 2 orders of magnitude.
Although the values of the pericenter shift seem to be small, they are actually detectable, which is shown in the recent work of GRAVITY collaboration [31]. The pericenter shift appears as a consequence of the PN equations of motion. That means that using the PN equations of motion to fit the trajectories implies detecting the pericenter shift effect. In the frame of field approach it means that we are able to measure the field energy density, which is positive and localizable.

Light Propagation
We must take into consideration several relativistic effects related to light propagation. Such as the relativistic Doppler effect, which is the reason we must transform modelled RV's to the observable RV's. Another significant light propagation effect is the gravitational redshift, which was successfully measured by Do et al. [30]. Taking both of these effects into account is described in [25] by the equation The detailed study of this effect is described in the work of Iorio [37]. It is also necessary to take into account the Rømer delay. The gravitational lensing and Shapiro time delay effects are negligible in this problem.

Parameters
Our aim is to build a model of the S-stars observations which depends on PPN parameters β and γ. After that, we will be able to use some of the optimization techniques to solve the inverse problem and estimate the parameters. So the question is how many parameters do we need to build a model.
To define the orbit and the position of the star, we need 6 parameters. For example, these could be Keplerian parameters or the components of the initial phase vector. The coordinate system that we use is defined as follows: the xy plane is the sky itself, where y axis points in the direction of declination, x axis has the reversed direction of right ascension, and the z axis is pointed from Sgr A* to the observer. The next 2 parameters that we need are the mass M of Sgr A* and the distance R 0 to the Galactic Center. The Sgr A* itself also has the velocity relative to our Solar System, so we have two more parameters of its initial position on the sky (which alongside with the distance R 0 form the vector of the Sgr A* position relative to the Solar System.) and 3 parameters, which define its velocity. The dataset that we use contains 3 different observational data arrays, that were obtained by different teams, and hence have their own different reference frames. The biggest observational dataset is the one from Gillessen et al. [4]. We will use it as the reference observational dataset aligned with our coordinate system. As we have two more datasets (from the papers [28,30]), we must add 4 more parameters, which are the RA and Dec offsets from the dataset of [28] or [30] and the reference dataset of [4]. And the final 2 of our parameters are the PPN parameters of β and γ that we want to constrain. In total we have 19 parameters if we consider 1 star, or 31 parameters if we consider 3 stars. It is a large amount of parameters. In the next section it is shown that the model construction implies very non-trivial dependency of the likelihood function from the parameters. Because of this along with the large amount of parameters, we decided to use the bayesian inference techiques, such as the MCMC sampling.

Constructing a model
The first step in the modelling is to find the trajectory. We want to integrate the equations of motion in a 2D plane. For this purpose, we need 4 parameters of initial phase vector components (x 0 , y 0 ,ẋ 0 ,ẏ 0 ). They are used as the initial conditions for the integrator, which is the Owren-Zennano 4/3-order method implemented in the DifferentialEquations package [38] of the Julia language [39]. It produces the array of phase vectors (x, y,ẋ,ẏ), which is a modelled trajectory. To rotate it, we need 2 angles: the longitude of the ascending node Ω and the inclination i, which are the remaining 2 parameters. To perform the rotations, one should pad the modelled array of phase vectors (x, y,ẋ,ẏ) to 6D phase space (x, y, z,ẋ,ẏ,ż), where z andż are zeros. After that one should multiply each triplet of coordinates (x, y, z) and velocities (ẋ,ẏ,ż) with the rotation matrices that contain Ω and i. As the result, one will obtain the modelled array of phase vectors in the reference frame that is needed. In total we have 6 parameters that define the orbit and the position of the star, as it was stated before.
The next step is getting RA and Dec from x and y of modelled trajectory. The Rømer time delay is taken into account at this step. We shall divide the coordinates by R 0 , but one important thing is the sign. We should inverse x, as RA scale is clockwise. And we also must take into account the proper motion of Sgr A*.
And finally, we have to obtain RV fromż, using the equation (9) due to the light propagation effects.

MCMC analysis
Bayesian inference is the statistical method of inference in which the probability of a hypothesis is updated with the Bayes's theorem where H stands for the hypothesis and E stands for evidence. The P(E|H) is a likelihood. It is a probability of observing the data (evidence) given hypothesis. P(E) = P(E|H)P(H) + P(E|¬H)P(¬H) is called model evidence or the marginalised likelihood. P(H) stands for the prior probability (the probability of hypothesis before taking the observed data into the account) and P(H|E) is the posterior probability (after taking the data into account). In the context of the bayesian parameter inference, the Bayes's theorem may be interpreted as where X is a sample of the observed data and θ is a vector of the parameters. P(θ) is a prior distribution of the parameters, P(X|θ) is the likelihood and P(θ|X) is a posterior distribution of the parameters. Once the model is defined, it is possible to simulate the data values and compare them with the real observable values.
To do that, one can define the chi-square residuals and the likelihood function. The latter can be used by the Markov Chain Monte Carlo sampler. We use the one implemented in the AffineInvariantMCMC package (documentation available at https://github.com/madsjulia/AffineInvariantMCMC.jl) for The Julia language [39]. Using the MCMC sampling, one can obtain the constraints of parameters of a defined model.
In the previous work [25], some of the parameters were considered to be fixed, so the model was simplified and had the less amount of the parameters. However, in the work [31] it is stated, that this simplification leads to the uncertainties. In our case, this is taken into account and the model now depends on 31 parameters. The process was launched with 300 independently computed chains (the so-called 'walkers'), where each chain consists of the parameters and corresponding log-likelihood values for 2,000 iterations.   The Figures 3, 4, 5 and 6 show us the best-fit orbit trajectories and RV plots. We can see that the trajectories are non-closed. That is not caused by the pericenter shift, but by the proper motion (transverse velocity) of the Galactic Center relative to Solar System.

Tidal disruption
Such a strong gravitational field also produces significant tidal forces. The question is how strong are these forces in relation to the objects orbiting the Sgr A*. We can calculate the Roche limit for the fluid satellite (the distance of its tidal disruption) near the Sgr A*.
where R * and M * are radius and mass of the star. We assumed them equal to the mass and radius of the Sun just for the estimation. For S2, S38 and S55 the pericenter is ∼ 200 AU. The tidal disruption is extremely unlikely. It is possible only for extended objects, which are weakly gravitationally bound. There was a suggestion that the observed object G2 [40] is a gas cloud. But after the pericenter passage the G2 object conserved without tidal disruption, and it became clear that it is a star.

Discussion & Conclusion
In this work the PPN parameters β and γ are constrained using the S-stars data. The difference from the previous work [25] is that the helpful remarks from GRAVITY collaboration [31] were taken into account. This time the parameters of M, R 0 , etc. were not fixed, so they were constrained along with the other parameters. The dataset reference systems difference was also taken into account. The offsets between the reference frames were considered as the varying parameters for the MCMC sampling.
The In the frame of the field-theoretical description of the gravitational interaction the observed pericenter shifts of the S-stars confirm the positive localizable energy density of the gravity field having value (∇ϕ N ) 2 /8πG. Future studies of the motion of the most close to SgrA* S-stars will test the gravity theories and give the crucial information on the nature of the central relativistic compact object SgrA*. We thank the anonymous referee for their helpful suggestions. We also thank the GRAVITY collaboration for their useful comments on S2 data.

Abbreviations
The following abbreviations are used in this manuscript:      6.4 6.5 6.6 6.7 6.8