Issue 
A&A
Volume 647, March 2021



Article Number  A59  
Number of page(s)  17  
Section  Astronomical instrumentation  
DOI  https://doi.org/10.1051/00046361/202040208  
Published online  09 March 2021 
Improved GRAVITY astrometric accuracy from modeling optical aberrations^{⋆}
^{1}
Max Planck Institute for Extraterrestrial Physics, Giessenbachstraße 1, 85748 Garching, Germany
^{2}
LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université de Paris, 5 place Jules Janssen, 92195 Meudon, France
^{3}
Max Planck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany
^{4}
1st Institute of Physics, University of Cologne, Zülpicher Straße 77, 50937 Cologne, Germany
^{5}
Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France
^{6}
Universidade de Lisboa – Faculdade de Ciências, Campo Grande, 1749016 Lisboa, Portugal
^{7}
Faculdade de Engenharia, Universidade do Porto, rua Dr. Roberto Frias, 4200465 Porto, Portugal
^{8}
European Southern Observatory, KarlSchwarzschildStraße 2, 85748 Garching, Germany
^{9}
European Southern Observatory, Casilla 19001, Santiago 19, Chile
^{10}
Sterrewacht Leiden, Leiden University, Postbus 9513, 2300 RA Leiden, The Netherlands
^{11}
Departments of Physics and Astronomy, Le Conte Hall, University of California, Berkeley, CA 94720, USA
^{12}
CENTRA – Centro de Astrofísica e Gravitação, IST, Universidade de Lisboa, 1049001 Lisboa, Portugal
^{13}
Department of Astrophysical & Planetary Sciences, JILA, Duane Physics Bldg., 2000 Colorado Ave, University of Colorado, Boulder, CO 80309, USA
^{14}
Department of Particle Physics & Astrophysics, Weizmann Institute of Science, Rehovot 76100, Israel
^{15}
Institute of Astronomy, Madingley Road, Cambridge CB3 0HA, UK
^{16}
Department of Physics, Technical University Munich, JamesFranckStraße 1, 85748 Garching, Germany
^{17}
Max Planck Institute for Radio Astronomy, Auf dem Hügel 69, 53121 Bonn, Germany
Received:
22
December
2020
Accepted:
27
January
2021
The GRAVITY instrument on the ESO VLTI pioneers the field of highprecision nearinfrared interferometry by providing astrometry at the 10−100 μas level. Measurements at this high precision crucially depend on the control of systematic effects. We investigate how aberrations introduced by small optical imperfections along the path from the telescope to the detector affect the astrometry. We develop an analytical model that describes the effect of these aberrations on the measurement of complex visibilities. Our formalism accounts for pupilplane and focalplane aberrations, as well as for the interplay between static and turbulent aberrations, and it successfully reproduces calibration measurements of a binary star. The Galactic Center observations with GRAVITY in 2017 and 2018, when both Sgr A* and the star S2 were targeted in a single fiber pointing, are affected by these aberrations at a level lower than 0.5 mas. Removal of these effects brings the measurement in harmony with the dualbeam observations of 2019 and 2020, which are not affected by these aberrations. This also resolves the small systematic discrepancies between the derived distance R_{0} to the Galactic Center that were reported previously.
Key words: Galaxy: center / Galaxy: fundamental parameters / instrumentation: interferometers / instrumentation: high angular resolution / methods: data analysis
GRAVITY is developed in a collaboration by the Max Planck Institute for extraterrestrial Physics, LESIA of Observatoire de Paris/Université PSL/CNRS/Sorbonne Université/Université de Paris and IPAG of Université Grenoble Alpes/CNRS, the Max Planck Institute for Astronomy, the University of Cologne, the CENTRA – Centro de Astrofisica e Gravitação, and the European Southern Observatory.
© GRAVITY Collaboration 2021
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Open Access funding provided by Max Planck Society.
1. Introduction
The distance to the Galactic Center (GC), R_{0}, can be measured directly from stellar orbits around Sgr A*, the radio source associated with the GC massive black hole (MBH) (see, e.g., Genzel et al. 2010 and BlandHawthorn & Gerhard 2016 for a recent overview of alternative methods). To this end, the stellar proper motion, given in angle per unit time, is compared to its radial velocity, obtained in absolute length per unit time from spectroscopic observations. The GC distance then directly follows as a scaling parameter between the two measurements. Best suited for measuring R_{0} is S2, a massive young mainsequence B star on a 16year orbit with a semimajor axis of a ≃ 125 mas and an apparent Kband magnitude m_{k} ≃ 14 (Ghez et al. 2003; Eisenhauer et al. 2005; Martins et al. 2008; Gillessen et al. 2009a, 2017; Habibi et al. 2017). During its pericenter passage in 2018, S2 was closely monitored in astrometry and spectroscopy (Gravity Collaboration 2018; Do et al. 2019). In particular, the GRAVITY instrument (Gravity Collaboration 2017) directly measured the distance between S2 and Sgr A* during the flyby at a high angular resolution of around 30 μas. The combination of ultrahigh astrometric precision from nearinfrared interferometry and the spectroscopic precision of ≲10 km s^{−1} allowed us to determine the GC distance at the unprecedented precision of < 1% (Gravity Collaboration 2019).
Operating in the Kband, GRAVITY combines the light from either the four Unit Telescopes (UTs) or Auxiliary Telescopes (AT) of the ESO Very Large Telescope Interferometer (VLTI). Fringetracking on a bright reference object enables minutelong integration times on the fainter science target and measuring differential complex visibilities. The extremely high angular resolution of ≃3 mas results in very accurate astrometry with error bars between 10 μas and 100 μas (Gravity Collaboration 2017). However, the latest R_{0} measurement in Gravity Collaboration (2020) indicates a possible systematic difference with earlier determinations (Gravity Collaboration 2018, 2019). While the shift is small, of only 𝒪(1%), it is nevertheless significant because of the high precision of the measurement.
The difference in the measured GC distance coincides with a change in the observing mode. GRAVITY observes the GC with two different methods, depending on the separation between Sgr A* and S2. Close to pericenter passage, that is, in 2017 and 2018, the sources were detected simultaneously in a single fiber pointing in the socalled singlebeam mode. In later epochs, their separation exceeded the fiber field of view (FOV), and S2 and Sgr A* were targeted individually. This is referred to as dualbeam mode.
In singlebeam mode, it is not possible to align the two sources with the fiber center. To further improve the GRAVITY astrometry, we therefore conducted an analysis of how optical aberrations affect the visibility measurement across the full FOV. A similar concept of fielddependent errors already exists in radio interferometry, where it is known as directiondependent effects (DDEs; see, e.g., Bhatnagar et al. 2008; Smirnov 2011; Smirnov & Tasse 2015; Tasse et al. 2018 and references therein). The DDEs can arise either at the instrument level from the antenna beam pattern or at the atmospheric level, for example, from the ionosphere. In particular for the latest generation of interferometers (e.g. VLA, MeerKAT, and LOFAR^{1}) with a wide FOV and a large fractional bandwidth, DDEs cannot be neglected. However, to our knowledge, there is no equivalent discussion in the context of optical/nearIR interferometry.
Our analysis indeed shows that small optical imperfections in the beam combiner induce fielddependent phase errors that are reflected in the inferred binary separation. We developed an analytical model to describe this effect, and verified it by application to a dedicated testcase observation. When it is applied to the GC observations, the model induces a shift in the S2 relative position of about 0.1−0.2 mas in 2018 and ∼0.5 mas in 2017 in both right ascension (RA) and declination (Dec). Despite being small, the change is nonnegligible at the high astrometric accuracy achieved by GRAVITY. We can show that the corrected 2017 and 2018 data agree with the dualbeam observations of 2019 and 2020. Furthermore, when the correction is retroactively applied to the data sets used in Gravity Collaboration (2018, 2019), the ensuing GC distance is fully consistent with the latest result (Gravity Collaboration 2020).
We introduce the analytical model in Sect. 2 and compare it to calibration measurements in Sect. 3. Verification from the binary test case and the improved S2 position are presented in Sect. 4, and we discuss the implications for the GC distance in Sect. 5. Finally, we conclude in Sect. 6.
2. Formal description of static aberrations and their effect on visibility measurements
Static aberrations along the optical path of the instrument affect the measured visibilities by introducing a complex fielddependent factor for each telescope. We express this gain in its polar representation and decompose it into a phase map ϕ_{i}(α) and an amplitude map A_{i}(α). Here, the index i labels the telescope, and α denotes positions in the image plane. Phase and amplitudemaps lead to a modification of the observed complex visibilities V^{obs} from the wellknown van CittertZernike theorem (cf. Eq. (23)). As we demonstrate in the following, they are given by
where b_{i, j} is the baseline vector between the two telescopes, and O(α) denotes the intensity distribution of the observed astronomical object.
In this section, we show that the phase and amplitudemaps follow from optical aberrations. To this end, we start from the overlap integral, which determines the electromagnetic field from a single telescope arriving at the beam combiner. Subsequently, we propagate the effect of static aberrations from the overlap integral to the measured complex visibility to arrive at a rigorous derivation of Eq. (1). Finally, we account for the superposition of static and turbulent aberrations to obtain a formalism that is applicable in realistic observation scenarios.
2.1. Static fielddependent aberrations at fiber injection
Singlemode fibers transport the light collected by each telescope E_{tel} to the beamcombiner instrument. The overlap integral between light and the fiber mode E_{fib} then determines the transmitted electric field (Neumann 1988),
We assumed a normalized fiber mode ∫dξE_{fib}(ξ)^{2} = 1 and expressed imageplane positions by twodimensional vectors, ξ and β. Following the description of Perrin & Woillez (2019), the overlap integral was converted into the pupil plane by the ParsevalPlancharel theorem,
where ℱ^{−1} denotes the inverse Fourier transform, that is, transformation from the image to the pupil plane, and E_{obj} is the light emitted by the astronomical object. The latter is connected to ℱ^{−1}(E_{tel}) by multiplication with the pupil function P(u), corresponding to a convolution in the image plane. In the most simple case of a single point source located at α_{0}, the light is described by a pure phase . The pupil and imageplane coordinates, ξ and u, respectively, are Fourierconjugate to each other and chosen to be dimensionless. That is, any length scale in the pupil plane is given by λu, where λ refers to the wavelength and u = u. To discuss our results, we converted the dimensionless image plane coordinates ξ into the corresponding angular separation in UT observations. In an aberrationfree scenario, the pupil function of a spherical telescope with diameter 2r_{tel} and central obscuration 2r_{cent} simply is
Optical aberrations multiply the pupil function by a positiondependent complex phase, and we considered the case of purely static aberrations. These are characterized by an optical path difference (OPD) d_{pup}(u) in the pupil plane that can be expanded in terms of Zernike polynomials ,
We adopted the convention that is dimensionless and the coefficient corresponds to the root mean square of the term over the unit circle. When the turbulencefree complex fiber mode apodized by the pupil function is defined as
the overlap integral reads
The overlap integral obviously depends on the fiber profile, which for a perfectly aligned ideal singlemode fiber is
GRAVITY was designed for optimal fiber injection (Pfuhl et al. 2014), which is obtained for (Wallner et al. 2002). Here, the parameter ϵ is of order unity and describes the pupil shape.
From the comparison of the model predictions and the calibration measurements in Sect. 3.2, we find that pupilplane distortions alone are not sufficient to describe the observed aberration pattern. We also need to account for optical errors in the focal plane. Misalignment of the optical fiber as well as higher order aberrations at fiber injection introduce a complex phase in Eq. (8) and can distort the amplitude of the fiber profile.
To illustrate the effect of focal plane aberrations, we can consider the three types of misalignment depicted in Fig. 1: (A) Lateral misplacement of the fiber by (δx,δy), which in the pupil plane produces a phase slope ξ_{fib} = (δx/f, δy/f), where f is the focal length. (B) Fiber tilt by an angle φ_{fib} = (φ_{1}, φ_{2}) with respect to the optical axis of the system, which shifts the backpropagated fiber mode by u_{fib} = φ ⋅ f/λ. Finally (C), a defocus or axial fiber misplacement by δz that introduces an additional phase curvature exp[π iδzλ/f^{2} u^{2}]. When all three effects are taken into account, the generalized fiber profile, projected onto the pupil, is (Wallner et al. 2002)
Fig. 1. Schematic depiction of the pupil and focal plane aberrations that enter the overlap integral. Both effects in combination are required to describe the aberration patterns observed in calibration measurements. The lowest order aberrations in the pupil function are shown explicitly, which are (A) lateral fiber misplacement, (B) fiber tilt, and (C) defocus. Their effect is further explained in the text. 
By rearranging the phase term in the pupil plane, we can decompose it into a piston, tiptilt, and defocus
The phase terms in Eqs. (10)–(12) thus affect the overlap integral in the same way as the lowest order aberrations in d_{pup}(u). For the coordinate shift of the Gaussian profile, on the other hand, there is no such correspondence, and it alters the way in which the optical fiber scans the pupilplane aberrations.
During GRAVITY observations, the misplacement term (A) depends on the performance of the fiber tracker, but also on the uncertainty of the source position. In particular for exoplanet observations, the latter can be sizable. Fiber tilt (B) is controlled by the GRAVITY pupil tracker, and the adaptive optics calibration is one example that affects the defocus (C).
While lateral misplacement (A) and defocus (C) describe the misplacement of a pointlike fiber entrance, fiber tilt (B) accounts for the alignment of the fiber surface. This surface can exhibit irregularities beyond a simple tilt, which lead to a positiondependent OPD in the focal plane, d_{foc}(x), as illustrated in Fig. 1. Generally, aberrations from optical elements not conjugated to the pupil are fielddependent and known as Seidel aberrations. In this context, d_{foc}(x) arising in the focal plane constitutes an extreme example. It is still possible to decompose the focal plane distortions into a series of Zernike polynomials, in analogy to Eq. (5). In this representation, axial fiber offset (C) and fiber tilt (B) simply correspond to the lowest order coefficients, and higher order terms amount to a generalization of Wallner et al. (2002). Again, the phase terms introduced in by higher order aberrations are degenerate with d_{pup}(u), but the amplitude distortions need to be modeled explicitly by themselves.
Finally, for a single point source located at α_{0} in the image plane, the overlap integral averaged over a timescale much longer than the source coherence time ⟨...⟩_{obj} is
Evaluation of the Fourier transform as a function of α_{0} results in a twodimensional complex map. We show several examples of such maps in Fig. 2, assuming different Zernike coefficients to determine d_{pup}(u). The perfect Airy pattern, obtained in the limit of zero aberrations, exhibits zero phase in the central part and a phase jump by 180° at α ≃ 1.22 λ/(2r_{tel}). Antisymmetric terms, such as tilt, coma, and trefoil (not shown), only alter the location and shape of the phase jump, while defocus (not shown), astigmatism, and higher order terms produce smooth phase gradients. For a general choice of d_{pup}(u) and in the absence of focalplane aberrations, there is a saddle point where the phase maps average to zero, but significant phase shifts are encountered at larger radii.
Fig. 2. Example phase screens (top) and amplitude maps (bottom) in the image plane induced by loworder Zernike aberrations in the pupil plane at a wavelength of λ_{0} = 2.2 μm. From left to right: the considered aberrations are a perfect Airy pattern, a vertical tilt of 0.4 μm RMS, a vertical astigmatism of 0.2 μm RMS, a vertical coma of 0.2 μm RMS, and a combination of astigmatism, coma, and trefoil (with RMS 0.2 μm, 0.2 μm, and 0.1 μm, respectively). The rightmost panel also considers an additional fiber tilt with 0.2 μm RMS. 
Focalplane aberrations break the radial symmetry of the fiber profile. If the perturbations are small enough, however, the phase maps show a saddle point, but its value differs from zero and its location may be shifted. In any case, the transmitted amplitude is deformed and/or misplaced from the perfect Airy case. Pupilplane aberrations typically widen the amplitude, while imageplane aberrations have the opposite effect. They lead to a widening of the fiber in the pupil plane, and correspondingly, to a narrower imageplane profile. The exact scaling relation for the position of the Airy ring remains true only approximately in the presence of higher order aberrations such that maps at two different wavelengths, λ_{1} and λ_{2}, can be related by
2.2. Effect on visibility measurements and astrometry
The overlap integral defines the electromagnetic wave transmitted to the beam combiner from each of the four telescopes. After pairwise beam combination, the complex visibilities are obtained from the inference pattern I_{i, j},
where i and j denote the telescopes involved in the measurement, and I is the intensity. The complex pupil function enters each of these terms. Focusing on the singletelescope component first, we find from Eq. (7)
where the ⊗operator denotes autocorrelation, and O(α) = E_{obj}(α)^{2} is the brightness distribution of the observed astronomical object, which obeys
Similarly, the inference term is given by
where b_{i, j} is the baseline vector.
All optical aberrations discussed previously are encoded in the backprojected apodized pupil, which is a complex fielddependent function. Expressing the pupil function in its polar representation,
we refer to A_{i} as the telescopedependent amplitude map and to ϕ_{i} as the phase map. These quantities are closely related to the photometric and the interferometric lobes, and
respectively.
From the measured inference pattern, the complex visibilities are obtained as
By contrast, in an ideal aberrationfree setting, the vanCittertZernike theorem relates the complex visibilities to the object brightness distribution
Comparison of Eqs. (22) and (23) readily suggests that static aberrations at fiber injection distort both the measured visibility phases and amplitudes. We therefore need to adapt the interferometric equation accordingly. To make this effect even more explicit, we first considered the case of a single unresolved object at position α_{0},
In the aberrationfree case, the phase and amplitudemaps of either telescope are given by the perfect Airy pattern, shown in the leftmost panel of Fig. 2, and ϕ_{i/j}(α_{0}) equals zero or 2π. The presence of static aberrations introduces a phase shift by ϕ_{i}(α_{0}) − ϕ_{j}(α_{0}). For an interferometric binary with positions α_{1}, α_{2} and flux ratio f^{bin}, the measured visibility becomes
Finally, for a generic extended object with an intensity distribution O(α), the vanCittertZernike theorem generalizes to the expression stated at the beginning of this section in Eq. (1),
Singlepoint sources typically are observed at the fiber center, where fiber injection is highest and the phase distortions are close to zero. In situations where a very precise alignment is not possible, for example, in exoplanet observations, the visibilities can pick up some small contribution from the phase maps. For binaries with a separation comparable to the fiber width, a configuration in which the phase and amplitudemaps are irrelevant cannot be obtained in principle. In this case, the effect of static aberrations needs to be modeled and corrected for in the data analysis.
2.3. Interplay with turbulent aberrations
We did so far not consider the effect of timevarying phase aberrations. These are introduced by atmospheric turbulence or timevarying imperfections in the optical system such as tiptilt jitter from the adaptive optics. Their effect is to multiply the static pupil function by another timedependent phase,
To determine how timedependent aberrations affect the visibility measurement, we briefly recall the arguments of Perrin & Woillez (2019). When we assume that the detector integration time exceeds the coherence time of phase fluctuations by far, the longtime average ⟨…⟩_{turb} over the telescope lobes is
where D_{ϕ}(u) is the structure function of the turbulent phase (Roddier 1981), which saturates to 2σ_{ϕ} on large scales. Two assumptions underlie these expressions. First, that the fluctuations are stationary, and second, that the baseline between the telescopes is long enough for the respective apertures to become uncorrelated. As in Perrin & Woillez (2019), we assumed both to be fulfilled.
In the case of GRAVITY observations, atmospheric phase variations across the telescope apertures are corrected for by the adaptive optics system, and the turbulent aberrations are dominated by tiptilt jitter. Thus, the turbulent phase is
where the two directions of t_{i}(t) are independent and follow a Gaussian distribution with zero mean and variance . The structure function then becomes D_{t}(u) = (2πσ_{t}u)^{2}, and the photometric lobe is given by
where ⊛ denotes convolution. In case of the interferometric lobe, we further assume that the jitter is uncorrelated between telescopes, which yields
These turbulent lobes replace the static expressions of the previous sections in the prediction of the observed visibility, that is, in Eqs. (1), (24), and (25). The tiptilt jitter acts like a Gaussian convolution kernel on the static maps, which is applied to the amplitude map squared in case of the photometric lobe, but to the full complex map in the case of the interferometric lobe.
3. Measurement and characterization of aberrations for the GRAVITY beam combiner
GRAVITY observes the GC in its socalled dualfield mode, which requires the presence of a bright reference target (IRS 16C) within 2″ of the actual science targets, Sgr A* and S2. The field at each telescope is split, and reference and science source are separately injected into the fringetracking (FT) and sciencechannel (SC) fibers. Short detector integration times on the FT allow for the optical path delay to be constantly adjusted for atmospheric turbulence in order to maintain a high fringe contrast. The science channel then measures a differential visibility phase with respect to the fringe tracker on each baseline.
Phase and amplitudemaps are inherently singlefield effects in the sense that they individually affect the fringe tracker and the science channel for each telescope separately. Based on the optical layout of the fiber coupler (Pfuhl et al. 2014), there is no reason to expect equal aberrations on the SC and FT. However, the fringetracking object is a bright unresolved source that is actively tracked by the fiber center in closed loop, such that the phase distortions introduced from static aberrations are small. Moreover, any possible phase distortion from the fringe tracker cancels in the analysis of closure phases or induces a global shift without affecting the binary separation in the analysis of visibility phases. However, a description of the SC phase and amplitudemaps is essential for robustly measuring a binary separation in the science channel.
Here we report on measurements with the GRAVITY Calibration Unit (Blind et al. 2014) and on our subsequent analysis to extract SC phase and amplitudemaps. We then fitted the staticaberration model from Sect. 2.1 to these maps in order to demonstrate its validity and to obtain a compressed representation of the aberrations in form of a small number of Zernike coefficients.
3.1. Phase map measurements with the calibration unit
The GRAVITY Calibration Unit, which we use for the measurement of static aberrations, is directly attached to the beam combiner and creates the light of an artificial science and fringetracker star. By modulating the voltage on the GRAVITY positioning mirror, the position of that star relative to the fiber can be changed. We scan the FOV out to ∼70 mas in a pattern of in and outspiral, which is applied simultaneously to the FT and SC on one single telescope at a time, see Fig. 3.
Fig. 3. Examples of the scanning pattern applied in the Calibration Unit measurements. SC aberration maps were obtained with a slow modulation frequency (left). For the corresponding FT measurement, a faster scanning was used, and the right panel only shows a single iteration of in and outspiral. 
In normal observation mode, GRAVITY controls the differential OPD between science channel and fringe tracker by its laser metrology and the common path to the telescopes by fringe tracking. During the phase map calibration measurement, however, fringetracking is not possible because the fringes are lost at the margins of the scanning region. Instead, the common path from the telescope to the instrument drifts in time. Thus the determination of the aberration pattern from the absolute FT and SC phase requires a drift correction. On the FT, the short detector integration time with maximum sampling frequency of 1 kHz allows us to resolve fast modulation of the source position, and the full FOV can be scanned within ∼15 s. Over this short time span, the drift is well described by a constant velocity, which we fitted and subtracted from the data. On the SC, in contrast, the minimum detector integration time is 0.13 s and a full scan of the FOV takes 2−3 min, which is too long to model the drift by a simple polynomial fit. Instead, we obtained the science channel aberrations via a detour and first analyzed the differential driftfree SCFT phase. The pure science channel aberrations then follow from knowledge of the absolute fringe tracker phase.
The data were reduced by the standard GRAVITY pipeline, and we obtained the correlated flux in six FT spectral channels (ranging from 1.99−2.38 μm) and in medium resolution for the SC (233 wavelength bins in the range 1.97−2.48 μm). With the chosen setup, where the source position is varied on only one of the two beams forming a baseline, the measured correlated flux is given by
Thus, the measurement directly scanned the phase and amplitudemaps on the modulated channel. Potential offsets in the accompanying nonmodulated beam, ϕ_{j}(0) ≠ 0, can only cause a global phase shift, which we fitted and removed in the subsequent analysis. Finally, we considered the amplitude maps normalized to their maximum value, such that A_{j}(0) has no effect on our result.
In summary, we applied the following analysis steps to obtain the FT and the differential SCFT phase and amplitudemaps:

We fitted and subtracted a linear time drift from the phases measured in each spectral channel and on each baseline.

Phases and amplitudes were binned on a spatial grid with resolution 1 mas and averaged over all available periods of in and outspiral.

The image plane coordinates do not align perfectly with the amplitude maximum, that is, the source position for which the coupling to the fiber is most efficient. We corrected for this effect by fitting a Gaussian profile and shifting the coordinate origin to its maximum.

Interpolation over the gridded data provides the phase and amplitudemap per spectral channel and baseline.

All spectral channels were combined into a single map at reference wavelength λ_{0} = 2.2 μm by applying the approximate coordinate scaling from Eq. (14). We verified that the individual maps are consistent over the full spectral range. Crossvalidation with simulated maps showed that the error introduced by the approximate scaling relation is small, except at the very margins of the map. It furthermore cancels out between channels above and below λ_{0} to a very good degree.

From consideration of all baselines, three maps are available for each telescope. We again verified their consistency and averaged them into a single phase and amplitudemap.
This method results in an FT and a differential SCFT map for each telescope. Subtracting the former from the latter, we finally arrived at the desired SC phase map, which is shown in Fig. 4. The amplitude map on the SC, on the other hand, was measured directly.
Fig. 4. Science channel phase maps reconstructed by the procedure of Sect. 3.1 from the Calibration Unit measurement on 3 March 2020 for all four GRAVITY beams. 
The Calibration Unit measurement was performed twice with a fourmonth break, in late 2019 and early 2020, and we use the data to construct two independent sets of maps. These agree very well in the qualitative features and structures they display. On the quantitative level, the maps display moderate differences of ∼10°, which are smaller at the center and increase toward the map margins.
3.2. Representation in the pupil plane
Analyzing the Calibration Unit measurement as described in the previous subsection, we obtain the phase and amplitudemaps on a grid discretizing the image plane. We used this result to infer the underlying pupilplane and fiber aberrations, d_{pup}(u) and d_{foc}(u), in their Zernike representation. To this end, we developed a simulation tool that creates complex maps of imageplane distortions from a set of Zernike coefficients according to Eqs. (5), (6), and (13).
For the fit we considered the two Calibration Unit measurements from 2019 and 2020 separately and combined the phase and amplitudemaps for each telescope into a complex map. We then minimized the square absolute difference to the model prediction summed over all pixels with respect to the input coefficients. Because of the nature of the approximate coordinate scaling (step 5 of the analysis pipeline), at the edge of a map, only the smallest wavelengths contribute. We limited the radius to which the data are considered in the fit to α_{max} × λ_{low}/λ_{high}. Here, α_{max} is the size of the full map, and λ_{low} and λ_{high} are the wavelength of the lowest and highest channel, respectively. This choice ensures equal participation of all channels in the fit.
The optical layout of observations with the Calibration Unit has some important differences with the onsky situation, for which the phase maps will be applied later: the lack of a central obscuration and an enlarged outer stop r_{GCU} = 9.6 m/2 alter the shape of the pupil defined in Eq. (4). As a consequence, the Calibration Unit pupil illuminates imageplane aberrations out to a slightly larger radius. We chose to normalize the Zernike polynomials by r_{tel} = 8.0/2 m, that is, the telescope area covered by the secondary mirror, to optimize our parameterization for the onsky case. Image plane distortions, on the other hand, are normalized over the imageplane fiber width at λ_{0}, ,
Of the different types of maps we constructed, the fringe tracker provides the cleanest system and thus gives an important benchmark point for the agreement between model and data. We thus used the FT maps to determine the order n_{max} to which Zernike polynomials in the pupil and focalplane aberrations are considered. Successively increasing the fit order, we found that pupilplane aberrations with n_{max} = 6 and focalplane aberrations with n_{max} = 2 provide satisfactory model consistency while still allowing for manageable convergence times. Increasing the Zernike order in the pupil plane is especially important to reduce phase residuals at larger radii, while the central part of the maps can also be described by polynomials of lower order. Fits without focalplane aberrations reproduce the phase structure to a satisfactory degree, but show poor consistency between the phase and the amplitudedata. Finally, an additional parameter accounts for the overall amplitude scaling between measured and predicted maps, such that each fit constrains at least 34 degrees of freedom. The phase RMS achieved for the fringe tracker fits is ∼1° for all beams and data sets; extrapolation of the fit result to the full map radius yields an RMS of a few degrees.
In principle, it is possible to directly fit the SC maps by the same procedure employed for the FT. However, by further refining the analysis, we can remove additional systematic effects from the SC maps. When we created the maps, we corrected for misalignment of the image plane coordinates with the amplitude maximum (step 3 in the analysis pipeline). This shift, however, is not guaranteed to be identical on SC and FT, and as a result, there can be a small offset between the FT phase entering the differential SCFT measurement. To describe this effect, we fitted a differential map, predicted from two sets of Zernike coefficients, to the SCFT maps. The latter of these two sets of parameters was largely fixed to the previously obtained FT coefficients, and only the tiptilt terms were allowed to vary. The SC parameters, on the other hand, were all left free, such that the fit eventually determined the desired SC maps and the offset between the two channels.
From the bestfit coefficients of the differential SCFT fit, which we summarize in Appendix A, we constructed a complex SC map. Its phase is displayed in Fig. 5. As expected, the structure agrees very well with the maps obtained by direct evaluation of the Calibration Unit measurement in Fig. 4. Residuals between measured and fitted SCFT map, shown in Fig. 6, are low over the full radius considered for the fit. We obtain a bestfit RMS of 1° −2° for most beams and data sets and two slightly worse results with RMS ∼3° and ∼5°. At larger radii, the disagreement between fit and data starts to increase. This can either be caused by wavelengthdependent errors or by higher order aberrations beyond those considered for the fit. When we optimized n_{max}, we noted indeed that every increase improved the extrapolation to large separations. However, at such large offaxis distances, fiber damping becomes very significant, resulting in a poor signaltonoise ratio. We therefore consider the Zernike decomposition up to sixth order sufficient for our applications.
Fig. 5. Science channel phase maps obtained from fits to the differential SCFT maps, measured on 3 March 2020 for all four GRAVITY beams. 
Fig. 6. Phase residuals of the fit to the differential SCFT map measured on 3 March 2020 for all four GRAVITY beams. Only the data within the dashed circle are considered in the fit; the cancellation of wavelengthdependent scaling errors is not guaranteed at larger radii. 
4. Application to GRAVITY observations
Static fielddependent aberrations affect the visibility measurement whenever the size of an observed object is comparable to the fiber FOV. We applied the formalism developed in Sect. 2 alongside the characterization of aberrations from Sect. 3 to observations of two different binary systems. First we considered as a proof of concept a testcase binary observed with the Auxiliary Telescopes (ATs), where the system position in the FOV was systematically varied and thus screened over the phase and amplitudemaps. Second we applied the aberrationcorrection to GC observations with the UTs from 2017 and 2018. During these epochs, which are around the pericenter passage of S2, S2 and Sgr A* where observed simultaneously in a single fiber pointing.
The data considered in either analysis consisted of visibility amplitudes, squared visibilities, and closure phases with a relative weighting of (1:1:2). To infer the source separations, we fit a binary model based on Eq. (25), which we extended to account for the effect of finite spectral resolution and for a homogeneous background with flux ratio f^{bkg} relative to the first binary component.
Phase distortions enter this expression via the OPD correction . Furthermore, the pointsource visibility averaged over a spectral channel is
The spectral bandpass P(λ) is given by a top hat function. The source positions α_{1} and α_{2}, the flux ratios f^{bin} and f^{bkg}, and the spectral index of the central component (ν_{1}) and the background flux (ν_{bkg}) were free fit parameters, while the spectral slope of the companion was fixed to ν_{2} = 3.
Finally, , , and in Eq. (34) refer to the phase maps, amplitude maps, and the photometric lobes as they are encountered in onsky observations. These have two important differences to the Calibration Unit measurement. First, while the pupilplane representation of the aberrations is the same for both settings, the presence of a central obscuration and the smaller outer stop affects the realization of the maps in the image plane. This was conveniently captured by using the Zernike coefficients found in Sect. 3.2 to create a new set of maps with adjusted pupil configuration. Second, the maps are subject to turbulent smoothing according to Eqs. (30) and (31).
4.1. Verification for a binary test case
The testcase observations, carried out with the ATs in astrometric configuration, targeted HIP 41426, a binary with a Kband magnitude m_{K} ≃ 5.393 at RA = 8 : 26 : 57.75 h, Dec = −52 : 42 : 17.8 (Cutri et al. 2003). The system has an approximate separation of 200 mas. Its position relative to the GRAVITY fiber was kept fixed for three of the four telescopes and was varied in 24 steps between ±400 mas on AT2. At each offset, ten frames with a 6 s integration time were taken. The setup is illustrated in Fig. 7, which shows the two binary components relative to the fiber profile on all four telescopes. The shift was applied along the xaxis in the frame of the GRAVITY pupil, whose rotation with respect to the field results in a diagonal movement on the sky.
Fig. 7. Illustration of the AT binary test observations, showing the position of the two binary components (circles and diamonds, respectively) relative to the fiber profile (gray shading). Color gradients are chosen in accordance with Fig. 8. For this test, the fiber position was only varied on AT2 and was kept fixed on the other three telescopes. 
We used the Zernike coefficients obtained for the SC in Sect. 3.2 to produce phase and amplitudemaps tailored to observations with the ATs. In this case, the pupil, cf. Eq. (4), is defined by r_{tel} = 1.82 m/2 and r_{cent} = 0.14 m/2. After beam collimation, ATs and UTs illuminate the same section on the GRAVITY mirrors, such that the pupilplane phase screen can simply be scaled to the AT radius, that is, r_{tel} = 1.82 m/2 also applies in the Zernike decomposition of Eq. (5). To authenticate the effect of correct aberration modeling, we compared our results to a second nomap analysis. In this latter scenario, we set all phase maps to zero and all amplitude maps to one, , .
For fiber offsets that are too large, the signaltonoise ratio on AT2 is poor because of the high fiber damping, and we consequently discarded these data. The remaining pointings are shown in Fig. 7, and the corresponding separation, measured from a binary fit to the data according to Eq. (34), is given in Fig. 8.
Fig. 8. Binary separation inferred for a varying fiber offset on AT2 with (right panel) and without (left panel) application of the phase and amplitudemaps. Each data point shows the average over two polarization states, and the range of offsets corresponds to ±200 mas, approximately. 
The AT binary test case clearly validates our aberration corrections. Different configurations yield consistent results only if phase and amplitudemaps are considered in the analysis. Including the correct aberration model in the analysis clearly shifts the result and reduces the scatter. Even more importantly, however, the separation found in the nomap analysis systematically depends on the fiber position; it is largest for positive fiberoffsets and smallest for offsets in the negative direction. When the aberration correction is applied, this systematic is largely removed.
We considered the binary testcase observations primarily as a proof of concept and therefore did not apply a full analysis of the systematic measurement error as carried out for the GC. These uncertainties arise from the accuracy to which the phase maps can be determined and from the uncertainty of the atmospheric smoothing kernel. Furthermore, there can be minor differences in the phase and amplitudemaps between AT und UT observations, and our treatment is optimized for the UT scenario.
As the shift in its central value indicates, the binary separation is large enough that even at perfect fiber pointing, at least one source lies in a region of the FOV where aberrationinduced phase errors are significant. Accurate astrometry is thus not a question of precise fiber alignment, but is only possible with a consistent treatment of the pupilplane distortions in the analysis.
4.2. Separation between S2 and Sgr A*
After verifying our approach to correct for aberrationinduced systematic errors, we also applied it to GC observations with GRAVITY. During 2017 and 2018, that is, close to pericenter passage, S2 and Sgr A* were observed simultaneously in a single fiber pointing. In particular during 2017, when the offaxis distance of S2 was larger, the aberration correction improves the inferred binary separation. In 2019, in contrast, the Sgr A*S2 separation exceeded the single telescope beam size of about 60 mas, and GRAVITY observed both sources separately in socalled dualbeam mode. Their separation is then obtained by calibrating Sgr A* with S2 and fitting a pointsource model to its visibilities (see Gravity Collaboration 2020 for details). In this configuration, each source can be well aligned with the fiber center, such that fielddependent aberrations do not affect the measurement.
To derive the aberrationinduced shift of the S2 position, we examined a subset of the GRAVITY data used in Gravity Collaboration (2019). In particular, we applied stricter quality cuts and demanded a high signaltonoise ratio. Phase and amplitudemaps were generated from the coefficients obtained in Sect. 3.2 by accounting for the specific geometry of UT observations, that is, r_{tel} = 8.0 m/2 and r_{cent} = 0.96 m/2. The residual turbulent tiptilt is between 10 mas and 15 mas per axis (Perrin & Woillez 2019). In total, we considered four different realizations of the aberration maps that are given by the independent analysis of the two calibration measurements in 2019 and 2020 each convolved with the minimum and maximum smoothing assumption. A representative example for the phase maps applied in the GC analysis is shown in Fig. 9 in relation to the orbit of S2.
Fig. 9. Orbit of S2 relative to the phase maps as applied for the GC analysis (measurement from 3 March 2020, σ_{t} = 10 mas). Dots indicate the position of S2 on 2017.2, 2017.6, 2018.2, and 2018.7, respectively, while the cross marks Sgr A*. 
Our main result, the difference in the position of S2 with and without aberration corrections averaged per month, is shown in Fig. 10. As expected, the correction is largest early in 2017 and smallest around pericenter passage in May 2018. Furthermore, the mean corrections per epoch obtained with the four different realizations of the aberration maps are consistent over the full observational period.
Fig. 10. Difference in the position of S2 obtained from an analysis with and without the aberration corrections. Colored dots indicate the epochwise mean for different realizations of the phase and amplitudemaps, and gray dots show the results for individual observations. From these, we determine a mean position correction as a function of time with a corresponding upper and lower limit as indicated by the solid black line and the gray band. The thin dashed line finally represents the correction applied in Gravity Collaboration (2019). 
As the orbit of S2 smoothly scans over the phase and amplitudemaps (see Fig. 9), we also expect a smooth variation in the positioncorrection. The timedependence in Fig. 10 is indeed well described by a secondorder polynomial fit,
where τ = t/years − 2018.4 refers to the shifted observation date in years.
In addition to the mean correction per epoch, Fig. 10 also shows the individual filebyfile results as gray dots. These give some insight into the uncertainty of the aberration correction. When we fit the orbit of S2, any such uncertainty must to be propagated as a source of systematic error. We constructed a upper and a lower estimate of the correction, containing 68% of the files per epoch. This is shown in Fig. 10 as a gray band.
In addition to the systematic error, we also need to account for the statistical uncertainty of the position of S2. That is, as the phase and amplitudeerror changes when the position of S2 is varied within its error bars, we need to propagate this effect to the final correction. To this end, we took the position error of the original uncorrected data point, from which we drew 100 realizations, and shifted the aberration maps by this. We then derived the correction from each realization independently and used their scatter to estimate the statistical error of the correction for the position of S2. The resulting mean statistical uncertainty per epoch is small, between 10 μas and 30 μas, but we nevertheless also accounted for it in the orbit fitting.
A further check was to determine which correction created a best match of the 2017 and 2018 GRAVITY positions to the remaining S2 data. To this end, we included a scaling factor f_{corr} in the correction we applied, such that f_{corr} = 1 was our best correction and f_{corr} = 0 was no correction. This parameter was then included in the orbit fit (see Sect. 5.1). The best fit yields f_{corr} = 0.99 ± 0.06, which is identical to the correction we derived purely from calibration data. This independently confirms our concept and the resulting aberration correction: Our correction yields the most consistent S2 orbit.
The aberration correction presented here constitutes a further refinement of the analysis in Gravity Collaboration (2020). There we directly applied the measured aberration maps as shown in Fig. (4) instead of the fitted decomposition in terms of pupilplane Zernike polynomials. To account for the widening of the maps, which occurs when we project from the enlarged stop on the Calibration Unit to the telescope pupil, we applied a wider smoothing kernel of σ_{t} = (19±5) mas. The resulting best estimate for the correction is depicted in Fig. 10 as a dashed line. Both methods give consistent results, which confirms the robustness of the approach. The only sizable deviation is in 2017.2, when S2 was observed at a separation comparable to the maximum radius for which we obtained the calibration measurement (see Fig. 9). This case shows the strength of the Zernike decomposition, which allows a welldefined extrapolation.
5. Results
5.1. Determination of the S2 orbit
In the following we evaluate the effect of the aberration correction on the S2 orbit. The data used are similar to those of Gravity Collaboration (2020) and are described in detail in Appendix B. We employed the same fitting procedure as in Gravity Collaboration (2020), using a 13parameter postNewtonian orbit model. Six of these parameters describe the Kepler orbit (a, e, i, ω, Ω, and t_{peri}), and another six describe the reference frame relative to the AO spectroscopy and assumed correction for the local standard of rest (LSR) (x_{0}, y_{0}, R_{0}, ẋ_{0}, ẏ_{0}, ż_{0}). Here, R_{0} is the distance to the GC, the prime focus of this work, and M_{•} is the central mass. The bestfit parameters are given in Table 2.
Contribution to the systematic errors affecting the measurement of R_{0}, for details see Gravity Collaboration (2019).
Orbital parameters of S2 with their statistical uncertainties.
To determine the systematic uncertainty, we followed the approach in Gravity Collaboration (2019) of varying our assumptions and tracing the associated changes in R_{0}. Compared to our earlier work, we also included the uncertainty due to the aberration correction, as given by the gray band in Fig. 10. The individual contributions are given in Table 1. The aberration correction is the dominant contributor to the systematic error. The total systematic uncertainty is 33 pc when the contributions are added quadratically.
Our best estimate of the GC distance thus is
5.2. Comparison to previous results
Our previous determinations of the GC distance in Gravity Collaboration (2018, 2019), and Gravity Collaboration (2020) were biased by the fielddependent aberrations. Taking them into account brings all our measurements into agreement, as shown in Fig. 11 and Table 3. In contrast to Gravity Collaboration (2020), we here also applied a correction for the 2018 data, when S2 and Sgr A* were close to each other and close to the field center. We note that, the small aberration corrections lead to a small upward correction of R_{0} of around 30 pc. However, this shift is comparable to the systematic error. Further, we find that the orbit is particularly sensitive to the pericenter data. This causes a strong decrease in the statistical uncertainty with time, while the systematic uncertainty even increases slightly because varying the assumptions then leads to stronger variations in the fit result.
Fig. 11. Measurements of the GC distance over time with a focus on studies of the S2 orbit. Blue points show results obtained with the SINFONI, NACO, and GRAVITY data with (dark blue) and without (light blue) application of the aberration corrections. Gray R_{0} determinations are based on data from the Keck observatory. For comparison, we show in black results that are based on the statistical parallax of the nuclear star cluster (Chatzopoulos et al. 2015) and from modeling the Milky Way dynamics based on observations of molecular masers (Reid et al. 2019). BlandHawthorn & Gerhard (2016) finally reported the GC distance based on a combination of various methods. 
Published values of R_{0} (bold) and the corresponding values when the aberrations are taken into account (right column).
5.3. Comparison with further S2based results
We estimate that the accuracy of our VLTbased result is at the 40 pc level. However, it deviates significantly from the Keckbased value reported in Do et al. (2019), with the difference being at the 300 pc level. Because both works use the orbit of S2 around Sgr A* to determine R_{0}, it is important to investigate the cause of the discrepancy, and we address this in Appendix C. Overall, we conclude that the combination of a difference in the radial velocity data and a modest offset of the Keck coordinate system in the declination direction might explain the discrepancy. Each effect contributes roughly 50%.
About 20% of the radial velocity difference can be attributed to the Doppler formula in StarKit that was implicitly used by Do et al. (2019). The remaining 80% are unexplained and might lie either in the Keck or the VLT data.
The origin of the coordinate system offset is unclear as well. Trying to explain the offset with a shift in the VLT coordinate systems is much harder than imposing a shift in the Keck system because the precision of the GRAVITY data is so high.
6. Conclusions
GRAVITY delivers highresolution astrometry, which in combination with spectroscopic data, allows a very precise determination of the distance to the Galactic Center. The values inferred from different epochs (Gravity Collaboration 2018, 2019, 2020) show a small discrepancy at the 1% level, which nevertheless is significant because of the high precision of the measurement.
We were able to relate this shift to optical aberrations introduced in the instrument, which lead to a fielddependent distortion of the visibility phase. Their effect is stronger the farther offaxis an object lies within the FOV. In particular, GC observations close to the S2 pericenter passage are affected, when S2 and Sgr A* are detected simultaneously in a single fiber pointing but at a separation comparable to the FOV. In earlier and later epochs, in contrast, we employed the socalled dualbeam method and targeted each source individually. In this case, as for most other GRAVITY science observable, each source can be well centered and aberration corrections become irrelevant. The dualbeam observation mode was also assumed to derive the astrometric error budget in Lacour et al. (2014), who did not include the effect of phase maps for this precise reason.
The full analytical description we developed here allowed us to propagate the effect of optical aberrations at fiber injection to the measured visibilities. Fitting this model to dedicated calibration measurements confirmed its validity and enabled us to account for the effect in the data analysis. We further verified the approach with dedicated testcase observations.
The formalism we developed is applicable beyond GRAVITY to any optical/nearIR interferometer in which aberrations are introduced in the pupil or the focal plane. In several cases in the literature, more than one object lay in the interferometer FOV, for example, some Keck (Colavita et al. 2013), CHARA (Ten Brummelaar 2005), or NPOI (Armstrong et al. 1998) results on binary stars. How severely aberrations affect an observation depends not only on their strength for a particular instrument, however, but also on the offaxis distance considered and on the statistical noise in the measurement. In the example of GRAVITY on the UTs, the mean phase error introduced at 20 mas separation is 4−5 degrees per telescope and increases to 14−20 degrees at 50 mas. While a binary test case as presented in Sect. 4.1 can serve as a general strategy to diagnose whether aberrationinduced systematics are a problem, dedicated calibration measurement are required for their correction in the analysis for each individual instrument.
With the results from the GRAVITY Calibration Unit measurements and our refined analysis scheme, we were able to further improve the separation between S2 and Sgr A* in 2017 and 2018, introducing shifts up to 0.5 mas caused by the phase aberrations. In Fig. 12 we show a detailed view of the S2 orbit in 2017, in which we also included two dualbeam measurements that do not suffer from phase aberrations. The improved data agrees very well with these positions.
Fig. 12. Detailed view of the S2 orbit in 2017. Dualbeam points do not suffer from aberrationrelated systematic errors and agree very well with our corrected data points. 
Of all orbital parameters, the distance to the GC R_{0} is most strongly affected by the change in the S2 position. This can be easily understood when R_{0} is viewed as the scaling factor between angular and proper velocity. The fielddependent phase errors discussed in this work then fully explain the shift between earlier R_{0} measurements with GRAVITY data. Applying the analysis scheme developed here lifts any such discrepancies (see Sect. 5.2). In particular, Fig. 11 demonstrates that later corrected data sets of earlier publications give fully consistent results whose accuracy increases with time.
References
 Armstrong, J. T., Mozurkewich, D., Rickard, L. J., et al. 1998, ApJ, 496, 550 [Google Scholar]
 Bhatnagar, S., Cornwell, T. J., Golap, K., & Uson, J. M. 2008, A&A, 487, 419 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 BlandHawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Blind, N., Eisenhauer, F., Haug, M., et al. 2014, Proc. SPIE, 9146, 91461U [Google Scholar]
 Boehle, A., Ghez, A. M., Schödel, R., et al. 2016, ApJ, 830, 17 [Google Scholar]
 Chatzopoulos, S., Fritz, T. K., Gerhard, O., et al. 2015, MNRAS, 447, 948 [Google Scholar]
 Chu, D. S., Do, T., Hees, A., et al. 2018, ApJ, 854, 12 [Google Scholar]
 Colavita, M. M., Wizinowich, P. L., Akeson, R. L., et al. 2013, PASP, 125, 1226 [Google Scholar]
 Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, VizieR Online Data Catalog: II/246 [Google Scholar]
 Do, T., Hees, A., Ghez, A., et al. 2019, Science, 365, 664 [NASA ADS] [CrossRef] [Google Scholar]
 Eisenhauer, F., Genzel, R., Alexander, T., et al. 2005, ApJ, 628, 246 [Google Scholar]
 Genzel, R., Eisenhauer, F., & Gillessen, S. 2010, Rev. Mod. Phys., 82, 3121 [Google Scholar]
 Ghez, A. M., Duchêne, G., Matthews, K., et al. 2003, ApJ, 586, L127 [Google Scholar]
 Ghez, A., Salim, S., Weinberg, N. N., et al. 2008, ApJ, 689, 1044 [Google Scholar]
 Gillessen, S., Eisenhauer, F., Fritz, T. K., et al. 2009a, ApJ, 707, L114 [Google Scholar]
 Gillessen, S., Eisenhauer, F., Trippe, S., et al. 2009b, ApJ, 692, 1075 [Google Scholar]
 Gillessen, S., Plewa, P. M., Eisenhauer, F., et al. 2017, ApJ, 837, 30 [Google Scholar]
 Gravity Collaboration (Abuter, R., et al.) 2017, A&A, 602, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gravity Collaboration (Abuter, R., et al.) 2018, A&A, 615, L15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gravity Collaboration (Abuter, R., et al.) 2019, A&A, 625, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gravity Collaboration (Abuter, R., et al.) 2020, A&A, 636, L5 [CrossRef] [EDP Sciences] [Google Scholar]
 Habibi, M., Gillessen, S., Martins, F., et al. 2017, ApJ, 847, 120 [Google Scholar]
 Lacour, S., Eisenhauer, F., Gillessen, S., et al. 2014, A&A, 567, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lindegren, L., & Dravins, D. 2003, A&A, 401, 1185 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Martins, F., Gillessen, S., Eisenhauer, F., et al. 2008, ApJ, 672, L119 [NASA ADS] [CrossRef] [Google Scholar]
 Neumann, E. G. 1988, Singlemode Fibers [Google Scholar]
 Perrin, G., & Woillez, J. 2019, A&A, 625, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pfuhl, O., Haug, M., Eisenhauer, F., et al. 2014, Proc. SPIE, 9146, 914623 [Google Scholar]
 Plewa, P. M., & Sari, R. 2018, MNRAS, 476, 4372 [Google Scholar]
 Plewa, P. M., Gillessen, S., Eisenhauer, F., et al. 2015, MNRAS, 453, 3234 [Google Scholar]
 Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2019, ApJ, 885, 131 [Google Scholar]
 Roddier, F. 1981, Prog. Opt., 19, 281 [Google Scholar]
 Smirnov, O. M. 2011, A&A, 527, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Smirnov, O. M., & Tasse, C. 2015, MNRAS, 449, 2668 [NASA ADS] [CrossRef] [Google Scholar]
 Tasse, C., Hugo, B., Mirmont, M., et al. 2018, A&A, 611, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ten Brummelaar, T. A., et al. 2005, ApJ, 628, 453 [Google Scholar]
 Wallner, O., Leeb, W. R., & Winzer, P. J. 2002, J. Opt. Soc. Am. A, 19, 2445 [Google Scholar]
