Relativistic Effects in Orbital Motion of the S-Stars at the Galactic Center

Abstract

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.

1. 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). Therefore, we are able to explore them in the frame of Post-Newtonian formalism. For recent overviews of General Relativity (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 Parameterized 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 $\beta$ and $\gamma$, 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 $\gamma$ parameter (with the accuracy of $2\times {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 $\gamma$ and $\beta$ were also measured with the EPM2017 ephemerides [11,12]. It turns out that $\beta$ is at the ${10}^{-5}$ level as well, and that also the $\gamma$ measurement is reaching the Cassini experiment level. In particular, [13] delivers ${\mathrm{err}}_{\beta }=3.9\times {10}^{-5}$ and [14] yields ${\mathrm{err}}_{\beta }=1.8\times {10}^{-5}$. Iorio [11], relying upon recent 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].

Post-Newtonian equations of motion also can be used for a possible solution of the well-known problem of pseudo-tensor of the energy–momentum tensor of the gravitational field in GR (Landau and Lifshitz [18] (Par. 101 “The energy-momentum pseudotensor”, p. 304) and review Baryshev [19]). The point is that in the frame of the Feynman’s relativistic quantum-field approach to gravitational interaction [20] there is ordinary (not “pseudo”) energy–momentum tensor $T^{\mu \nu }$ of the gravity field. So there is positive localizable energy density of the gravitational field $T_g^{00}$, which in the case of the PN approximation equals ${(\nabla {\phi }_N)}^2/8\pi G(erg/{\mathrm{cm}}^3)$. Fortunately, this energy–mass distribution around central massive body gives 16.7% contribution to the relativistic pericenter shift and together with linear part gives the standard value of the Mercury perihelion shift [21].

The S-stars were being used in a big amount of different investigations. The most interesting in the context of this work are investigations of the motion of the S-stars. The possibility of detecting the Post-Newtonian effects in the motion of S2 was being researched some time ago [22]. It was shown that the perturbations produced by these effects are significant and have a specific time dependence. In the work Parsa et al. [23], the authors used the PN approximation to calculate the stellar orbits and investigate how do they change if the fitting is employed in different parts of the orbit. The important effects that can be seen in the S-stars data are the effects related to the light propagation. It is possible to investigate the relativistic redshift using the spectroscopic data. The authors of work Saida et al. [24] show that there is a double-peak-appearance in the time evolution of the defined GR evidence.

The purpose of this work is to consider the relativistic Post-Newtonian motion of several S-stars, using the data from the work [25], and to attract attention to the theoretical uncertainty in analysis of PN motion of S-stars due to coordinate dependence of the orbital shapes of the PN test body motion in GR. Additionally in the frame of the Feynman’s field-theoretical approach (FGT) we estimate the energy density of the static gravitational field for S-stars conditions.

The effects related to spin of the central massive object, such as Lense–Thirring precession, are also detectable via observations of the most recently discovered S-stars [5,26,27]. It is important that we consider only the Schwarzschild problem and do not take the spin-related effects into consideration, which we shall study in separate paper.

2. Observational Data

The S-stars cluster has been 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 are 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. What 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 consist of radial velocities, corrected with respect to the Local Standard of the Rest (the so-called VLSR correction) of the S-stars. We cannot use the modeled velocities directly to fit the observed radial velocities. We shall interpret the latter as the redshift values. The modeled 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].

3. Post-Newtonian Effects

3.1. 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 [32] (1991), the Post-Newtonian equations of motion of a test particle in the static gravitational field of central massive body are given by

$$\ddot{\mathbf{x}}=-\nabla {\phi }_N\left(1+(4-2\alpha )\frac{{\phi }_N}{c^2}+(1+\alpha )\frac{{\dot{\mathbf{x}}}^2}{c^2}-3\alpha \frac{{(\mathbf{x}·\dot{\mathbf{x}})}^2}{c^2{x}^2}\right)+(4-2\alpha )\left(\nabla {\phi }_N·\frac{\dot{\mathbf{x}}}{c}\right)\frac{\dot{\mathbf{x}}}{c},$$

(1)

where $\alpha$ denotes the coordinate system choice. For example, $\alpha =1$ for the Standard (Schwarzschild) coordinates and for either isotropic or harmonic coordinates the parameter $\alpha =0$. Note that, though the secular effects do not depend on the parameter $\alpha$, the shape of the orbit is different for different $\alpha$. 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, for the first time, to measure directly the shape of the S-star orbit, which means to test the value of the parameter $\alpha$ itself. This raises a question about the role of the coordinate system choice in GR.

For isotropic or harmonic coordinates, the parameter $\alpha =0$, so the corresponding PN equations of motion in the frame of GR are

$$\ddot{\mathbf{x}}=-\nabla {\phi }_{\mathrm{N}}\left(1+4\frac{{\phi }_{\mathrm{N}}}{c^2}+\frac{{\dot{\mathbf{x}}}^2}{c^2}\right)+4\left(\nabla {\phi }_{\mathrm{N}}·\frac{\dot{\mathbf{x}}}{c}\right)\frac{\dot{\mathbf{x}}}{c}.$$

(2)

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\frac{{\phi }_{\mathrm{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

$$T_g^{00}={ϵ}_g={(\nabla {\phi }_N)}^2/8\pi G(erg/c{m}^3).$$

(3)

This additional mass distribution gives correction to the 00-component of the tensor gravitational potential ${\psi }^{\mu \nu }$ in the form:

$${\psi }^{00}=\varphi ={\phi }_N+\frac{1}{2}\frac{{\phi }_{\mathrm{N}}^2}{c^2},$$

(4)

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 by

$$\ddot{\mathbf{x}}=-\nabla {\phi }_{\mathrm{N}}\left(1+2(\beta +\gamma )\frac{{\phi }_{\mathrm{N}}}{c^2}+\gamma \frac{{\dot{\mathbf{x}}}^2}{c^2}\right)+2(\gamma +1)\left(\nabla {\phi }_{\mathrm{N}}·\frac{\dot{\mathbf{x}}}{c}\right)\frac{\dot{\mathbf{x}}}{c}.$$

(5)

Substituting $\beta =\gamma =1$ leads to GR PN equations in the frame of the isotropic or harmonic coordinate systems. We can see that in the limit $c\to \infty$ these equations reduce to Newtonian equations of motion

$$\ddot{\mathbf{x}}=-\nabla {\phi }_{\mathrm{N}}.$$

(6)

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. Therefore, Brumberg [32] used to have 3 PPN parameters $(A,B,K)$, instead of the 2 usual PPN parameters $(\beta ,\gamma )$.

3.2. 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]

$$\Delta \omega =\frac{6\pi r_{\mathrm{g}}}{a(1-e^2)},$$

(7)

where a is a semi-major axis, e is an eccentricity, and $r_{\mathrm{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:

$$\Delta \omega =\frac{6\pi r_g}{a(1-e^2)}=\frac{7\pi r_g}{a(1-e^2)}-\frac{\pi r_g}{a(1-e^2)}$$

(8)

The first term with a factor of $7\pi$ 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 $\pi$ 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 $\frac{{\phi }_N}{c^2}$ in pericenter predicted for S-stars are presented in the

Table 1. lS-stars pericenter shift.

Star	S2	S38	S55
$\Delta \omega$	${12}^{\prime }$	${7.1}^{\prime }$	${6.7}^{\prime }$
$\dot{\omega }$	${45}^{″}/\mathrm{yr}$	${22}^{″}/\mathrm{yr}$	${31}^{″}/\mathrm{yr}$
$\dot{\omega }·100\mathrm{yrs}$	$1^{\circ }{15}^{\prime }$	${37}^{\prime }$	${52}^{\prime }$
${\left(\frac{{\phi }_N}{c^2}\right)}_{\mathrm{per}}$	$-3.48\times {10}^{-4}$	$-1.99\times {10}^{-4}$	$-1.69\times {10}^{-4}$

.

The observable effect is similar to one with the 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 $\frac{{\phi }_N}{c^2}$ is $-3.21\times {10}^{-8}$, and for binary pulsar PSR 1913 + 16 it is $-2.7\times {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.

3.3. 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 modeled RVs to the observable RVs. 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

$${\mathrm{RV}}_{\mathrm{obs}}/c={\mathrm{RV}}_{\mathrm{model}}/c+\frac{{\phi }_{\mathrm{N}}}{c^2}+\frac{v^2}{2c^2}.$$

(9)

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.

4. Modeling

4.1. Parameters

Our aim was to build a model of the S-stars observations which depends on PPN parameters $\beta$ and $\gamma$. After that, we will be able to use some of the optimization techniques to solve the inverse problem and estimate the parameters. 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 $\mathfrak{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]. The final 2 of our parameters are the PPN parameters of $\beta$ and $\gamma$ 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. As a result of this, along with the large amount of parameters, we decided to use the Bayesian inference techniques, such as the MCMC sampling.

4.2. Constructing a Model

The first step in the modeling 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,{\dot{x}}_0,{\dot{y}}_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,\dot{x},\dot{y})$, which is a modeled trajectory. To rotate it, we need 2 angles: the longitude of the ascending node $\mathsf{\Omega }$ and the inclination i, which are the remaining 2 parameters. To perform the rotations, one should pad the modeled array of phase vectors $(x,y,\dot{x},\dot{y})$ to 6D phase space $(x,y,z,\dot{x},\dot{y},\dot{z})$, where z and $\dot{z}$ are zeros. After that, one should multiply each triplet of coordinates $(x,y,z)$ and velocities $(\dot{x},\dot{y},\dot{z})$ with the rotation matrices that contain $\mathsf{\Omega }$ and i. As the result, one will obtain the modeled 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 modeled 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 the RA scale is clockwise. We also must take into account the proper motion of Sgr A*.

Finally, we have to obtain RV from $\dot{z}$, using the Equation (9) due to the light propagation effects.

5. MCMC Analysis

Bayesian inference is the statistical method of inference in which the probability of a hypothesis is updated with the Bayes’s theorem

$$P\left(H\right|E)=\frac{P\left(E\right|H\left)P\right(H)}{P\left(E\right)},$$

(10)

where H stands for the hypothesis and E stands for evidence. The $P\left(E\right|H)$ is a likelihood. It is a probability of observing the data (evidence) given hypothesis. $P\left(E\right)=P\left(E\right|H\left)P\right(H)+P(E|\neg H)P(\neg H)$ is called model evidence or the marginalized likelihood. $P\left(H\right)$ stands for the prior probability (the probability of hypothesis before taking the observed data into the account) and $P\left(H\right|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

$$P\left(\mathbf{\theta }\right|\mathbf{X})=\frac{P\left(\mathbf{X}\right|\mathbf{\theta }\left)P\right(\mathbf{\theta })}{P\left(\mathbf{X}\right)}\propto P\left(\mathbf{X}\right|\mathbf{\theta })P\left(\mathbf{\theta }\right),$$

(11)

where $\mathbf{X}$ is a sample of the observed data and $\mathbf{\theta }$ is a vector of the parameters. $P\left(\mathbf{\theta }\right)$ is a prior distribution of the parameters, $P\left(\mathbf{X}\right|\mathbf{\theta })$ is the likelihood, and $P\left(\mathbf{\theta }\right|\mathbf{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-squared 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 a lesser 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 2000 iterations.

The posterior distributions of the ${\beta }_{\mathrm{PPN}}$ and ${\gamma }_{\mathrm{PPN}}$ are presented in Figure 1. The right and top plots show us the histograms of the posterior distributions of the considered parameters. The bottom-left plot is a ‘2D-histogram’ of these parameters, where the color shows us the amount of points in a bin. This ‘2D-histogram’ is actually a projection of our 31-dimensional posterior distribution of the parameters to the 2D plane of ${\beta }_{\mathrm{PPN}}$ and ${\gamma }_{\mathrm{PPN}}$ parameters. We can see that the PPN parameters do not make a significant contribution to the model. Therefore, these parameters have relatively high errors. It is possible to constrain the PPN parameters with this approach. The errors are much smaller than in the previous work [25].

The other important parameters constrained by our model are the mass of Sgr A* in the units of ${10}^6{M}_{\odot }$ and the distance to the Galactic Center in kpc. The posterior distributions of these parameters are presented in Figure 2. In this posterior distribution, it is clearly visible that there exists a correlation between these parameters, as the posterior distribution points line up. Orbital parameters of S2, S38, S55, and reference frame offset values are given in Appendix A.

Figure 3, Figure 4, Figure 5 and Figure 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.

6. 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*.

$$d=2.44R_{*}{\left(\frac{M_{\mathrm{SgrA}*}}{M_{*}}\right)}^{\frac{1}{3}}\approx 2.44R_{\odot }{\left(\frac{M_{\mathrm{SgrA}*}}{M_{\odot }}\right)}^{\frac{1}{3}}=2\mathrm{AU},$$

(12)

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. After the pericenter passage, the G2 object conserved without tidal disruption, and it became clear that it is a star.

7. Discussion and Conclusions

In this work the PPN parameters $\beta$ and $\gamma$ 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 $\mathfrak{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 estimates of PPN parameters in the vicinity of the supermassive relativistic compact object Sgr A* are $\beta ={1.01}_{-0.21}^{+0.29}$ and $\gamma ={0.99}_{-0.29}^{+0.19}$. The PPN scaling parameter of $f_{\mathrm{SP}}=(2+2\gamma -\beta )/3$ is estimated to be ${1.01}_{-0.17}^{+0.20}$, which is comparable with the values obtained by GRAVITY collaboration [31]. The result is consistent with the General Relativity predictions. The formal estimation of mass of Sgr A* and the distance to the Galactic Center are $\mathfrak{M}={4.50}_{-0.03}^{+0.04}·{10}^6{M}_{\odot }$ and $R_0={8.53}_{-0.03}^{+0.03}\mathrm{kpc}$. We note that the question about possible contribution of the choice of coordinate system in PN equations of motion needs further investigations.

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 ${(\nabla {\phi }_N)}^2/8\pi G$. Future studies of the motion of the closest to SgrA* S-stars will test the gravity theories and give the crucial information on the nature of the central relativistic compact object SgrA*.