LETTER TO THE EDITOR

# Exploring the presence of a fifth force at the Galactic Center

GRAVITY Collaboration*: K. Abd El Dayem¹, R. Abuter⁴, N. Aimar¹⁰,⁷, D. P. Amaro Seoane¹⁴,²,¹⁸, A. Amorim⁸,⁷, D. J. P. Berger³,⁴, D. H. Bonnet⁴, G. Bourdarot², W. Brandner⁵, D. V. Cardoso⁷,¹⁵, D. Y. Clénet¹, D. R. Davies², D. P. T. de Zeeuw¹⁹, D. A. Drescher², D. A. Eckart⁶,¹³, D. F. Eisenhauer²,¹⁷, H. Feuchtgruber², D. G. Finger², N. M. Förster Schreiber², D. A. Foschi¹,², **, D. P. Garcia¹⁰,⁷, D. E. Gendron¹, D. R. Genzel²,¹¹, S. Gillessen², D. M. Hartl², X. Haubois⁹, D. F. Haussmann², T. Henning⁵, D. S. Hippler⁵, D. M. Horrobin⁶, L. Jochum⁹, L. Jocou³, D. A. Kaufer⁹, D. P. Kervella¹, D. S. Lacour¹,⁴, D. V. Lapeyrère¹, J.-B. Le Bouquin³, D. P. Léna¹, D. D. Lutz², D. F. Mang², D. N. More², J. Osorno¹, D. T. Ott², T. Paumard¹, D. K. Perraut³, D. G. Perrin¹, D. S. Rabien², D. C. Ribeiro², M. Sadun Bordoni², D. S. Scheithauer⁵, D. J. Shangguan²⁰, D. T. Shimizu², D. J. Stadler¹²,², O. Straub²,¹⁶, D. C. Straubmeier⁶, D. E. Sturm², D. L. J. Tacconi², I. Urso¹, D. F. Vincent¹, D. S. D. von Fellenberg¹³,², D. E. Wieprecht², and J. Woillez⁴

*(Affiliations can be found after the references)*

Received 21 March 2025 / Accepted 1 May 2025

## ABSTRACT

**Aims.** We investigate the presence of a Yukawa-like correction to Newtonian gravity at the Galactic Center, leading to a new upper limit on the intensity of such a correction.

**Methods.** We performed a Markov chain Monte Carlo (MCMC) analysis using the astrometric and spectroscopic data of star S2 collected at the Very Large Telescope by GRAVITY, NACO, and SINFONI instruments, covering the period from 1992 to 2022.

**Results.** The precision of the GRAVITY instrument allows us to derive the most stringent upper limit at the Galactic Center for the intensity of the Yukawa contribution ($\propto ae^{-\lambda}$) of $|a| < 0.003$ for a scale length of $\lambda = 3 \cdot 10^{13}$ m (~200 AU). This is an improvement on all estimates obtained in previous works by roughly one order of magnitude.

**Key words.** gravitation – celestial mechanics – Galaxy: center

## 1. Introduction

General relativity (GR) is the most widely recognized theory of gravity today. Its predictions have been extensively tested on Solar System scales and using gravitational waves emission by black holes (BHs) and binary pulsars (Will 2014, 2018a; Nitz et al. 2021). Until now, no significant deviation from GR has been detected in any of these observations. However, it is also known that beyond the regime we are currently able to test experimentally, GR is generally considered insufficient.

Cosmological observations indicate an expanding Universe whose acceleration can only be explained by introducing a cosmological constant ad hoc, which lacks a solid theoretical explanation and raises several issues (Weinberg 1989; Peebles & Ratra 2003). Other observational evidence, such as the rotational curve of galaxies (van Albada et al. 1985; Salucci 2019) or gravitational lensing effects (Massey et al. 2010) indicate the presence of a dark massive component of the Universe, whose nature is still unknown. Furthermore, it is well known that GR lacks a quantum description at high-energy scales and several attempts have been made to create a theory valid at all scales (for a review on the state-of-the-art approaches, see, e.g., Esposito 2011; Kiefer 2023).

One way to address these inconsistencies between theory and experiments is to directly modify GR, giving rise to a plethora of possible extended theories of gravity (ETG). In particular, a Yukawa-like interaction emerges quite naturally in the weak field limit of several ETGs; for instance, scalar-tensor-vector theories (Moffat 2006), massive gravity theories (Visser 1998; Hinterbichler 2012), theories in higher dimensions with Kaluza-Klein compactification (Bars & Visser 1986; Hoyle et al. 2001), massive Brans-Dicke theories (Perivolaropoulos 2010; Alsing et al. 2012), or $f(R)$ theories (Capozziello et al. 2015). However, the so-called “fifth-force” scenario has also made an appearance in a number of specific dark matter models (Frieman & Gradwohl 1991; Gradwohl & Frieman 1992; Carroll et al. 2009).