Appendix A: List of Zernike coefficients
The Zernike coefficients obtained by fitting the Calibration Unit measurements from late 2019 and early 2020 are summarized in Tables A.1 and A.2, respectively. We provide the science channel results for all for GRAVITY beams (GV1 to GV4) in units of μm according to the definitions in Eqs. (5) and (33), where labels pupilplane aberrations and those in the focal plane.
Zernike coefficients for sciencechannel aberrations fit to the calibration measurement on 3 November 2019.
Zernike coefficients for sciencechannel aberrations fit to the calibration measurement on 3 March 2020.
Appendix B: Data
We used the data set presented in Gravity Collaboration (2020) with the following changes. Each singlebeam astrometric position was corrected according to Eq. (37), and we added the statistical error of this correction in quadrature, which increased the individual uncertainties by around 15 μas. Further, we corrected the radial velocity of the epoch 2018.1277, which was 13 km s^{−1} too high in the previous data set. Finally, we were able to add one interferometric position measurement of S2 from early March 2020. Like in 2019, the separation between S2 and Sgr A* exceeds the fiber FOV, and hence a dualbeam measurement needed to be employed.
Because the observability of the GC is limited in early March and we expected observations in the following months, we did not attempt to observe Sgr A* in March 2020, but only pointed to S2 and to our usual calibrator star R2, with the aim of testing the stability of the GRAVITY astrometry. Pointings to Sgr A* were planned for later in the year. They had to be canceled due to the pandemicrelated closure of the VLT(I) from midMarch on. To still determine the S2 – Sgr A* separation vector from this observation, we needed to proceed in two steps and first measured the S2 – R2 distance, then we referenced R2 to Sgr A*.
The distance between S2 and R2 was measured with the dualbeam method (Sect. 4.2), where we calibrated the S2 files with R2. In addition to the 2020 measurement, this separation is also available for 56 epochs in 2017, 2018, and 2019. It can be measured very precisely because the two stars are bright. Because the S2 – Sgr A* vectors were previously determined in Gravity Collaboration (2020), we can also refer R2 to Sgr A* in these earlier epochs. We then fit a simple quadratic function for the time evolution of the R2 coordinates relative to Sgr A* and extrapolated it to March 2020. Given the large number of data and the small time range to extrapolate for, the additional uncertainty that this introduces is well below the 100 μas level.
We derived the S2 position in 2020 from the four scientifically usable exposures as their mean. We assigned an error of 150 μas to each coordinate for this data point, which reflects both the smaller number of files compared to what we typically had available in 2019 and the additional uncertainty due to the additional step of referencing via R2. The new data point falls well onto the expected orbit, but its error bar is too large to have a significant effect on the fitted parameters.
Our data set consists of 128 AObased astrometric points, 58 GRAVITYbased astrometric points, and 97 radial velocities, of which the first three before 2003 are from Do et al. (2019).
Appendix C: Analysis of the difference between R_{0} determinations from Keck and VLT data sets
We expect that our determination of R_{0} is accurate to the 40 pc level, but we note that the value published in Do et al. (2019) is discrepant at the 300 pc level. Both teams used the orbit of the star S2 around Sgr A* to determine R_{0}, and it is therefore natural to ask what caused the differences.
C.1. Data
In addition to our (VLT) data set (Appendix B), we used the Keck data set published in Do et al. (2019). We applied the NIRC2 radial velocity offset of +80 km s^{−1} as determined in Do et al. (2019) to the NIRC2 data, that is, we added 80 km s^{−1} to these radial velocities. Unlike Do et al. (2019), we then did not fit for this offset. Furthermore, we dropped the last astrometric data point (epoch 2018.67148268), as suggested by the authors in a private communication. The data set consisted of 45 astrometric points and 116 radial velocities, of which 41 are actually from the VLT data set between 2003 and 2016. The published table also includes one radial velocity from the epoch 2019.3567, which may not have been part of the data set used in Do et al. (2019).
C.2. Difference in R_{0}
We fit the orbit with a simple 13parameter model: The six orbital elements of the star (corresponding to the initial conditions of the star in phase space), six parameters for the position and velocity of the MBH, and the mass of the MBH. The fits were made using the relativistic corrections as in Gravity Collaboration (2020), that is, we fixed f_{RS} = f_{SP} = 1. For this nonKeplerian motion, the meaning of the orbital elements is that they osculate at a reference epoch, for which we chose T = 2010.35, which is close to the apocenter passage time of S2.
To fit the VLT data set, we used the same approach as in Gravity Collaboration (2020): for the GRAVITY data, we assumed that the astrometry directly refers the S2 positions to the mass center and we directly measure the separation vector between the two objects interferometrically. For the NACO (AOimaging based) data, we allowed a coordinate system offset, on which we set priors following Plewa et al. (2015), and we included the NACO flare positions as an additional constraint to locate the mass. This fit yields
where a is the semimajor axis, i is the inclination, and Ω is the position angle of the ascending node of the S2 orbit, and the errors are the statistical fit uncertainties. The VLT astrometry is dominated by the GRAVITY points, as illustrated by dropping all AO data points, which results in R_{0} = 8276 ± 10 pc.
Fitting the Keck data set with the same 13parameter model as used for Eq. (C.1) yields
This is not the same number as in Do et al. (2019), who reported R_{0} = 7959 ± 59 pc. The small (and statistically insignificant) difference is most likely due to the noise model that Do et al. (2019) included in their analysis, which we do not have readily available. Applying the noise model at hand (Plewa & Sari 2018; Gravity Collaboration 2019) yields R_{0} = 7965 ± 56 pc. This shows that the value reported by Do et al. (2019) lies between the two numbers we obtain by refitting their data. In the following, we use for simplicity and for equal treatment of the data the value and approach as in Eq. (C.2). We therefore have a difference of ΔR_{0} = 340 ± 45 pc.
C.3. Comparing, combining, and adjusting the astrometry
Gillessen et al. (2009a) have noted that a simple attempt of comparing the astrometric data sets by plotting them on top of each other fails. We need to allow for an offset and a drift between the two coordinate systems (i.e., four parameters Δx, Δy, Δv_{x}, and Δv_{y}). This yields a 17parameter fit. Comparing the bestfitting parameters in Eqs. (C.1) and (C.2) shows that they differ significantly in Ω. This parameter is fully degenerate with the angular orientation (called β here) of the coordinate system. The difference in Ω thus suggests that the two astrometric data sets are rotated with respect to each other.
Therefore we extend the combination scheme by an additional fifth parameter, Δβ, resulting in an 18parameter fit. With this, we fit the two data sets simultaneously, omitting the 41 VLT radial velocities from the Keck data set while also dropping the three Keck velocities in the VLT data set. This fit matches the two coordinate systems ideally onto each other and results in
The value of Δβ matches the difference ΔΩ. We conclude that the Keck and VLT data are indeed rotated with respect to each other. The other parameters are very similar to Eq. (C.1), which is due to the considerably smaller astrometric uncertainties of the GRAVITY data compared to the adaptive optics data.
With the bestfit coordinate system difference in hand, we can transform the Keck astrometric data into the VLT coordinate system and vice versa. We chose to do the former because the VLT data set is more directly calibrated by the interferometric data. After applying the coordinate system difference to the Keck data, we fit them again with a 13parameter model. This yielded the same bestfit parameters as in Eq. (C.2) (with the exception of Ω). This shows that transforming the astrometry does not change the more fundamental differences between the two orbits, while a direct comparison is now feasible. The value of Ω can be omitted in the following.
C.4. Discrepancy in the radial velocity data
Chu et al. (2018) have investigated the consistency of the radial velocity data between the Keck and VLT data sets for 2000 to 2016 and concluded that the data agree with each other. We repeated the exercise, now also extending into the time of the pericenter passage in 2018 (Fig. C.1). To our surprise, the radial velocities differ systematically from ∼2011 on, and the difference increases as the radial velocity continues to increase. The difference reaches ∼50 km s^{−1} in 2018, just before the star passed pericenter. Moreover, there is one obvious outlier in the Keck data, the earliest 2018 point. We verified that dropping this measurement does not change the Keck fit result in any significant way.
Fig. C.1. Comparison of the radial velocity data sets. Blue points are data from the VLT data set, and red points are data from the Keck data set. Top left: radial velocity as a function of time for the VLT fit (Eq. (C.1)). Top right: yearly averages of the residua of the two data sets to the fit from Eq. (C.1). By construction, the VLT data thus scatter around 0. The Keck data deviate systematically from 2011 on, and the discrepancy increases in later years. Bottom left: same as the left panel, but zooming in to the period 2015–2020, and showing all individual data points. The bestfit Keck orbit corresponding to Eq. (C.2) is the red line. The difference is apparently largest when the radial velocity is highest (in 2018 at pericenter passage). Bottom right: both data sets show a clear peak in radial velocity in 2018 compared with the Keplerian part of the VLT fit (Eq. (C.1)), i.e., both data sets clearly detect the redshift term. 
The obvious question is what the radial velocity effect is on R_{0}. To answer this, we swapped the radial velocities between the two data sets. The VLT set together with the Keck astrometry yields
Conversely, the Keck radial velocities together with the VLT astrometry yield R_{0} = 8214 ± 14 pc. Because the Keck radial velocity set contains 35% of the VLT radial velocities, the fit in Eq. (C.4) is the cleaner test. We thus explain roughly half of the difference in R_{0} with the radial velocity data, that is, 159 pc.
So far, we could only find an explanation for ∼20% of the radial velocity difference: We applied the fitting based on the stellar atmosphere model with the StarKit package used in Do et al. (2019) also to the VLT spectroscopy. We found a significant difference for high radial velocities, which we were able to trace back to the Doppler formula used by the StarKit package. While both Do et al. (2019) and Gravity Collaboration (2020) stated that the spectroscopic observable is v_{r} = z c, that is, the redshift of a given spectrum, the StarKit package applies a Doppler formula that includes the longitudinal relativistic correction: . In this form, the Doppler formula ignores the (significant) tangential motion v_{t} of S2. In order to apply a relativistic correction, the full Doppler formula (Lindegren & Dravins 2003) is required. For this correction, however, the spectroscopic information is not sufficient. In general, a spectrum cannot be Dopplercorrected in a relativistic way without knowledge of the other motion component. Furthermore, even if the full correction were applied, it would be impossible to fit for the relativistic redshift later.
The difference between the two formulae is small at velocities far lower than the speed of light, but it becomes large close to pericenter, when S2 reaches a velocity of nearly 8000 km s^{−1}. It amounts to ∼25 km s^{−1} at most, however, and thus is smaller than the observed difference in Fig. C.1. This difference is also visible in Fig. 1 of Do et al. (2019): The plotted model spectra are slightly more redshifted than the underlying data suggest. Changing the Keck radial velocities accordingly yields a fit with R_{0} = 7972 ± 44 pc, which accounts for 37 pc of the 159 pc.
Further checks did not yield any clues for the remaining significant difference in the radial velocities. In this context, we verified whether the time stamps were assigned consistently between the two data sets, and did not find a difference. We also note that in the bottom right panel in Fig. C.1 the redshift peak around pericenter is clearly visible for both data sets.
C.5. Discrepancy in the astrometry
Comparison of the fits in Eqs. (C.1) and (C.2) shows that they differ not only in R_{0}, but also in the size of the semimajor axis a. We find Δa/a = 1.28 ± 0.22%. This does not hold for the semiminor axis, however, Δb/b is consistent with 0. Interestingly, the projected ellipses as given by the astrometric data in the plane of sky agree in both semimajor and semiminor axes to within 0.17%. Hence, the inclinations i need to differ, which Eqs. (C.1) and (C.2) confirm. We find in accordance with the above 1 − sin(i_{VLT})/sin(i_{Keck}) ∼ 1.3%.
The inclination of the ellipse determines where the projected center of mass is located. The orientation of the S2 orbit and the disagreement in a but not in b indicate an offset of the center of mass in the declination direction. We show that introducing an offset to either y or v_{y} (the mass position and velocity in declination) can indeed explain the remaining discrepancy. Starting from the fit of the transformed Keck data set, we fixed v_{y} to its bestfit value of −0.15 mas yr^{−1}. All other parameters were again left free for a subsequent fit. Additionally using the VLTI velocities in this fit instead of the Keck velocities yields
This fit therefore yields the same value for R_{0} from the Keck astrometry as the VLT fit. Moreover, the semimajor axis a and inclination i have moved to the VLT values when v_{y} is forced to have an offset. Because the mass position is parameterized with a time origin at T = 2000.0, the bestfit y also changes, from −0.972 mas to 1.234 mas. The systematic uncertainty on y and v_{y} estimated by Do et al. (2019) is 1.16 mas and 0.066 mas yr^{−1}, respectively. This means that the difference that needs to be enforced is within ∼2σ of the systematic uncertainty, and the residuals in Fig. C.2 (left) appear to be acceptable. Essentially the same can be achieved by forcing an offset to y and leaving v_{y} free instead.
Fig. C.2. Comparison of the astrometric residual after forcing an offset in declination such that the fit to the Keck data set matches that to the VLT (left) and such that the fit to VLT data set matches that to the Keck data (right). Lighter blue corresponds to AO data from the VLT data set, and darker blue to the GRAVITY data. 
This argument might be reversed and a similar offset be applied to the VLT data in order to lower the VLTbased value of R_{0}. In a first attempt, we applied the same offset to the VLT AO data. However, even an offset 10 x larger (i.e., 1.2 mas yr^{−1}) changed R_{0} by only ∼30 pc. This is expected because the VLT astrometry is completely dominated by the GRAVITY data. We therefore instead tried to vary v_{y} and y for the GRAVITY data, giving up the assumption that the GRAVITY source itself is the mass center. We also exchanged the VLT radial velocities for the Keck velocities. We find that we need to change v_{y} by −1.4 mas yr^{−1} in order to obtain a similar distance to the Keck value:
The fit achieves the lower R_{0} by tilting the orbit similar to the fit from Eq. (C.2). The enforced change of v_{y} is unrealistically large (12× larger than what was needed for the Keck data). The GRAVITY data also show very strong and systematic residuals of up to 0.5 mas (right panel in Fig. C.2), and the reduced χ^{2} of the fit increased from 1.50 to 2.63.
All Tables
Contribution to the systematic errors affecting the measurement of R_{0}, for details see Gravity Collaboration (2019).
Published values of R_{0} (bold) and the corresponding values when the aberrations are taken into account (right column).
Zernike coefficients for sciencechannel aberrations fit to the calibration measurement on 3 November 2019.
Zernike coefficients for sciencechannel aberrations fit to the calibration measurement on 3 March 2020.
All Figures
Fig. 1. Schematic depiction of the pupil and focal plane aberrations that enter the overlap integral. Both effects in combination are required to describe the aberration patterns observed in calibration measurements. The lowest order aberrations in the pupil function are shown explicitly, which are (A) lateral fiber misplacement, (B) fiber tilt, and (C) defocus. Their effect is further explained in the text. 

