Mass distribution in the Galactic Center based on interferometric astrometry of multiple stellar orbits

Stars orbiting the compact radio source Sgr A* in the Galactic Center serve as precision probes of the gravitational ﬁeld around the closest massive black hole. In addition to adaptive optics-assisted astrometry (with NACO / VLT) and spectroscopy (with SINFONI / VLT, NIRC2 / Keck and GNIRS / Gemini) over three decades, we have obtained 30–100 µ as astrometry since 2017 with the four-telescope interferometric beam combiner GRAVITY / VLTI, capable of reaching a sensitivity of m K = 20 when combining data from one night. We present the simultaneous detection of several stars within the di ﬀ raction limit of a single telescope, illustrating the power of interferometry in the ﬁeld. The new data for the stars S2, S29, S38, and S55 yield signiﬁcant accelerations between March and July 2021, as these stars pass the pericenters of their orbits between 2018 and 2023. This allows for a high-precision determination of the gravitational potential around Sgr A*. Our data are in excellent agreement with general relativity orbits around a single central point mass, M • = 4 . 30 × 10 6 M (cid:12) , with a precision of about ± 0 . 25%. We improve the signiﬁcance of our detection of the Schwarzschild precession in the S2 orbit to 7 σ . Assuming plausible density proﬁles, the extended mass component inside the S2 apocenter ( ≈ 0 . 23 (cid:48)(cid:48) or 2 . 4 × 10 4 R S ) must be (cid:46) 3000 M (cid:12) (1 σ ), or (cid:46) 0 . 1% of M • . Adding the enclosed mass determinations from 13 stars orbiting Sgr A* at larger radii, the innermost radius at which the excess mass beyond Sgr A* is tentatively seen is r ≈ 2 . 5 (cid:48)(cid:48) ≥ 10 × the apocenter of S2. This is in full harmony with the stellar mass distribution (including stellar-mass black holes) obtained from the spatially resolved luminosity function.


Introduction
The GRAVITY instrument on the Very Large Telescope Interferometer (VLTI) has made it possible to monitor the positions of stars within 0.1" from Sgr A*, the massive black hole (MBH) at the Galactic Center (GC), with a precision of ≈ 50µas (Gravity Coll. 2017). The GRAVITY data taken in 2017-2019 together with the adaptive optics (AO) and Speckle data sets obtained since 1992 (at ESO telescopes), or since 1995 at the Keck telescopes have delivered exquisite coverage of the 16-year highly elliptical orbit of the star S2, which passed its most recent pericenter in May 2018. Besides the direct determinations of the mass of Sgr A* (M • ) and the distance to the GC (R 0 ), these interferometric data have provided strong evidence for general relativistic effects caused by the central MBH on the orbit of S2, namely, the gravitational redshift and the prograde relativistic precession (Gravity Coll. 2018b, 2019, 2021Do et al. 2019).
Due to its short period and brightness, S2 is the most prominent star in the GC, but ever higher quality, high-resolution imaging and spectroscopy of the nuclear star cluster over almost three decades have delivered orbit determinations for some 50 stars (Schödel et al. 2002Ghez et al. 2003Ghez et al. , 2008Eisenhauer et al. 2005;Gillessen et al. 2009Gillessen et al. , 2017Meyer et al. 2012;Boehle et al. 2016). The motions of these stars show that the gravitational potential is dominated by a compact, central mass of ≈ 4.3 × 10 6 M that is concentrated within S2's (3D) pericenter distance of 14 mas (or 120 AU) and 1400 times the event horizon radius R S of a Schwarzschild (non-rotating) MBH for a distance of 8.28 kpc (Gravity Coll. 2019, 2021. S2 passes its pericenter with a mildly relativistic orbital speed of 7700 km/s (β = v/c = 0.026). Based on the monitoring of the star's radial velocity and motion on the sky from data taken prior to and up to two months after pericenter, Gravity Coll. (2018b) were able to detect the first post-Newtonian effects of general relativity (GR), namely: the gravitational redshift, along with the transverse Doppler effect of special relativity. The combined effect for S2 shows up as a 200 km/s residual centered on the pericenter time, relative to the Keplerian orbit with the same parameters. Gravity Coll. (2019) improved the statistical robustness of the detection of the gravitational redshift to 20σ. Do et al. (2019) confirmed these findings from a second, independent data set mainly from the Keck telescope. While the redshift occurs solely in the wavelength space, the superior astrometry of GRAVITY sets much tighter constraints on the orbital geometry, mass, and distance, all the while decreasing the uncertainty on the redshift parameter more than three times relative to data sets from single telescopes.
The precession due to the Schwarzschild metric is predicted to lead to a prograde rotation of the orbital ellipse in its plane of ∆ω = 12.1 per revolution for S2, corresponding to a shift in the milli-arcsec regime of the orbital trace on sky; hence using interferometry is particularly advantageous in this case. Gravity Coll. (2020) detected the Schwarzschild precession at the 5σ level. The uncertainties on the amount of precession can then be interpreted as limits on how much extended mass (leading to retrograde precession) might be present within the S2 orbit.
Here, we expand our analysis by two more years, up to 2021.6. We combine GRAVITY data from four stars, alongside with the previous AO data. Section 2 presents the new data and Section 3 describes our analysis. In Section 4, we show the combined fits, improving the accuracy of the measured post-Newtonian parameters of the central black hole and the limits on the extended mass (Section 5). In combination with earlier measurements of stars with larger apocenters, we study the mass distribution out to ≈ 3". Section 6 summarizes our conclusions.