Due to the importance that a modification of Newtonian gravity would have on our understanding of the Universe, the presence of a Yukawa-like contribution has been repeatedly investigated in the past. The fifth-force intensity is well constrained at Solar System scales via the motion of planets (Konopliv et al. 2011; Hees et al. 2014; Bergé 2017; Will 2018b; Shankaranarayanan & Johnson 2022), from the Lunar Laser Ranging experiment (Hofmann & Müller 2018) and also by making use of the planetary ephemerides (Mariani et al. 2023; Fienga & Minazzoli 2024). Recent constraints have been obtained from asteroids tracking (Tsai et al. 2023, 2024) and test of the weak equivalence principle (Touboul 2022).

The discovery of orbiting stars around the Galactic Center (GC) (Eckart & Genzel 1996; Schödel et al. 2002; Ghez et al.

<footnote>* GRAVITY is developed in collaboration by MPE, LESIA of Paris Observatory/CNRS/Sorbonne Université/Univ. Paris Diderot, and IPAG of Université Grenoble Alpes/CNRS, MPIA, Univ. of Cologne, CENTRA – Centro de Astrofísica e Gravitação, and ESO.</footnote>

<footnote>** Corresponding author: arianna.foschi@obspm.fr</footnote>
2003; Gillessen et al. 2009a,b; Sabha 2012), all located within one arcsecond distance from the supermassive black hole (SMBH) Sagittarius A* (Sgr A*), allows us to test GR in a completely different environment from the Solar System.

The importance of looking for a fifth force in the GC lies in the fact that many ETGs that predict a Yukawa-like term also display a screening mechanism that suppresses the fifth force contribution at Solar System scales and prevents its detection. This would explain why it would be yet unobserved, while its effect may be different around SMBHs.

The current constraints on the intensity of a fifth force in the GC come from the analysis of SgrA**'s shadow by the Event Horizon Telescope (Vagnozzi 2023), from the measurement of the Schwarzschild precession in S2 motion (GRAVITY Collaboration 2020; Jovanović et al. 2023, 2024a,b) and from the analysis of S-stars publicly available (or mock) data (Borka et al. 2013, 2021; Capozziello et al. 2014; Zakharov et al. 2016, 2018; de Martino 2021; Della Monica et al. 2022), while also including the presence of a (expected) bulk mass distribution around Sgr A* (Jovanović et al. 2021). For a complete and comprehensive review on astrophysical and theoretical constraints at the GC, we refer to de Laurentis et al. (2023).

In the context of S-stars, a previous work by Hees et al. (2017), with a full analysis of the S2 data, showed that the intensity of such a contribution cannot exceed $\alpha \sim 0.01$ at scales comparable to the S2-SgrA* distance.

In this Letter, we use the astrometric and spectroscopic measurements of the star S2 collected at the Very Large Telescope (VLT) by GRAVITY, NACO, and SINFONI to constrain the intensity of a possible Yukawa correction at the GC. Although we do not expect our results to consistently deviate from the estimates obtained in the aforementioned literature, a complete analysis of S2 including GRAVITY data, which dominate the $\chi^2$ due to their very small uncertainties, is still lacking. As we show here, the precision of the GRAVITY instrument allows us to place a significantly stronger constraint than the previous estimates.

## 2. Observations

The set of available data, $D$, can be organized according to the following criteria:

(a) Astrometric data Dec, RA

- 128 data points collected using both the SHARP camera at the New Technology Telescope between 1992 and 2002 (~10 data points, accuracy of ≈4 mas) and the NACO imager at the VLT between 2002 and 2019 (118 data points, accuracy of ≈0.5 mas);

- 76 data points collected by GRAVITY at the VLT between 2016 and April 2022 (accuracy of ≈50 µas).

(b) Spectroscopic data $V_R$

- 102 data points collected by SINFONI at the VLT (100 points) and NIRC2 at Keck (2 points) collected between 2000 and March 2022 (accuracy in good conditions of ≈10–15 km/s).

## 3. Yukawa correction to Newtonian force

### 3.1. Model

The potential we aim to test takes the following form:

$$U = -\frac{GM}{r} (1 + |\alpha|e^{-r/\lambda}), \quad (1)$$

where $\alpha$ represents the strength of interaction and $\lambda$ is a scale parameter, which depends on the specific theory considered. For example, when new massive fields are included in the theory, $\lambda$ represents the Compton wavelength of the field, which is related to the mass by $m_\phi = h/c\lambda$, where $h$ is the Planck constant.

In comparison to previous works, here the Schwarzschild precession in the S2 orbit is also included, since it has been formally detected at a $10\sigma$ confidence level by the GRAVITY Collaboration (GRAVITY Collaboration 2020, 2024).

Although a formal parametrized post-Newtonian (PN) treatment is not possible when a massive field is included in the action (Alsing et al. 2012; Poisson & Will 2012), we can still derive the equations of motion of a test particle, assuming the parametrized PN parameters to be $\gamma = \beta = 1$. This latter assumption is valid for ETGs that are indistinguishable from GR at 1PN order and it is supported by different experimental observations, including those carried out for the GC (see, e.g., Will 2018a; Hofmann & Müller 2018; GRAVITY Collaboration 2020; Iorio 2024). On the other hand, non-metric theories of gravity or specific subclasses of f(R) and scalar-tensor theories, with parametrized PN parameters that significantly stray from unity, are excluded from this analysis (Will 2014).