In the text 
Fig. 2. Example phase screens (top) and amplitude maps (bottom) in the image plane induced by loworder Zernike aberrations in the pupil plane at a wavelength of λ_{0} = 2.2 μm. From left to right: the considered aberrations are a perfect Airy pattern, a vertical tilt of 0.4 μm RMS, a vertical astigmatism of 0.2 μm RMS, a vertical coma of 0.2 μm RMS, and a combination of astigmatism, coma, and trefoil (with RMS 0.2 μm, 0.2 μm, and 0.1 μm, respectively). The rightmost panel also considers an additional fiber tilt with 0.2 μm RMS. 

In the text 
Fig. 3. Examples of the scanning pattern applied in the Calibration Unit measurements. SC aberration maps were obtained with a slow modulation frequency (left). For the corresponding FT measurement, a faster scanning was used, and the right panel only shows a single iteration of in and outspiral. 

In the text 
Fig. 4. Science channel phase maps reconstructed by the procedure of Sect. 3.1 from the Calibration Unit measurement on 3 March 2020 for all four GRAVITY beams. 

In the text 
Fig. 5. Science channel phase maps obtained from fits to the differential SCFT maps, measured on 3 March 2020 for all four GRAVITY beams. 

In the text 
Fig. 6. Phase residuals of the fit to the differential SCFT map measured on 3 March 2020 for all four GRAVITY beams. Only the data within the dashed circle are considered in the fit; the cancellation of wavelengthdependent scaling errors is not guaranteed at larger radii. 

