Constraints on graviton mass from Schwarzschild precession in the orbits of S-stars around the Galactic Center

In this paper we use a modification of the Newtonian gravitational potential with a non-linear Yukawa-like correction, as it was proposed by C. Will earlier to obtain new bounds on graviton mass from the observed orbits of S-stars around the Galactic Center (GC). This phenomenological potential differs from the gravitational potential obtained in the weak field limit of Yukawa gravity, which we used in our previous studies. We also assumed that the orbital precession of S-stars is close to the prediction of General Relativity (GR) for Schwarzschild precession, but with a possible small discrepancy from it. This assumption is motivated by the fact that the GRAVITY Collaboration in 2020 and in 2022 detected Schwarzschild precession in the S2 star orbit around the Supermassive Black Hole (SMBH) at the GC. Using this approach, we were able to constrain parameter $\lambda$ of the potential and, assuming that it represents the graviton Compton wavelength, we also found the corresponding upper bound of graviton mass. The obtained results were then compared with our previous estimates, as well as with the estimates of other authors.


Introduction
Here we use a phenomenological modification of the Newtonian gravitational potential with a non-linear Yukawa-like correction, as provided in [1,2], with the aim of obtaining new constraints on graviton mass from the observed orbits of S-stars around the Galactic Center.
Graviton is the gauge boson of the gravitational interaction and in GR theory is considered as massless, moving along null geodesics at the speed of light c.On the other hand, according to a class of alternative theories, known as theories of massive gravity, gravitational interaction is propagated by a massive field, in which case graviton has some small non zero mass [13][14][15][16][17][18][19][20][21][22][23][24].This approach was first introduced in 1939 by Fierz and Pauli [13].
The LIGO and Virgo collaborations considered a theory of massive gravity to be an appropriate approach and presented their estimate of the mass of the graviton m g < 1.2 × 10 −22 eV in the first publication about the gravitational wave detection from binary black holes [25].Analyzing signals observed during the three observing runs collected in The third Gravitational-wave Transient Catalog (GWTC-3) the LIGO -Virgo -KAGRA collaborations obtained a stringer constraint of m g < 1.27 × 10 −23 eV [26].Various experi- mental constraints on the graviton mass are given in [27].
This paper is organized as follows: in Section 2 we presented the orbital precession in Yukawa-like gravitational potential.In Section 3 we introduced PPN equations of motion, together with the other important expressions that we used for analysis of the stellar orbits around Sgr A* in Yukawa gravity.We then performed an analysis for the potential from [2] that we have previously developed for some other modified potentials, and obtained the results for upper bound on graviton mass in the case of more different S-stars.These results are presented and discussed in section 4. Section 5 is devoted to the concluding remarks.

Orbital precession in Yukawa-like gravitational potential
In order to derive the expression for orbital precession in the gravitational potential (3), we assume that it does not differ significantly from the Newtonian potential Φ N (r) = . Namely, it is well known that orbital precession ∆ϕ per orbital period, induced by small perturbations to the Newtonian gravitational potential which are described by the perturbing potential V(r) = Φ(r) − Φ N (r), could be evaluated as (see e.g.[65] and references therein): where r is related to z via r = L 1 + ez and L = a 1 − e 2 is the semilatus rectum of the orbital ellipse.Approximate formula for orbital precession ∆ϕ Y can be obtained by performing the power series expansion of the perturbing potential V(r) and by inserting the first order term of this expansion into Eq.( 5), which results with: Note that a similar expression for orbital precession was obtained in [2] and it was used for bounding the Compton wavelength and mass of the graviton by the Solar System data.