The total acceleration experienced by the star is

$$a_{\text{TOT}} = a_{\text{New}} + a_{\text{Yuk}} + a_{\text{1PN}}, \quad (2)$$

where $a_{\text{New}} + a_{\text{Yuk}}$ are derived from the potential in Eq. (1) and

$$a_{\text{1PN}} = \frac{GM}{c^2 r^2} \left[ \left( \frac{4GM}{r} - v^2 \right) \frac{\mathbf{r}}{r} + 4i\mathbf{v} \right], \quad (3)$$

with $\mathbf{r} = r\hat{\mathbf{r}}, \mathbf{v} = (i\hat{\mathbf{r}}, r\hat{\theta}\hat{\phi}, r\hat{\phi}\sin\theta\hat{\phi})$ and $\mathbf{v} = |\mathbf{v}|$.

The above expression coincides with the 1PN acceleration derived in Alves et al. (2024) for the two-body problem in massive Brans-Dicke theory. Section 3.4 will be devoted to a comparison of Eq. (3) with results developed in the literature when extra massive degrees of freedom are included in the theory.

### 3.2. Method

The numerical integration of the equations of motion was performed using a Runge-Kutta 4(5) method, with further details reported in Appendix A. During the fit of the S2 data, the Yukawa length scale, $\lambda$, was kept fixed, with values set between $10^{12} \le \lambda \le 10^{15}$ m, while the intensity, $\alpha$, was allowed to vary together with other parameters describing the system.

Specifically, the set of parameters is given by:

$$\Theta_i = \{e, a_{\text{sma}}, \Omega_{\text{orb}}, i_{\text{orb}}, \omega_{\text{orb}}, t_p, R_0, M, x_0, y_0, v_{x_0}, v_{y_0}, v_{z_0}, \alpha\}, \quad (4)$$

where $e$ is the eccentricity and $a_{\text{sma}}$ the semi major axis of the star S2, while $\Omega_{\text{orb}}$, $i_{\text{orb}}$, and $\omega_{\text{orb}}$ are the three angles used to project the star's orbital frame into the observer reference frame using the procedure reported in Appendix B.1. Then, $t_p$ is the time of pericenter passage, while $M$ and $R_0$ are the SMBH mass and the GC distance, respectively. The additional parameters $\{x_0, y_0, v_{x_0}, v_{y_0}, v_{z_0}\}$ characterize the NACO/SINFONI data reference frame with respect to Sgr A* (Plewa et al. 2015).

To fit the S2 data, we perform a Markov chain Monte Carlo (MCMC) analysis using the Python package EMCEE (Foreman-Mackey et al. 2013). The log-likelihood is given by