In the text 
Fig. 7. Illustration of the AT binary test observations, showing the position of the two binary components (circles and diamonds, respectively) relative to the fiber profile (gray shading). Color gradients are chosen in accordance with Fig. 8. For this test, the fiber position was only varied on AT2 and was kept fixed on the other three telescopes. 

In the text 
Fig. 8. Binary separation inferred for a varying fiber offset on AT2 with (right panel) and without (left panel) application of the phase and amplitudemaps. Each data point shows the average over two polarization states, and the range of offsets corresponds to ±200 mas, approximately. 

In the text 
Fig. 9. Orbit of S2 relative to the phase maps as applied for the GC analysis (measurement from 3 March 2020, σ_{t} = 10 mas). Dots indicate the position of S2 on 2017.2, 2017.6, 2018.2, and 2018.7, respectively, while the cross marks Sgr A*. 

In the text 
Fig. 10. Difference in the position of S2 obtained from an analysis with and without the aberration corrections. Colored dots indicate the epochwise mean for different realizations of the phase and amplitudemaps, and gray dots show the results for individual observations. From these, we determine a mean position correction as a function of time with a corresponding upper and lower limit as indicated by the solid black line and the gray band. The thin dashed line finally represents the correction applied in Gravity Collaboration (2019). 

In the text 
Fig. 11. Measurements of the GC distance over time with a focus on studies of the S2 orbit. Blue points show results obtained with the SINFONI, NACO, and GRAVITY data with (dark blue) and without (light blue) application of the aberration corrections. Gray R_{0} determinations are based on data from the Keck observatory. For comparison, we show in black results that are based on the statistical parallax of the nuclear star cluster (Chatzopoulos et al. 2015) and from modeling the Milky Way dynamics based on observations of molecular masers (Reid et al. 2019). BlandHawthorn & Gerhard (2016) finally reported the GC distance based on a combination of various methods. 