Stellar orbits in extended/modified PPN formalisms
We simulated the orbits of S-stars around GC using the parameterized post-Newtonian (PPN) equations of motion (for more details about PPN approach see e.g.[107] and references therein).However, it is well known that Yukawa-like potentials could not be entirely represented by the standard PPN formalism and thus require its extension/modification (see the related references in [54]).This is also valid for the potential (3) and its corresponding metric (4).Moreover, since Yukawa gravity is indistinguishable from GR up to the first post-Newtonian correction [57], in addition to the standard PPN equations of motion rGR in GR, PPN equations of motion rY in potential (3) also include an additional term rλ with exponential correction due to the perturbing potential V(r).In this extended PPN formalism (denoted here as PPN Y ), the equations of motion are: where the rN is the Newtonian acceleration, rPPN is the first post-Newtonian correction and rλ is additional Yukawa correction given by the following expressions, respectively: The additional Yukawa correction rλ becomes negligible when λ → ∞, and then rY → rGR , i.e.PPN equations of motion rY in potential (3) reduce to the standard PPN equations of motion in GR.Therefore, the orbits of S-stars in Yukawa gravity and GR can be then simulated by numerical integration of the corresponding expressions (7).The orbital precession in GR is given by the well known expression for the Schwarzschild precession [108]: where a is semi-major axis and e eccentricity of the orbit.Recently the GRAVITY Collaboration detected the orbital precession of the S2 star around the SMBH at GC and showed that it was close to the above prediction of GR [95].For that purpose, they introduced an ad hoc factor f SP in front of the first post-Newtonian correction of GR in order to parametrize the effect of the Schwarzschild metric.In this modified PPN formalism (denoted here as PPN SP ), the equations of motion are given by: The corresponding modified expression for the Schwarzschild precession is [95]: For f SP = 1 the expression (10) reduces to the standard PPN equations of motion rGR in GR given in Eq. ( 7), while the expression (11) also reduces to the corresponding GR prediction from Eq. ( 9).The parameter f SP shows to which extent some gravitational model is relativistic and it is defined as , where β and γ are the post-Newtonian parameters, and in the case of GR the both of them equals to 1 (and thus f SP = 1 in this case).For f SP = 0 the Newtonian case is recovered.However, in the case of S2 star the value of f SP = 1.10 ± 0.