$$\ln \mathcal{L} = \ln \mathcal{L}_{\text{pos}} + \ln \mathcal{L}_{\text{vel}}, \quad (5)$$
<table><caption>**Table 1.** Uniform priors used in the MCMC analysis.</caption><thead><tr><th>Parameter</th><th>Θ<sup>0</sup><sub>i</sub></th><th>Lower bound</th><th>Upper bound</th></tr></thead><tbody><tr><td>e</td><td>0.88441</td><td>0.83</td><td>0.93</td></tr><tr><td>a<sub>sma</sub> ["]</td><td>0.12497</td><td>0.119</td><td>0.132</td></tr><tr><td>i<sub>orb</sub> [°]</td><td>134.69241</td><td>100</td><td>150</td></tr><tr><td>ω<sub>orb</sub> [°]</td><td>66.28411</td><td>40</td><td>90</td></tr><tr><td>Ω<sub>orb</sub> [°]</td><td>228.19245</td><td>200</td><td>250</td></tr><tr><td>t<sub>p</sub> [yr]</td><td>2018.37902</td><td>2018</td><td>2019</td></tr><tr><td>M [10<sup>6</sup> M<sub>☉</sub>]</td><td>4.29950</td><td>4.1</td><td>4.8</td></tr><tr><td>R<sub>0</sub> [10<sup>3</sup> pc]</td><td>8.27795</td><td>8.1</td><td>8.9</td></tr></tbody></table>

<table><caption>**Table 2.** Gaussian priors used in the MCMC analysis.</caption><thead><tr><th>Parameter</th><th>Θ<sup>0</sup><sub>i</sub></th><th>ξ</th><th>σ</th></tr></thead><tbody><tr><td>x<sub>0</sub> [mas]</td><td>-0.244</td><td>-0.055</td><td>0.25</td></tr><tr><td>y<sub>0</sub> [mas]</td><td>-0.618</td><td>-0.570</td><td>0.15</td></tr><tr><td>v<sub>x<sub>0</sub></sub> [mas/yr]</td><td>0.059</td><td>0.063</td><td>0.0066</td></tr><tr><td>v<sub>y<sub>0</sub></sub> [mas/yr]</td><td>0.074</td><td>0.032</td><td>0.019</td></tr><tr><td>v<sub>z<sub>0</sub></sub> [km/s]</td><td>-2.455</td><td>0</td><td>5</td></tr></tbody></table>

**Notes.** $\xi$ and $\sigma$ represent the mean and the standard deviation of the Gaussian distributions, respectively, which come from Plewa et al. (2015).

where

$$ \ln \mathcal{L}_{\text{pos}} = - \sum_{i=1}^{N} \left[ \frac{(\text{Dec}_i - \text{Dec}_{\text{model},i})^2}{\sigma_{\text{Dec}_i}^2} + \frac{(\text{RA}_i - \text{RA}_{\text{model},i})^2}{\sigma_{\text{RA}_i}^2} \right], \quad (6) $$

and

$$ \ln \mathcal{L}_{\text{vel}} = - \sum_{i=1}^{N} \frac{(V_{R,i} - V_{\text{model},i})^2}{\sigma_{V_{R,i}}^2}. \quad (7) $$

The priors are listed in Tables 1 and 2. Uniform priors were used for the physical parameters, that is, we only imposed physically motivated bounds, while Gaussian priors were implemented for the offset parameters, since the latter have been well constrained by an independent previous work and were not expected to change (Plewa et al. 2015). The initial points $Θ_i^0$ in the MCMC are chosen to be the most recent best fit parameters of S2 orbit reported in the literature (GRAVITY Collaboration 2024).

In the sampling phase of the MCMC implementation, we used 64 walkers and $10^5$ iterations. The burning-in phase was skipped and the last 80% of the chains was used to compute the mean and standard deviation of the posterior distributions of the parameters. The convergence of the MCMC analysis was ensured by means of the autocorrelation time $\tau_c$, that is, we ran $N$ iterations such that $N \gg 50 \tau_c$.

### 3.3. Results

We can classify three different regimes in the posterior distributions $P(|\alpha|D)$, according to the value of $\lambda$ with respect to the orbital range of S2, which is $1.7 \cdot 10^{13} \text{ m} \lesssim r_{s2} \lesssim 1.5 \cdot 10^{14} \text{ m}$. When $\lambda \ll r_{s2}$, the acceleration is no longer dependent on the parameter $\alpha$ and no meaningful constraints can be obtained in this regime. The small difference in the 95% upper limit on

$|\alpha|$ with the UCLA group resides in the different model implemented to fit the data.

In particular, Hees et al. (2017) considered no Schwarzschild precession but an extended mass with power law distribution instead, which is absent from our work. The reason behind this choice comes from the fact that the presence of a spherically symmetric mass distribution around SgrA* has been extensively tested by the GRAVITY Collaboration (GRAVITY Collaboration 2018, 2022, 2024), finding a stringent upper limit of $M_{\text{ext}} \lesssim 1200 M_{\odot}$. Taking into account this current upper limit, we decided to neglect the presence of an extended mass, as it would not alter our results.

When $\lambda \sim r_{s2}$, the best constraints on $|\alpha|$ can be obtained, determining the most stringent upper limit of $|\alpha| < 0.003$ for $\lambda = 3 \cdot 10^{13} \text{ m} \sim 200 \text{ AU}$. This limit improves the previous estimate of Hees et al. (2017), which reported $|\alpha| < 0.016$ for $\lambda = 150 \text{ AU} \sim 2.2 \cdot 10^{13} \text{ m}$, stressing the importance of the precision of the GRAVITY instrument.

Finally, when $\lambda \gg r_{s2}$, the only component left in the equations of motion is the monopolar term $M(1+\alpha)/r$. This corresponds to a simple rescaling of the mass term and, hence, in this regime $M$ and $\alpha$ are completely degenerate and they cannot be constrained separately. To obtain the upper limit on $\alpha$ in this regime, the bounds on $M$ in Table 1 have been extended to $M \in (10^{-4}, 10^4) \cdot 10^6 M_{\odot}$, and the same bounds have been used for $\alpha$.

A summary of the above results is reported in Figure 2, where the 95% confidence interval on $|\alpha|$ as function of $\lambda$ is shown. Those confidence intervals are estimated as three times the standard deviation when the posterior distributions are normal or by computing the upper limit that corresponds to 95% of the area below the curve when the distributions have different shapes. We note that GRAVITY data produce an overall improvement of roughly one order of magnitude over the entire parameter space tested.

If we assume that the gravitational interaction is mediated by a massive boson as in massive gravity theories (where $\alpha = 1$), the length scale $\lambda$ corresponds to the Compton wavelength of the particle; hence, an upper limit on the graviton's mass can be derived. Since $\alpha = 1$ is excluded at 95% confidence level for $\lambda \lesssim 8 \cdot 10^{14} \text{ m}$, this lower bound on the wavelength, $\lambda$, can be translated into an upper limit on the graviton mass, corresponding to $m_g \lesssim 2.5 \cdot 10^{-22} \text{ eV}$.

In the regime where $\lambda \sim r_{s2}$, when no 1PN acceleration is included, the presence of the Yukawa term in the equation of motion induces a prograde precession comparable to the Schwarzschild one. This is shown in Figure 1 for $\lambda = 10^{13} \text{ m}$, where we can see that the posterior distribution $P(\alpha|D)$ is a Gaussian with mean around $\alpha \sim 0.0026$, as opposed to the 1PN posterior.

Following Adkins & McDonnell (2007), we can compute the precession angle induced by a potential in a full orbit as

$$ \Delta\phi_p = -\frac{2L}{GMe^2} \int_{-1}^{1} \frac{dz \, z}{\sqrt{1-z^2}} \frac{dU(z)}{dz}, \quad (8) $$

where $U(z)$ is the perturbing potential evaluated at radius $r = L/(1+ez)$ with $L = a_{\text{sma}}(1-e^2)$. For $\alpha = 0.0026$, this corresponds to $\Delta\phi_p \sim 0.13^\circ$. Taking the most up-to-date value of the Schwarzschild precession reported in GRAVITY Collaboration (2024), $\Delta\phi_{\text{Sch}} = 12.1' \times (0.911 \pm 0.13) = (0.18 \pm 0.03)^\circ$, we can see that the precession angle induced by the Yukawa potential is compatible, within the $2\sigma$ uncertainties, with $\Delta\phi_{\text{Sch}}$.
<figure><img src="image_1.png" /><figcaption>**Fig. 1.** Comparison of the posterior distributions $P(|\alpha||D)$ between the Keplerian model (black curve) and the 1PN model (red curve).</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>**Fig. 3.** Posterior probability density of $|\alpha|$ using the 1PN expression in Eq. (3) (black curve) versus the 1PN expansion for $f(R)$ gravity derived in Tan & Lu (2024) (red curve), when $\lambda = 3 \cdot 10^{13}$ m.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** 95% confidence level on $|\alpha|$ obtained in this work (red dots), compared with previous estimates by Hees et al. (2017) (blue dots). The dotted vertical line represents the minimum of the curve, in correspondence to $\lambda \sim 3 \cdot 10^{13}$ m ~ 200 AU.</figcaption></figure>

### 3.4. Comparison with theoretical estimates

As stated in the previous section, the inclusion of the 1PN acceleration in the equations of motion implies an additional assumption; namely, that any correction to the GR expression is subleading with respect to Eq. (3) and, hence, negligible in our fit. To show this, we compared our results obtained using Eq. (3) with the analytical expressions derived in the literature for some specific theories.

Alves et al. (2024) derived the 1PN acceleration in massive Brans-Dicke theory, which results in exactly the same expression as Eq. (3), as long as we identify $\alpha = a_1/2$, where $a_1 = 2\varphi/(2\omega_0 + 3)$, and set the reduced mass to $\eta = 0$, which is clearly a good approximation for the SgrA\*-S2 system.

In Tan & Lu (2024) an analytical expression for the 1PN acceleration in $f(R)$ gravity is obtained and reported in their Eq. (17). We used this expression to show that our upper limits on $\alpha$ are not affected by this difference in the 1PN expansion, at least when the uncertainty on $\alpha$ is the smallest (i.e., for $\lambda = 3 \cdot 10^{13}$ m). In Fig. 3, the comparison between the posteriors is reported, showing that the upper limit on $\alpha$ is only changed by a factor of 2 when the full expression for $f(R)$ gravity is used.

Tan & Lu (2024) also showed that the use of a multiple-star fit (specifically including S2, S29, and S55) with the 1PN acceleration for $f(R)$ gravity derived in Eq. (17) could potentially break the degeneracy between $M$ and $\alpha$ in the large $\lambda$ limit. This could produce a stringent upper limit on the fifth force intensity for $\lambda > 10^{15}$ m as well.

In addition, Losada et al. (2025) demonstrated that combining S2 motion with S62 lensing observations can also potentially break the degeneracy between the parametrized PN parameters $\beta$ and $\gamma$. However, we leave a multi-star analysis for a future work.

## 4. Conclusions

In this paper, we update the current constraints on the fifth force intensity at the GC using GRAVITY data for S2 from 2017 to 2022, including the pericenter passage. These data have allowed us to significantly improve previous estimates on the same effect, giving a 95% confidence level curve that is one order of magnitude below the previous estimates (see Figure 2). Specifically, three different behaviors are found in the posterior distribution of $\alpha$, according to the value of the length scale of the new interaction $\lambda$ compared to S2 orbital range.

The minimum value of $\alpha$ is found for $\lambda = 3 \cdot 10^{13}$ m (~200 AU) where $|\alpha| < 0.003$. For comparison, in correspondence with the minimum found by Hees et al. (2017), $\lambda = 150$ AU ~ $2.2 \cdot 10^{13}$ m, we found $|\alpha| < 0.0031$. We also showed that the 1PN expansion used in this work coincides with the expression developed for massive Brans-Dicke theory and that additional terms proportional to $\alpha$ in the 1PN acceleration for $f(R)$ theories are subdominant with respect to the expression used in this work and thereby negligible, as the upper limit on $|\alpha|$ is unaffected.

A complete analysis including all S-stars, specifically those with apocenter passage farther away than S2, is left for future work, with the aim of further improving the confidence level curve. This would allow us to possibly expand the range of $\lambda$ and obtain additional and meaningful constraints at $\lambda \ge 8 \cdot 10^{14}$ m.

*Acknowledgements.* A.F. and F.V. would like to thank Aurélien Hees and Laura Bernard for useful comments and fruitful discussions during the preparation of this work. We are very grateful to our funding agencies (MPG, ERC, CNRS [PNCG, PNGRAM], DFG, BMBF, Paris Observatory [CS, PhyFOG], Observatoire des Sciences de l'Université de Grenoble, and the Fundação para a
Ciência e a Tecnologia), to ESO and the Paranal staff, and to the many scientific and technical staff members in our institutions, who helped to make NACO, SINFONI, and GRAVITY a reality. This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 101007855. We acknowledge the financial support provided by FCT/Portugal through grants 2022.01324.PTDC, PTDC/FIS-AST/7002/2020, UIDB/00099/2020 and UIDB/04459/2020. J.S. acknowledge the National Science Foundation of China (12233001) and the National Key R&D Program of China (2022YFF0503401).

## References

Adkins, G. S., & McDonnell, J. 2007, Phys. Rev. D, 75, 082001

Alsing, J., Berti, E., Will, C. M., & Zaglauer, H. 2012, Phys. Rev. D, 85, 064041

Alves, M. F. S., Toniato, J. D., & Rodrigues, D. C. 2024, Phys. Rev. D, 109, 044045

Bars, I., & Visser, M. 1986, Phys. Rev. Lett., 57, 25

Bergé, J. 2017, in 52nd Rencontres de Morion on Gravitation 191

Borka, D., Jovanović, P., Jovanović, V. B., & Zakharov, A. F. 2013, JCAP, 11, 050

Borka, D., Jovanović, V. B., Capozziello, S., Zakharov, A. F., & Jovanović, P. 2021, Universe, 7, 407

Capozziello, S., Borka, D., Jovanović, P., & Jovanović, V. B. 2014, Phys. Rev. D, 90, 044052

Capozziello, S., Harko, T., Koivisto, T. S., Lobo, F. S. N., & Olmo, G. J. 2015, Universe, 1, 199

Carroll, S. M., Mantry, S., Ramsey-Musolf, M. J., & Stubbs, C. W. 2009, Phys. Rev. Lett., 103, 011301

Catanzarite, J. H. 2010, arXiv e-prints [arXiv:1008.3416]

de Laurentis, M., De Martino, I., & Della Monica, R. 2023, Rept. Prog. Phys., 86, 104901

de Martino, I., & della Monica, R., & de Laurentis, M., 2021, Phys. Rev. D, 104, L101502

Della Monica, R., de Martino, I., & de Laurentis, M. 2022, MNRAS, 510, 4757

Eckart, A., & Genzel, R. 1996, Nature, 383, 415

Esposito, G. 2011, arXiv e-prints [arXiv:1108.3269]

Fienga, A., & Minazzoli, O. 2024, Living Rev. Rel., 27, 1

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Frieman, J. A., & Gradwohl, B.-A. 1991, Phys. Rev. Lett., 67, 2926

Ghez, A. M., Duchêne, G., Matthews, K., et al. 2003, ApJ, 586, L127

Gillessen, S., Eisenhauer, F., Trippe, S., et al. 2009a, ApJ, 692, 1075

Gillessen, S., Eisenhauer, F., Fritz, T. K., et al. 2009b, ApJ, 707, L114

Gradwohl, B.-A., & Frieman, J. A. 1992, ApJ, 398, 407

GRAVITY Collaboration (Abuter, R., et al.) 2018, A&A, 615, L15

GRAVITY Collaboration (Abuter, R., et al.) 2020, A&A, 636, L5

GRAVITY Collaboration (Abuter, R., et al.) 2022, A&A, 657, L12

GRAVITY Collaboration (Abd El Dayem, K., et al.) 2024, A&A, 692, A242

Grould, M., Vincent, F. H., Paumard, T., & Perrin, G. 2017, A&A, 608, A60

Hees, A., Do, T., Ghez, A. M., et al. 2017, Phys. Rev. Lett., 118, 211101

Hees, A., Folkner, W. M., Jacobson, R. A., & Park, R. S. 2014, Phys. Rev. D, 89, 102002

Hinterbichler, K. 2012, Rev. Mod. Phys., 84, 671

Hofmann, F., & Müller, J. 2018, Class. Quant. Grav., 35, 035015

Hoyle, C. D., Schmidt, U., Heckel, B. R., et al. 2001, Phys. Rev. Lett., 86, 1418

Iorio, L. 2024, General Post-Newtonian Orbital Effects: From Earth's Satellites to the Galactic Centre (Cambridge University Press)

Jovanović, P., Borka, D., Borka Jovanović, V., & Zakharov, A. F. 2021, Eur. Phys. J. D, 75, 145

Jovanović, P., Jovanović, V. B., Borka, D., & Zakharov, A. F. 2023, JCAP, 03, 056

Jovanović, P., Jovanović, V. B., Borka, D., & Zakharov, A. F. 2024a, Phys. Rev. D, 109, 064046

Jovanović, P., Borka Jovanović, V., Borka, D., & Zakharov, A. F. 2024b, Symmetry, 16, 397

Kiefer, C. 2023, arXiv e-prints [arXiv:2302.13047]

Konopliv, A. S., Asmar, S. W., Folkner, W. M., et al. 2011, Icarus, 211, 401

Losada, V. d. M., Della Monica, R., de Martino, I., & De Laurentis, M. 2025, A&A, 694, A280

Mariani, V., Fienga, A., Minazzoli, O., Gastineau, M., & Laskar, J. 2023, Phys. Rev. D, 108, 024047

Massey, R., Kitching, T., & Richard, J. 2010, Rept. Prog. Phys., 73, 086901

Moffat, J. W. 2006, JCAP, 03, 004

Nitz, A. H., Capano, C. D., Kumar, S., et al. 2021, ApJ, 922, 76

Peebles, P. J. E., & Ratra, B. 2003, Rev. Mod. Phys., 75, 559

Perivolaropoulos, L. 2010, Phys. Rev. D, 81, 047501

Plewa, P. M., Gillessen, S., Eisenhauer, F., et al. 2015, MNRAS, 453, 3234

Poisson, E., & Will, C. 2012, Gravity: Newtonian, Post-Newtonian, Relativistic (Cambridge: Cambridge University Press), 1

Reid, M. J., & Brunthaler, A. 2020, ApJ, 892, 39

Sabha, N., et al. 2012, A&A, 545, A70

Salucci, P. 2019, Astron. Astrophys. Rev., 27, 2

Schödel, R., Ott, T., Genzel, R., et al. 2002, Nature, 419, 694

Shankaranarayanan, S., & Johnson, J. P. 2022, Gen. Relativ. Gravitation, 54, 44

Tan, Y., & Lu, Y. 2024, Phys. Rev. D, 109, 044047

Touboul, P., et al. 2022, Phys. Rev. Lett., 129, 121102

Tsai, Y.-D., Wu, Y., Vagnozzi, S., & Visinelli, L. 2023, JCAP, 04, 031

Tsai, Y.-D., Farnocchia, D., Micheli, M., Vagnozzi, S., & Visinelli, L. 2024, Commun. Phys., 7, 311

Vagnozzi, S., et al. 2023, Class. Quant. Grav., 40, 165007

van Albada, T. S., Bahcall, J. N., Begeman, K., & Sancisi, R. 1985, ApJ, 295, 305

Visser, M. 1998, Gen. Rel. Grav., 30, 1717

Weinberg, S. 1989, Rev. Mod. Phys., 61, 1

Will, C. M. 2014, Living Rev. Rel., 17, 4

Will, C. M. 2018a, Theory and Experiment in Gravitational Physics (Cambridge University Press)

Will, C. M. 2018b, Class. Quant. Grav., 35, 17LT01

Zakharov, A. F., Jovanovic, P., Borka, D., & Jovanovic, V. B. 2016, JCAP, 05, 045

Zakharov, A. F., Jovanović, P., Borka, D., & Borka Jovanović, V. 2018, JCAP, 04, 050

<footnote n="1">1 LIRA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université de Paris, 5 place Jules Janssen, 92195 Meudon, France</footnote>

<footnote n="2">2 Max Planck Institute for Extraterrestrial Physics, Giessenbachstraße 1, 85748 Garching, Germany</footnote>

<footnote n="3">3 Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote n="4">4 European Southern Observatory, Karl-Schwarzschild-Straße 2, 85748 Garching, Germany</footnote>

<footnote n="5">5 Max Planck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany</footnote>

<footnote n="6">6 1st Institute of Physics, University of Cologne, Zülpicher Straße 77, 50937 Cologne, Germany</footnote>

<footnote n="7">7 CENTRA - Centro de Astrofísica e Gravitação, IST, Universidade de Lisboa, 1049-001 Lisboa, Portugal</footnote>

<footnote n="8">8 Universidade de Lisboa - Faculdade de Ciências, Campo Grande 1749-016, Lisboa, Portugal</footnote>

<footnote n="9">9 European Southern Observatory, Casilla 19001, Santiago 19, Chile</footnote>

<footnote n="10">10 Faculdade de Engenharia, Universidade do Porto, rua Dr. Roberto Frias, 4200-465 Porto, Portugal</footnote>

<footnote n="11">11 Departments of Physics & Astronomy, Le Conte Hall, University of California, Berkeley, CA 94720, USA</footnote>

<footnote n="12">12 Max Planck Institute for Astrophysics, Karl-Schwarzschild-Straße 1, 85748 Garching, Germany</footnote>

<footnote n="13">13 Max Planck Institute for Radio Astronomy, auf dem Hügel 69, 53121 Bonn, Germany</footnote>

<footnote n="14">14 Institute of Multidisciplinary Mathematics, Universitat Politècnica de València, València, Spain</footnote>

<footnote n="15">15 Center of Gravity, Niels Bohr Institute, Blegdamsvej 17, 2100 Copenhagen, Denmark</footnote>

<footnote n="16">16 ORIGINS Excellence Cluster, Boltzmannstraße 2, 85748 Garching, Germany</footnote>

<footnote n="17">17 Department of Physics, Technical University of Munich, 85748 Garching, Germany</footnote>

<footnote n="18">18 Higgs Centre for Theoretical Physics, Edinburgh, UK</footnote>

<footnote n="19">19 Leiden University, 2311 EZ Leiden, The Netherlands</footnote>

<footnote n="20">20 The Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, China</footnote>
# Appendix A: Details of the numerical integration

The numerical integration of the equation of motion is performed using the Python library `scipy.integrate.solve_ivp` with a Runge-Kutta 5(4) algorithm, which means that the steps are evaluated using a fifth-order method, while the error is controlled assuming the accuracy of the fourth-order method. The convergence of the integration is ensured by looking at the conservation of energy over the entire integration period (almost two orbits in ~ 30 years gives $\Delta E/E \sim O(10^{-10})$).

Kepler's equation can be solved instead using a Python root finder (`scipy CURVE FITTING`) which implements the Newton-Raphson method. The latter solves the equation with a precision of $O(10^{-16})$.

# Appendix B: Coordinates transformations and inclusion of relativistic effects.

## B.1. Coordinate transformation

The transformation from the orbital reference frame to the observer reference frame can be achieved by using the following conversion:

$$ 
\begin{aligned}
x' &= Ax_{BH} + Fy_{BH}, & v_{x'} &= Av_{xBH} + Fv_{yBH}, \\
y' &= Bx_{BH} + Gy_{BH}, & v_{y'} &= Bv_{xBH} + Gv_{yBH}, \\
z_{\text{obs}} &= -(Cx_{BH} + Hy_{BH}), & v_{z_{\text{obs}}} &= -(Cv_{xBH} + Hv_{yBH}),
\end{aligned} 
\quad (B.1) $$