In the text 
Fig. 12. Detailed view of the S2 orbit in 2017. Dualbeam points do not suffer from aberrationrelated systematic errors and agree very well with our corrected data points. 

In the text 
Fig. C.1. Comparison of the radial velocity data sets. Blue points are data from the VLT data set, and red points are data from the Keck data set. Top left: radial velocity as a function of time for the VLT fit (Eq. (C.1)). Top right: yearly averages of the residua of the two data sets to the fit from Eq. (C.1). By construction, the VLT data thus scatter around 0. The Keck data deviate systematically from 2011 on, and the discrepancy increases in later years. Bottom left: same as the left panel, but zooming in to the period 2015–2020, and showing all individual data points. The bestfit Keck orbit corresponding to Eq. (C.2) is the red line. The difference is apparently largest when the radial velocity is highest (in 2018 at pericenter passage). Bottom right: both data sets show a clear peak in radial velocity in 2018 compared with the Keplerian part of the VLT fit (Eq. (C.1)), i.e., both data sets clearly detect the redshift term. 

In the text 
Fig. C.2. Comparison of the astrometric residual after forcing an offset in declination such that the fit to the Keck data set matches that to the VLT (left) and such that the fit to VLT data set matches that to the Keck data (right). Lighter blue corresponds to AO data from the VLT data set, and darker blue to the GRAVITY data. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.