Results: constraints on the Compton wavelength and mass of the graviton
The constraints on parameter λ under which the orbital precession in the gravitational potential (3) deviates from the Schwarzschild precession in GR by a factor f SP , can be obtained provided that the total orbital precession in Yukawa gravity, given by the sum ∆ϕ GR + ∆ϕ Y , is close to the observed precession ∆ϕ SP obtained by the GRAVITY Collaboration: Taking into account the third Kepler law (since the orbits are almost Keplerian): then from Eqs. ( 12) and ( 13) one can obtain the following relation between λ and f SP : The above condition can be used for constraining the Compton wavelength λ of the graviton by the observed values of f SP only in the cases when f SP is larger than 1, since ∆ϕ Y always gives a positive contribution to the total precession in Eq. ( 12).The corresponding constraints on the graviton mass m g can be then found by inserting the obtained value of λ into Eq.( 2).The relative error of the parameter λ (and thus of the graviton mass m g ) in this case can be found by differentiating the logarithm of the above expression (14): It can bee seen that the potential (3) results with the same relative errors as the corresponding Yukawa potential derived in the frame of f (R) theories of gravity (see e.g.[55]).We first used three estimates for f SP obtained by the GRAVITY Collaboration in order to find the corresponding constraints on the Compton wavelength λ of the graviton and its mass m g in the case of S2 star.These are the values of f SP detected by GRAVITY Table 1.The Compton wavelength of the graviton λ, its mass m g , as well as their relative and absolute errors, calculated for three different values of f SP in the case of S2 star.  1, from which it can be seen that the most reliable result was obtained in the case of f SP with the lowest uncertainty, resulting from the fit with the four star data.In that case the relative error for λ and m g was the lowest.In contrast, the worst constraint with unrealistically high relative error was obtained in the second case with the lowest value of f SP = 1.01, due to the fact that it is too close to the corresponding prediction of GR.By comparing the results obtained in the first case with our previous corresponding estimates from Table I in [55], obtained for Yukawa-like potential derived from f (R) theories of gravity, it can be seen that the upper bound on graviton mass m g and its absolute error ∆m g were improved for ∼ 30% in the case of the phenomenological potential (3), although the relative error remained the same.Comparison between the simulated orbits of S2 star during one orbital period, obtained by numerical integration of the equations of motion given by Eq. ( 10) in PPN SP formalism for f SP = 1.10 (blue dotted line) and those given by Eq. ( 7) in PPN Y formalism for λ = 66361.5AU (red dashed line), which corresponds to f SP = 1.10 according to Eq. ( 14).Top right: The differences between the corresponding x and y coordinates in PPN Y and PPN SP formalisms as the functions of time.Middle: Comparison between the simulated orbits of S2 star, near pericenter (left) and near apocenter (right) in two PPN formalisms.Bottom: The same as in top panel, but for five orbital periods.
In the first case of f SP = 1.10 ± 0.19 we also graphically compared the simulated orbits of S2 star, obtained by numerical integration of equations of motion in PPN SP formalism given by Eq. ( 10) with those in PPN Y formalism given by Eq. ( 7) for λ = 66361.5AU which corresponds to this f SP according to Eq. ( 14).The comparisons during one and five orbital periods are presented in the top left and bottom left panels of Fig. 1, respectively, while in the corresponding right panels we presented the differences between the corresponding x and y coordinates in these two PPN formalisms as the functions of time.In order to see more clear the difference between the simulated orbits of S2 star in two PPN formalisms, we show in middle panel of Fig. 1 part of its orbit near the pericenter (left panel) and near the apocenter (right panel).
As it can be seen from Fig. 1, the differences between the simulated orbits of S2 star in these two PPN formalisms are very small, and the maximum discrepancy between the corresponding coordinates during the first orbital period is only ∼ 2 AU at the pericenter.This discrepancy slowly increases with time during the successive orbital periods due to neglecting of the higher order terms in power series expansion of the perturbing potential V(r) in Eq. ( 5).One should also note that the two PPN formalisms give close, but not exactly the same epochs of pericenter passage, as it can be seen from top right panel of Fig. 1.
This sufficiently small difference between the simulated orbits in the two studied PPN formalisms confirms that relation ( 14) could be used for obtaining the constraints on the Compton wavelength λ of the graviton and its mass m g from the latest estimates for f SP obtained by the GRAVITY Collaboration.
Taking the above considerations into account, we then estimated the Compton wavelength λ of the graviton, its mass m g , and also their relative and absolute errors for all S-stars from Table 3 in [85] except for S111 star.For that purpose, and in order to avoid the cases when f SP < 1, we adopted the same strategy as in [55] and assumed that the recent GRAVITY estimates for f SP are very close to the value in GR of f SP = 1, so that these estimates represent a confirmation of GR within 1σ.Therefore, we constrained the graviton mass m g in the particular cases when f SP = 1 + ∆ f SP ± ∆ f SP , i.e. for f SP = 1.19 ± 0.19, f SP = 1.16 ± 0.16 and f SP = 1.144 ± 0.144.As before, we used the expressions ( 14), (15)  and (2) for this purpose, and the obtained results are presented in Tables 2 and 3. Besides, Table 2 also contains the results for f SP = 1.10 ± 0.19 since, as shown in Table 1, this estimate is sufficiently larger than 1.
Using data from these tables, in Fig. 2 we give the comparison of the estimated Compton wavelength λ of the graviton, as well as for graviton mass upper bound, for four stars (S2, S29, S38, S55) which the GRAVITY collaboration used for the newest estimation of f SP .As it can be seen from Fig. 2, all constraints in the case of S2, S38 and S55 star are approximately of the same order of magnitude (λ ∼ 10 5 AU and m g ∼ 10 −22 eV).The only exception is S29 star since it results with an order of magnitude larger values of λ, and hence an order of magnitude smaller estimates for upper bound on graviton mass m g .This is not surprising because S29 star has much longer orbital period of ∼ 101 yr with respect to the orbital periods of the other three stars which are between ∼ 13 − 20 yr (see Table 3 from [85]).
If one compares these results with our previous corresponding estimates from Tables I and III in [55], obtained for Yukawa-like potential derived from f (R) theories of gravity, it can be noticed that both upper bound on graviton mass m g and its absolute error in the case of Yukawa-like potential (3) were further improved for ∼ 30%, in a similar way as for S2 star in the first case from Table 1 (i.e. for f SP = 1.10 ± 0.19).
Although the current GRAVITY estimates of f SP = 1.10 ± 0.19 (from [95]) and f SP = 0.85 ± 0.16 and f SP = 0.997 ± 0.144 (from [109]) can improve our previous constraints on the upper bound of graviton mass for about ∼ 30% (these results can be compared with our previous corresponding estimates from Tables I and III in [55] for the corresponding Sstar), we have to stress that we take assumption that f SP has been measured for all S-stars orbits already to a given precision.In reality, it is expected that the orbits of different S-Table 2. The Compton wavelength of the graviton λ, its mass m g , as well as their relative and absolute errors, calculated for f SP = 1.10 ± 0.19 and f SP = 1.19 ± 0.19 in the case of all S-stars from Table 3  stars should result with slightly different measured values and accuracies of f SP .Because of that, our assumption that f SP is the same for all S-stars (see Tables 2 and 3) probably does not hold.
The values of the graviton Compton wavelength λ, SMBH mass M, distance R to the GC and the osculating orbital elements a, e, i, Ω, ω, P, T which correspond to the minimum of χ 2 red are found using the differential evolution optimization method, implemented as Python Scipy function scipy.optimize.differential_evolution.This is a population-based metaheuristic search technique of finding the global minimum of a multivariate function which is especially suitable in evolutionary computations, since it is stochastic in nature, does not use gradient descent to find the minimum, can search large areas of candidate space and seeks to iteratively enhances a candidate solution concerning a specified quality metric.In order to improve the minimization slightly, the final result of the differential evolution optimization is further polished at the end using Python Scipy function scipy.optimize.minimize.This also results with an approximation for the inverse Hessian matrix of χ 2 red , which on the other hand could be considered as a good estimation for the covariance matrix of the parameters.Therefore, the standard error for each fitted parameter can be calculated by taking the square root of the respective diagonal element of this covariance matrix.
A particular value of χ 2 red which corresponds to some specific combination of the mentioned adjustable parameters is calculated in the following way: 1.
First, a simulated orbit of S2 star in the extended PPN Y formalism is calculated by numerical integration of the equations of motion (7), starting from initial conditions (x 0 , y 0 , ẋ0 , ẏ0 ), where the first two components specify the initial position and the last two the initial velocity in the orbital plane.In our simulations the initial conditions correspond to the time of apocenter passage t ap preceding the first astrometric observation at t 0 : t ap = T − (2k − 1) P 2 , where T is the time of pericenter passage, P is the orbital period and k is the smallest positive integer (number of periods) for which t ap ≤ t 0 .Then, the initial conditions are: x 0 = −r ap , y 0 = 0, ẋ0 = 0 and ẏ0 = −v ap , where r ap = a(1 + e) is the apocenter distance and v ap = 2π a P 1 − e 1 + e is the corresponding orbital velocity at the apocenter.