Observations
Interferometric astrometry with GRAVITY has several distinct advantages over single-telescope, AO imaging ( Figure 1): First, the higher angular resolution yields an order of magnitude better astrometric precision for isolated sources.
Secondly, for crowded environments, such as the central arcsecond that has a surface density > 100 stars per square arcsecond to K < 17 (and more for fainter limits, Genzel et al. 2003;Baumgardt et al. 2018;Waisberg et al. 2018), interferometric data are much less affected (by a factor of several hundred) by confusion noise. In the context of the GC cluster imaging, this issue was recognized early on (Ghez et al. 2003(Ghez et al. , 2008Gillessen et al. 2017;Do et al. 2019;Gravity Coll. 2020). For "orbit crossings" of modest duration for individual brighter stars, this often means that data over a duration of a few years are affected. The situation is worse at the pericenter passage of S2 (2002,2018), when the star and the variable source Sgr A* are in the same diffraction beam of an 8−10 m class telescope (Ghez et al. 2008). For 2021/2022, our data show explicitly that in addition to Sgr A*, several stars are present in the central beam (see also Gravity Coll. 2021), making single-telescope astrometry even more uncertain or unusable.
Third, close to Sgr A*, astrometry with interferometry reduces to fitting the phases and visibilities with a multiple point source model in a single pointing of the interferometric fiber, which is straightforward, once the optical aberrations across the fiber field of view are corrected for (Gravity Coll. 2020, 2021. Interferometric measurements beyond the fibre field of view require double pointings, which, thanks to GRAVITY's metrology system, can be related astrometrically to each other (Appendix A.1.2). The GRAVITY positions are directly referring to Sgr A*, since it is visible in each exposure and since it is one of the point sources in the multi-source model. In contrast, AO astrometry relies on establishing a global reference frame by means of stars visible both at radio wavelengths (where Sgr A* is visible as well) and in the near-infrared. A few such SiO maser stars exist in the GC (Reid et al. 2007) and they reside at larger separations (≈ 20 ). Hence, to establish an AO-based reference frame, it is necessary to correct for the distortions of the imaging system (Plewa et al. 2015;Sakai et al. 2019;Jia et al. 2019).

Analysis
For a single-star fit, we typically fit for 14 parameters: six parameters describing the initial osculating Kepler orbit (a, e, i, ω, Ω, t 0 ), R 0 and M • , along with five coordinates describing the on-sky position and the 3D velocity of the mass (relative to the AO spectroscopic/imaging frame), and a dimensionless parameter encoding the non-Keplerian effect we are testing for. For the gravitational redshift, we used f gr , which is 0 for Newtonian orbits and 1 for GR-orbits. In Gravity Coll. (2018b) we found f gr = 0.90 ± 0.17; and in Gravity Coll. (2019), we found f gr = 1.04 ± 0.05. Do et al. (2019) reported f gr = 0.88 ± 0.17. For the Schwarzschild precession, we use the first-order post-Newtonian expansion for a massless test particle (Will 1993) and add a factor f SP in the equation of motion in front of the 1PN terms, where f SP = 0 corresponds to Keplerian motion and f SP = 1 to GR. In Gravity Coll. (2020), we found f SP = 1.10 ± 0.19.
Similarly, we parameterize an extended mass distribution by including a parameter f Pl in the normalization of the profile. Following Gillessen et al. (2017) and Gravity Coll. (2020), we assume a Plummer (1911) profile with scale length, a Pl , and total mass of f Pl M • . We use a Pl = 1.27a apo (S2) = 0.3 (Mouawad et al. 2005). The enclosed mass within R is M(≤ R) = f Pl M • (R/a Pl ) 3 (1 + R 2 /a 2 Pl ) −3/2 . We fit for the fraction of M • that is in the extended mass, f Pl .
Following Gravity Coll. (2018b, 2019, 2020), we find the best-fit values by fitting simultaneously all parameters, including prior constraints. Throughout our study, we used an outlier robust fitting (Sect 3.2 in Gravity Coll. 2020). The inferred uncertainties are affected and partially dominated by systematics, especially when combining data from three or more measurement techniques. Our parameterization keeps correlations between f SP or f Pl and M • or R 0 small, but the two parameters of interest show some degeneracy with the coordinate system offsets.
To check the formal fit errors, we carried out a (Metropolis-Hastings) Markov-Chain Monte Carlo (MCMC) analysis. Using 100 000 realizations, we found the distributions and parameter correlations of the respective dimensionless parameter, f SP or f Pl , with the other parameters and test whether they are well described by Gaussian distributions. For more details, see Gravity Coll. (2018b, 2019, 2020, 2021) and Appendix A.
In Gillessen et al. (2009Gillessen et al. ( , 2017, we consistently found, based on the AO data, that the basic parameters describing the gravitational potential (M • and R 0 ) and the extended mass (e.g., f Pl for an assumed Plummer distribution) are best constrained by the S2 data. Including other stars only moderately improved the fitting quality and uncertainties. This is because of the superior number and quality of the S2 data compared to those of the other stars. Since the higher resolution GRAVITY astrometry became available, S2 has completely dominated our knowledge about the central potential. Another reason is that it is only for S2 that we have data at or near pericenter, which are the most sensitive component of the data to the mass distribution, as the explicit analysis in Gillessen et al. (2017) andHeißel et al. (2021) shows.
This situation changes with the data set used here. We now have GRAVITY data of four stars with comparable pericenters at 12 mas (S29), 14 mas (S2), 26 mas (S38), and 29 mas (S55). Naturally, we need to fit for 4 × 6 orbital parameters (a i , e i , i i , ω i , Ω i , t 0,i ) in addition to the NACO/SINFONI zero points (x 0 , y 0 , vx 0 , vy 0 , vz 0 ), M • and R 0 , as well as f SP and/or f Pl . However, the inclusion of near-pericenter GRAVITY data of S29, S38, and S55 lessens the parameter correlations and uncertainties ( Figure 2); this is also because the orbits are oriented almost perpendicular to each other in at least one of the Euler angles. Furthermore, the orbits of the four stars probe the precession over a wider parameter range in semi-major axis and eccentricity.
A comparison of the results obtained with different fitting schemes and codes of the consortium (of MPE, Univ. Cologne and LESIA) underlines the importance of two aspects of the data analysis: First, AO-based astrometry of stars near pericenter is subject to confusion, with the emission from Sgr A* and other neighboring stars contributing to the centroid. These data now can be discarded in favor of better-resolution interferometric data. Second, the NACO-frame zero point and drift on the one hand, and the pro-or retrograde precession on the other hand, are degenerate. To reduce this degeneracy, we can use three sources of additional information: (i) NACO astrometry of flares, (ii) the prior from the construction of the AO reference frame (Plewa et al. 2015), or (iii) the orbits of further stars that have sufficient phase coverage to constrain the zero points. In Gravity Coll. (2020), we combined (i) and (ii). Avoiding the potentially confused flare positions (i), here we use the combination of (ii) and the orbits of 13 further stars (iii). We derive x 0 = 0.57±0.15 mas, y 0 = −0.06±0.25 mas (both epoch 2010.35), vx 0 = 63±7 µas/yr, vy 0 = 33 ± 2 µas/yr, which are consistent with the earlier estimates, but with smaller uncertainties.

Schwarzschild precession for S2
Repeating the analysis of Gravity Coll. (2020) (S2 alone but with the updated zero points) and solving for the Schwarzschild precession parameter we find f SP = 0.85 ± 0.16 (χ 2 r = 1.11). This is naturally very similar to our previous results 1 , but the new data have decreased the 1σ uncertainty from ±0.19 to ±0.16.
Next, we fit with the four star (S2, S29, S38, S55) data, and find f SP = 0.997 ± 0.144, with χ 2 r = 2.17, (Figure 2). Figure 3 shows the residuals of the best fit and from the corresponding Newtonian ( f SP = 0) orbit. The combination of the nearpericenter GRAVITY data of four stars improves the constraints on the common parameters. The contributions raising χ 2 r > 1 come from the NACO data of S29, S38, and S55 covering the outer parts of their orbits more affected by confusion.
Applying MCMC analysis we find the most likely values of f SP = 0.85 ± 0.18 (only S2) and f SP = 0.99 ± 0.15 (S2, S29, S38, S55). Figure B.1 shows the full set of parameter correlations, including the well-known degeneracy between M • and R 0 (Ghez et al. 2008;Boehle et al. 2016;Gillessen et al. 2009Gillessen et al. , 2017. All of the 32 parameters of the four-star fit are well constrained. As discussed by Gravity Coll. (2020), the impact of the high eccentricity of the S2 orbit (e = 0.88) is that most of the precession happens in a short time-frame around pericenter. Due to the geometry of the orbit most of the precession shows up in the RA-coordinate, and the change in ω after pericenter appears as a kink in the RA-residuals. The data are obviously in excellent agreement with GR. Compared to Gravity Coll. (2020), the significance of this agreement has improved from 5 to 7 σ, from the combination of adding two more years of GRAVITY data to the S2 data set and the expansion to a four-star fit. Table B.1 gives the best fit orbit parameters, zero points, M • , and R 0 .
As of 2021, S2 is sufficiently far away from pericenter, such that the Schwarzschild precession can now be seen as a 1 Following strictly Gravity Coll. (2020), i.e., using the S2 data plus the flare positions of Sgr A* with the zero point priors of Plewa et al. (2015) yields f SP = 1.23 ± 0.14 (χ 2 r = 1.70).
≈ 0.6 mas shift between the data sets in RA (and less so in Dec) between two consecutive passages of the star on the apocenterside of the orbit. This effect is obvious when comparing the 2021 GRAVITY data to the 2005 NACO data, exactly one period prior (Figure 3 right). This comparison illustrates that the Schwarzschild precession dominates the entire orbit and that there is no detectable retrograde (Newtonian) precession due to an extended mass component (see Heißel et al. 2021).

Limits on extended mass
In the following, we fix f SP = 1 at its GR value and allow now for an extended mass component parameterized by f Pl . We find f Pl = (2.7±3.5)×10 −3 from a single S2 fit and f Pl = (−3.8±2.4)×10 −3 for the four-star fit, in accordance with the findings of Heißel et al. (2021). The latter 1σ error is consistent, albeit three to four times smaller than that of Gillessen et al. (2017) and 1.7 times smaller than that of Gravity Coll. (2020), corresponding to 2400 M within the apocenter of S2. In Figure 4, we include the 3σ uncertainty as a conservative upper limit, indicating that the extended mass cannot exceed 7500M . As in Gillessen et al. (2017) and Gravity Coll. (2020), we find again that varying a Pl or replacing the Plummer distribution by a suitable power law changes this result trivially. Furthermore, we get a weaker limit by a factor of ≈ 2 when omitting the NACO astrometry, using only GRAVITY & SINFONI data. The impact of an extended mass is naturally largest near apocenter of an orbit (e.g. Heißel et al. 2021). Figure C.2 shows the result of adding various amounts of extended mass on top of the best-fit residuals with a point mass only. Our data are just commensurate with an additional f Pl = 0.25% of M • , but a larger mass is excluded by both the near-peri-and near-apocenter data. The apparent sensitivity of the near-pericenter data in Figure C.2 is the result of referring the residuals to the osculating Keplerian orbit at apocenter in 2010.35, such that the accumulating retrograde precession enters the near-pericenter data.  Dec (bottom) between the GRAVITY (cyan filled circles, and 1σ uncertainties) and NACO data (blue filled circles) and the best GR fit (red curve, f SP = 1, Rømer effect, plus special relativity, plus gravitational redshift and Schwarzschild precession), relative to the same orbit for f SP = 0 (Kepler/Newton, plus Rømer effect, plus special relativity, plus redshift). The orbital elements for non-Keplerian orbits (i.e., with f SP 0) are interpreted as osculating parameters at apocenter time, 2010.35. NACO and GRAVITY data are averages of several epochs. The grey bars denote averages over all NACO data near apocenter (2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013). Top right: the same for the residual orbital angle on the sky δφ = φ( f SP = 1) − φ( f SP = 0). Bottom right: Zoom into the 2005/2021 part of the orbit, plotted in the mass rest frame. The earlier orbital trace does not coincide with the current one due to the Schwarzschild precession.
A second, independent measurement of the dynamically inferred mass distribution comes from fitting for the central mass using 13 individual stellar orbits with a = 0.1 to 3.8 (Gillessen et al. 2017), with R 0 and zero-points fixed to the best fitting values of the four-star fit (Figure C.1). We then averaged the results in four groups of four stars with 0.11 < a < 0.22 , five stars with 0.27 < a < 0.4 , three stars with 0.55 < a < 1.6 , and two stars with 1.6 < a < 3.8 . Most of the stars in the first two groups are classical "S-stars" (mostly early-type B stars) with typically large eccentricities (e.g., Ghez et al. 2008;Gillessen et al. 2017), while most of the stars in the outer groups are O and B stars in the clock-wise disk (Paumard et al. 2006;Bartko et al. 2009;Lu et al. 2009), with modest eccentricities. The stars in the first two groups indicate that the mass is consistent with M • to within 0.3-0.6%. There is no indication for an extended mass larger than ≈ 25, 000M within 2 a apo (S2) ≈ 0.5 . The outer stars suggests an extended mass of 15, 000M (and conservatively a 3σ limit of 50, 000M ) within 5 − 10 a apo (S2). Figure 4 summarizes the mass distribution within 5 (≈ 20× the apocenter of S2). These estimates and limits are in excellent agreement with the distribution of stars (and stellar mass black holes and neutron stars) contained in this inner region around Sgr A*, as estimated from models and simulations (Alexander 2017;Baumgardt et al. 2018), or from observations of faint stars and diffuse stellar light (Figure 4, Genzel et al. 2010;Gallego-Cano et al. 2018;Schödel et al. 2018;Habibi et al. 2019).
In summary, several precise (O(0.1-0.3%), 1σ) determinations show that the mass distribution in the GC within 5" ≈ 5 × 10 5 R S of Sgr A* is dominated by a central, compact mass. This mass is definitely enclosed within the pericenter of S29 (12 mas, ≈ 1200R S ). Taking the gas motions at ≈ 3 − 5 R S (Gravity Coll. 2018a) and the mm-size of Sgr A* (Doeleman et al. 2008;Johnson et al. 2017;Issaoun et al. 2019) into account, the data are in excellent agreement with the MBH paradigm.

Conclusions
Here, we present GRAVITY data obtained at the VLTI in 2021. Within the central 20 mas, we observe the motions of four stars between March and July, illustrating the power of the interferometric, high spatial resolution. Using the novel astrometry of the stars S2, S29, S38, and S55, along with new radial velocities obtained with GNIRS, we update our orbital analysis.
The star S2 has now returned to the part of its 16-year orbit for which good NACO AO-assisted positions were obtained during its previous passage. A direct comparison of the positions confirms that the orientation of the orbital ellipse has indeed shifted in its plane by the 12.1' expected from the prograde Schwarzschild precession induced by the gravitational field of the MBH, as reported in Gravity Coll. (2020).
At K = 14.1, S2 is comparably bright. With its increased distance from Sgr A* in 2021, we are now able to map with GRAVITY the immediate vicinity of the MBH to significantly fainter objects. This provided accurate positions for S29, S38, and S55. These stars have previously measured NACO positions when they were further away from Sgr A*. Combining these with the GRAVITY data improves the orbital parameters of the three stars substantially. In particular, S29 is on a deeply plunging (e = 0.97) orbit with a period of ≈ 90 years and pericenter passage in 2021.41, with a space velocity of ≈ 8740 km/s at only 100 AU from Sgr A*. S2, S29, S38, and S55 orbit in the same gravitational potential, and combining their astrometry and radial velocity data improves the accuracy of the determination of the properties of the central MBH. This leads to a 14% measurement precision of the Schwarzschild precession, which is in full agreement with the prediction of GR. The best fit further yields R 0 = (8277 ± 9) pc and M • = (4.297±0.012)×10 6 M (statistical errors, see Gravity Coll. 2021 for a discussion of the systematics that are ≈ 30 pc for R 0 and ≈ 40, 000M for M • ). Any smooth extended mass distribution would lead to a retrograde precession of the S2 orbit relative to the relativistically precessing one and we can thus place a limit on a hypothetical mass distribution. The measurement errors leave room for at most 3000 M in extended mass out to 230 mas. We included 13 stars further out with earlier measurements to trace the mass as a function of radius. The data are fully consistent with a single point mass, and only at r 2.5 does the enclosed mass tentatively exceed M • , which is consistent with the theoretically expected stellar mass distribution. Inside the 100 AU pericenter of S29, the orbits of Sgr A* flares (Gravity Coll. 2018a), together with the coincidence of mass location and light centroid (Plewa et al. 2015;Reid & Brunthaler 2020) constrain the mass distribution further, excluding, for example, dark matter spikes, as proposed by Becerra-Vergara et al. (2020, 2021, and the presence of an intermediate mass black hole as well (Gravity Coll. 2020).
Our multi-epoch GRAVITY data also confirm that at any time, there are likely a few stars that are sufficiently close to Sgr A* on the sky to systematically influence its position derived with AO-assisted imaging on single telescopes. In addition, in 2022, two stars will pass the pericenters of their orbits at less than 100 mas distance (S38 and S42). The upgrade of GRAVITY to GRAVITY+ (Gravity+ Coll. 2021) will push the sensitivity limit to K > 20, which may reveal more stars with even smaller orbits. The 39 m ELT equipped with MICADO (Davies et al. 2021) and HARMONI (Thatte et al. 2021) might be the prime choice for obtaining radial velocities of such stars. Yet, GRAV-ITY+ will beat the ELT's angular resolution by a factor three, allowing continued < 50 µas astrometry and going even deeper than what we have demonstrated so far (Gravity Coll. 2021 The full width half maximum (FWHM) of the interferometric field of view (IFOV) of GRAVITY is 70 mas. In consequence, not all stars discussed in this paper are observable simultaneously. The star S2 has moved too far away from Sgr A* compared to 2018 to be observable simultaneously with Sgr A*, while the stars S55 and S29 (and others, see Gravity Coll. 2021) are always observed alongside Sgr A*. Depending on the separation, there are two methods for determining the positions of the stars relative to Sgr A*: single-beam and dual-beam astrometry. Single-beam positions are extracted from pointings in which more than one source is present in the IFOV. The distances between the stars are extracted by fitting a multi-source model to the visibility amplitude and closure phase, each of which are measured at ≈ ten spectral channels for six baselines. This yields positions of the sources with respect to each other. Since Sgr A* is visible in all our central frames, for those pointings the relative positions also are the absolute ones, that is, with respect to the mass center. If the stars are not observable in a single IFOV, we need to observe them separately and apply the dual-beam technique. For the case of two isolated stars, one interferometrically calibrates the first source with the second. The first source serves as a phase reference relative to which offsets of the second source can be measured.
Appendix A.1.1: Single-beam astrometry If a star is in the same IFOV as Sgr A* (of particular interest here in 2021 are S29 and S55), we determine the relative separation between the star and Sgr A* by interferometric model fitting to the visibility amplitude and closure phases in the Sgr A* pointing. This methodology is unchanged with regard to the way separations were determined in Gravity Coll. (2021). We thus take into account the effect of phase aberrations as well as bandwidth smearing (Gravity Coll. 2020).
Appendix A.1.2: Dual-beam astrometry If the stars are separated by more than the IFOV of GRAVITY, we measure the separation between the two sources by using one of them as the phase reference for the other target, namely, by calibrating the complex visibilities of the object of interest with the ones of the phase reference. We use S2 as the phase reference, which in its interferometric observables is consistent with a single point source. The separation between any star and S2 is determined by two vectors: the vector by which the IFOV was moved between S2 and the star. This vector is measured by the metrology system of GRAVITY monitoring the internal optical path differences the phase center offset in the S2-calibrated star observation.
It is determined by fitting the visibility phase. Also for the dual-beam analysis, we account for the effect of phase aberrations and bandwidth smearing when calculating the model visibility phase.
The separation is affected by inaccuracies and systematic uncertainties of the metrology. Such telescope-based errors are inherent to the dual-beam part of the measurement. Typically, we find more stars than just one in the IFOV. We thus need to take into account the signatures induced by the ad-ditional stars in the dual-beam measurement. This occurs, for example, for the Sgr A* pointings (where S29, S55, S62, and S300 are present), but also for the S38 pointing (with S60 and S63 being present in the IFOV). We thus fit the interstellar separations and the phase center offset simultaneously in order to take into account their degeneracies. However, the separation vectors are mostly sensitive to the visibility amplitude and the closure phase information from the visibility phase, while the phase center offsets mostly acts as an additional term in the visibility phase. In this way we can relate all positions to our calibrator source, S2. Hence, we can relate the positions also to Sgr A* by subtracting the star-to-S2 and S2-to-Sgr A* separations.
Telescope-based errors cancel out in the closure phase and therefore the relative positions of the sources are not affected by phase errors, but the visibility phases carry the information of how different IFOVs are located with respect to one another. We find that by fitting the closure phases and the visibility phases with equal weights, we minimize the effect of the telescopebased errors, while still being sensitive to the phase information.
In order to average out the phase errors, we calibrate all N frames of a given pointing with all M available S2 frames individually. For each of the N × M resulting data sets, we determine the phase center position and average the resulting phase center locations. This calibration uncertainty adds a systematic uncertainty of 60 µas, divided by the square root of the number of available calibrations.
We further improve the accuracy of our phase center measurement by determining the best fit fringe-tracker and science target separation by fitting the S2 observations with a drifting point source model. This takes into account our imperfect knowledge of the separation prior to the observation. Here, we follow the concepts set out in Gravity Coll. (2020).

Appendix A.2: GRAVITY: Deep imaging
To obtain deep, high-resolution images of the GC, we developed a new imaging code called GRAVITY-RESOLVE or G R which is drawn from RESOLVE (Arras et al. 2021), a Bayesian imaging algorithm formulated in the framework of information field theory (Enßlin 2019) that is custom-tailored to GRAVITY observations of the GC. Here, we briefly outline the main ideas, (see Gravity Coll. (2021) for a detailed description.
With a Bayesian forward modeling approach, we can address data sparsity and account for various instrumental effects that render the relation between image and measurement more complicated than the simple Fourier transform of the van-Cittert Zernike theorem. To this end, the algorithm formulates a prior model which permits to draw random samples, processes them with the instrumental response function, and evaluates the likelihood to compare the predicted visibilities with the actual measurement. This approach can handle the non-invertible measurement equation and allows us to work with non-linear quantities such as closure phases. The exploration of the posterior distribution is done with metric Gaussian variational inference (Knollmüller & Enßlin 2019), and infers the mean image with respect to the posterior jointly with an uncertainty estimate. There are already some imaging tools available for optical/near-IR interferometry that implement a forward-modeling approach such as MIRA (Thiébaut 2008) or SQUEEZE (Baron et al. 2010). Our code differs from those with regard to the details of the measurement equation, the prior model, and how the maximization and exploration of the posterior are performed.
In the measurement equation, we implemented all instrumental effects relevant for GRAVITY: coupling efficiency, aber-ration corrections (Gravity Coll. 2021), averaging over finite sized wavelength channels (also known as bandwidth smearing), and the practice in optical and near-IR interferometry to construct the complex visibility as the coherent flux over a baseline divided by the total flux of each of the two telescopes. The latter signifies that the visibility amplitude can be unity at most, but coherence loss can degrade the observed visibility from the theoretical expectation. This we account for by a self-calibration approach where we infer a time-and baseline-dependent calibration factor jointly with the image. An appropriate prior model is essential to address the data sparsity inherent to optical/near-IR interferometry, and we specifically tailor it to the GC observations. There, we see Sgr A* as a point source in addition to some relatively bright stars whose approximate positions are known from orbit predictions. For those objects, we directly infer the position and brightness using a Gaussian and a log-normal prior respectively. The variability and polarization of Sgr A* is accounted for by allowing for an independent flux value in each frame and polarization state observed. In the actual image itself, we expect to see few faint, yet unknown, point sources and thus impose the individual pixels to be independent with their brightness following an Inverse Gamma distribution. We note that all sources other than Sgr A*, that is all non-variable sources, could in principle also be attributed to the image. However, modeling them as additional point sources improves convergence and mitigates pixelization errors.

Appendix A.3: GNIRS: Determining radial velocities
In 2021 we had four successful observations with the long-slit spectrograph GNIRS using the AO system ALTAIR at the Gemini observatory. We used the long slit in the K-Band with the 10.44 l/mm grating. The slit was positioned such that we observed S2 and S29 simultaneously. To calibrate the data we used the daytime calibration from the day after the observation, which contains a set of dark frames to determine a bad pixel mask, flat frames, and a wavelength calibration. Additionally, a telluric star was observed right after the observation. To determine the velocity of the stars we used template fitting with a high SNR S2 spectrum, in the same way as we extracted the SINFONI velocities (see Gravity Coll. 2018b). We were able to detect a velocity for S2 in all four observing nights. As S29 is significantly fainter than S2 we needed excellent conditions to get a detection, which was only possible in two of the four nights.

Appendix B: Fit details
In Table B.1, we give the best-fit parameters of the four-star fitdetermining f SP , comparing also with similar fits from the literature. Figure B.1 gives the full posterior of the four-star fit in the form of a corner plot.

Appendix C: Additional figures
In Figure C.1 we show the orbital data of additional S-stars that were auxiliary in this work. Figure C.2 illustrates that our S2 data are compatible at most with an extended mass component of around 0.1% enclosed within the S2-orbit.

Appendix D: List of observations
In Tables D.1, D.2, and D.3 we list the observations used in this work. We note that most of the raw data are available from the observatories' archives, most easily found via the program ID numbers provided here. The derived astrometric and spectroscopic data set can be made available upon personal request and after discussing the terms of a collaboration on a case-by-case basis. The line M • [10 6 M ] 8277 pc gives the masses rescaled to a common distance of R 0 = 8277 pc, using M ∝ R 2 0 (Gillessen et al. 2017). the offsets x 0 , y 0 , the distance R 0 , the velocity offsets vx 0 , vy 0 , vz 0 , the precession parameter f SP , followed by four times six orbital elements for each of the four stars used, in the order a, e, i, Ω, ω, t peri for S2, S29, S38, S55. Orbital data for S9, S13, S18 and S21. Middle: S4, S12, S31 and S42. Right: S66, S67, S87, S96 and S97. These data are complemented by multi-epoch spectroscopy for the orbital fitting.  Figure 3, with the dashed red curve denoting the f SP = 1 GR curve for the best fitting orbit and mass. In addition, we show orbital models with the same central mass, distance and orbital parameters but now adding an extended mass component assumed to have a Plummer shape (Gillessen et al. 2017) showing the impact of adding a Plummer mass of M ext within the 0.25 apocenter radius of S2. Black, orange and blue solid curves show the changes expected if this extended Plummer mass is 0.1, 0.3 and 0.6% of M • (4.4 × 10 3 , 8.9 × 10 3 and 1.78 × 10 4 M within the apocenter of S2, R apo = 0.24 ). Formal fitting shows that such an extended mass greater than about ≈ 0.1% of M • is incompatible with the data.     S2 2003S2 -05-09 Eisenhauer et al. (2003 VLT, NACO AO spectr. rad. vel. S2 2003-06-12 Eisenhauer et al. (2003 VLT