where A, B, C, F, G, H are the Thiele-Innes parameters (Catanzarite 2010) defined as:

$$ 
\begin{aligned}
A &= \cos \Omega \cos \omega - \sin \Omega \sin \omega \cos i, \\
B &= \sin \Omega \cos \omega + \cos \Omega \sin \omega \cos i, \\
F &= -\cos \Omega \sin \omega - \sin \Omega \cos \omega \cos i, \\
G &= -\sin \Omega \sin \omega + \cos \Omega \cos \omega \cos i, \\
C &= -\sin \omega \sin i, \\
H &= -\cos \omega \sin i.
\end{aligned} 
\quad (B.2) $$

The Cartesian coordinates {$x_{BH}, y_{BH}, z_{BH}$} and velocities {$v_{xBH}, v_{yBH}, v_{zBH}$} are those obtained from the numerical integration. For a more detailed discussion on how the coordinate system {$x', y', z_{\text{obs}}$} and the above transformation are defined, we refer to Figure 1 and Appendix B of Grould et al. (2017).

## B.2. Relativistic effects and Rømer's delay

To produce a better fit, there are observational effects that must be included in the model.

Rømer's delay is the difference between the time of emission of the signal $t_{\text{em}}$ and the actual observational dates $t_{\text{obs}}$, due to the finite speed of light. To include this delay, we used the first order Taylor's expansion of the Røemer equation, which is expressed as:

$$ t_{\text{em}} = t_{\text{obs}} - \frac{z_{\text{obs}}(t_{\text{obs}})}{1 + v_{z_{\text{obs}}}(t_{\text{obs}})}. \quad (B.3) $$

The difference between the exact solution of Røemer equation and the approximated solution in (B.3) is at most ~ 4 s over the S2 orbit and therefore negligible. The Røemer effect affects both the astrometry and the spectroscopy, with an impact of ≈ 450 µas on positions and ≈ 50 km/s at periastron on radial velocities. Our results recover the previous estimates for this effect reported in Grould et al. (2017), GRAVITY Collaboration (2018).

Moreover, there are two relativistic effects that must be taken into account when S2 approaches the periastron: the relativistic Doppler shift and the gravitational redshift. Both induce a shift in the spectral lines of S2 that affects the radial velocity measurements. The former is given by

$$ 1 + z_D = \frac{1 + v_{z_{\text{obs}}}}{\sqrt{1 - v^2}}, \quad (B.4) $$

while the gravitational redshift is defined as

$$ 1 + z_G = \frac{1}{\sqrt{1 - 2U(r_{\text{em}})}}, \quad (B.5) $$

where $U(r_{\text{em}})$ is the potential in Eq. (1) evaluated at the time of emission $t_{\text{em}}$.

The two shifts can be combined using Eq. (D.13) of Grould et al. (2017) to obtain the total radial velocity

$$ V_R \approx \frac{1}{\sqrt{1 - \epsilon}} \cdot \frac{1 + v_{z_{\text{obs}}}}{\sqrt{1 - v^2/(1 - \epsilon)}} - 1. \quad (B.6) $$

where $\epsilon = 2U(r_{\text{em}})$.

In the total space velocity, $v = |v|$, we must also add a correction due to Solar System motion. We followed the most recent work of Reid & Brunthaler (2020) and take a proper motion of Sgr A* of

$$ 
\begin{aligned}
v_x^{\text{SSM}} &= -5.585 \text{ mas/yr} = 6.415 \cos(209.47^\circ) \text{ mas/yr}, \\
v_y^{\text{SSM}} &= -3.156 \text{ mas/yr} = 6.415 \sin(209.47^\circ) \text{ mas/yr}.
\end{aligned} 
\quad (B.7) $$