2.
The true orbit obtained in the first step, which depends only on a, e, P, T, was then projected to the observer's sky plane using the remaining geometrical orbital elements i, Ω, ω, in order to obtain the corresponding positions (x c , y c ) along the apparent orbit: where A, B, F, G are the Thiele-Innes elements: In addition, the radial velocities V rad are obtained from the corresponding true positions (x, y) and orbital velocities ( ẋ, ẏ) as (see e.g.[62] and references therein): where θ = arctan y x .

3.
Finally, χ 2 red is obtained according to Eq. ( 16), taking into account only those apparent positions (x c , y c ) which are calculated at the same epochs as the astrometric observations (x o , y o ).The obtained results of the orbital fitting in the case of S2 star are presented in Fig. 3, and the corresponding best-fit values of the parameters are given in Table 4.is a bit larger, but still within the error intervals of the corresponding values from Tables 2 and 3, obtained according to Eq. ( 14) from the detected values of f SP .Therefore, the results of direct orbital fitting are in agreement with our constraints on the graviton Compton wavelength λ and mass m g presented in Tables 2 and 3, which were obtained from the values of f SP estimated by the GRAVITY Collaboration.

Conclusions
Here we used the phenomenological Yukawa-like gravitational potential from [1,2] in order to obtain the constraints on the graviton mass m g from the detected Schwarzschild precession in the observed stellar orbits around the SMBH at GC.For that purpose we used two modified/extended PPN formalisms in order to derive the relation between the Compton wavelength λ of the graviton and parameter f SP which parametrizes the effect of the Schwarzschild metric, and which was obtained by the GRAVITY Collaboration from the observed stellar orbits at GC.The results from this study can be summarized as follows: 1.
We found the condition for parameter λ of the phenomenological Yukawa-like gravitational potential (3) under which the orbital precession in this potential deviates from the Schwarzschild precession in GR by a factor f SP ; 2.
The relation ( 14) derived from the phenomenological potential (3) in the frame of the two modified/extended PPN formalisms could be used for obtaining the reliable constraints on the graviton mass m g from the latest estimates for f SP by the GRAVITY Collaboration in the cases when f SP > 1; 3.
Both studied PPN formalisms result with close and very similar simulated orbits of S-stars, which practically overlap during the first orbital period and then begin to slowly diverge over time due to some assumed theoretical approximations; 4.
In most cases, the constraints on the upper bound on graviton mass m g and its absolute error ∆m g , obtained using the phenomenological potential (3) were improved for ∼ 30% in respect to our previous corresponding estimates from [55], obtained using a slightly different Yukawa-like potential derived in the frame of f (R) theories of gravity, although the relative errors in both cases remained the same; 5.
These results were also confirmed in the case of S2 star by fitting of its observed orbit in the frame of the extended PPN Y formalism, which resulted with the bestfit value for the graviton Compton wavelength λ within the error intervals of its corresponding estimates obtained according to Eq. ( 14) from the detected values of f SP ; 6.
The least reliable constraints with unrealistically high uncertainties were only obtained from the estimates for f SP which were very close to its value predicted by GR, being just slightly larger than 1; 7.
If one compares the results from Table 2 with those from Table 3, it can be seen that the upper bounds on graviton mass m g are very similar.In the case of S2 star m g < (1.5 ± 0.8) × 10 −22 eV and relative error is around 50%.We can conclude that the more precise future observations are required in order to further improve the upper graviton mass bounds.

Figure 1 .
Figure 1.Top left:Comparison between the simulated orbits of S2 star during one orbital period, obtained by numerical integration of the equations of motion given by Eq.(10) in PPN SP formalism for f SP = 1.10 (blue dotted line) and those given by Eq. (7) in PPN Y formalism for λ = 66361.5AU (red dashed line), which corresponds to f SP = 1.10 according to Eq. (14).Top right: The differences between the corresponding x and y coordinates in PPN Y and PPN SP formalisms as the functions of time.Middle: Comparison between the simulated orbits of S2 star, near pericenter (left) and near apocenter (right) in two PPN formalisms.Bottom: The same as in top panel, but for five orbital periods.

19 Figure 2 .
Figure 2. Left: Constraints on the Compton wavelength λ of the graviton from the orbits of S2, S29, S38 and S55 star, in case of the following f SP estimates: 1.1, 1.144, 1.16 and 1.19.Right: Upper bounds on graviton mass m g , for the same group of S-stars and the same f SP .

Figure 3 .
Figure 3. Left: Comparison between the best-fit orbit of the S2 star (blue solid line), simulated in the extended PPN Y formalism, and the corresponding astrometric observations from [85] (black circles with error bars).Right: The same for the radial velocity of the S2 star (top), as well as for its α (middle) and δ (bottom) offset relative to the position of Sgr A* at the coordinate origin.Red dots in the right panels denote the corresponding O-C residuals.
19was obtained by the GRAVITY Collaboration, indicating the orbital precession of f SP × 12.
′ 1 in its orbit around Sgr A* [95].Recently, this collaboration updated the above first estimate to f SP = 0.85 ± 0.16, also obtained from the detected Schwarzschild precession of S2 star orbit [109].Besides, they also presented the following estimate obtained from the fit with the four star (S2, S29, S38, S55) data with the 1σ uncertainty: f SP = 0.997 ± 0.144 [109].

Table 4 .
Best-fit values of the graviton Compton wavelength λ, SMBH mass M, distance R to the GC and the osculating orbital elements a, e, i, Ω, ω, P, T of the S2 star orbit.Concerning the graviton Compton wavelength λ obtained from S2 star orbit, it can be seen that its best-fit value of λ ≈ 8.2 × 10 4 AU from Table