© G. Heißel et al. 2022

1. Introduction

The motions of the S-stars, which make up the nuclear cluster at the centre of the Milky Way, have been monitored closely now for almost three decades with respect to astrometry and radial velocity (Eckart & Genzel 1996; Ghez et al. 1998, 2003, 2008; Schödel et al. 2002; Gillessen et al. 2009, 2017). Their trajectories revealed the existence of a dark compact ∼4 × 10⁶ M_⊙ body at the cluster centre, in coincidence with the location of the radio source Sgr A^* (Reid et al. 2009; Plewa et al. 2015; Reid & Brunthaler 2004, 2020). Observations are currently in best agreement with the notion that the compact body is a massive black hole (MBH) (Genzel et al. 2010, Sect. IV). The combination of the high-eccentricity and low-pericentre distance of the orbit of star S2, as well as its sufficiently bright magnitude of ∼14 in the K band, predestined this object as the prime probe for the task of constraining the MBH parameters and observing relativistic effects (Grould et al. 2017; Waisberg et al. 2018). Flares of radiation near the innermost stable circular orbit of the MBH are of similar importance in this respect (GRAVITY Collaboration 2018a). With regard to S2, one milestone has been the detection of the gravitational redshift together with the special relativistic transverse Doppler shift during the time of pericentre passage in 2018 (GRAVITY Collaboration 2018b; Saida et al. 2019; Do et al. 2019). More recently, the Schwarzschild precession of the stellar orbit could be measured at 12 arcmin (arcmin) per orbital period by GRAVITY Collaboration (2020). Schwarzschild precession thereby refers to the relativistic pericentre advance in the orbital plane due to the mass of a non-spinning black hole (Merritt 2013; Poisson & Will 2014).

In general, an in-plane orbital precession can also be caused by other influences. In particular, if the star is moving through a (dark) extended continuous mass distribution (in the following, simply dark mass or extended mass), for instance, one that consists of faint stars and stellar remnants or dark matter, then the acceleration imposed by this extended mass causes a pericentre shift even from a Newtonian perspective (Jiang & Lin 1985; Rubilar & Eckart 2001; Merritt 2013, Sect. 4.4). This phenomenon is commonly referred to as mass precession, and it has been demonstrated that it can be of the same order as or may even dominate the Schwarzschild precession for physically reasonable density distributions (Rubilar & Eckart 2001; Merritt 2013, p. 208). For the case of the galactic centre of the Milky Way, the uncertainties of current observational data allow one to give upper bounds on a possibly present extended mass within the apocentre of S2. For instance, a 0.3 arcsec (as) Plummer profile mass is limited to below 0.1% of the mass of the MBH, which translates into ∼4000 M_⊙ (GRAVITY Collaboration 2020, Sect. 4)¹. For completeness we note that the related search for compact (as opposed to extended) masses near the MBH is also of great interest, in particular regarding potential gravitational wave sources from extreme mass ratio inspirals for Laser Interferometer Space Antenna (LISA) (Amaro-Seoane et al. 2017; Gourgoulhon et al. 2019).

There are several motivations to continue to study, constrain, and if possible, detect an extended mass within the orbit of S2. Firstly, from this we can directly infer details about the nature of the matter content in the immediate vicinity of the MBH. Secondly, more data and a refined analysis may help to place constraints on or rule out central object models based on continuous dark matter distributions that have been proposed in addition, or as alternatives, to the MBH paradigm (Lacroix 2018; GRAVITY Collaboration 2019; Becerra-Vergara et al. 2020, 2021a,b; Nampalliwar et al. 2021; Argüelles et al. 2022). Finally, beyond the interest in its own right, it is also important to study the mass precession as perturbation for the measurement of other effects. On secular timescales, this concerns in particular other precession effects such as the Schwarzschild and Lense-Thirring precessions. The latter is caused by the frame-dragging due to the spin of the black hole (Wex & Kopeikin 1999; Merritt 2013; Zhang et al. 2015; Yu et al. 2016; Grould et al. 2017; Waisberg et al. 2018; Qi et al. 2021). In contrast to the Schwarzschild precession, the spin generally precesses the orbit within its orbital plane as well as out of its orbital plane. However, so does the mass precession when the distribution is not spherically symmetric (Merritt 2013). Because of this interference, a separation of these effects solely based on the measurement of the overall accumulated secular precession of a single stellar orbit may be difficult².

We restrict our considerations to spherically symmetric extended masses that cause a pure in-plane precession, just like the Schwarzschild precession. We study how both effects can be clearly separated from each other based on their distinct astrometric signatures within a single orbit, even though their secular precessions directly interfere with each other. To this end, we first study the functional form of the argument of pericentre and of other osculating orbital elements versus true anomaly over one orbit. We underline the differences in these functional forms and identify the orbital sections in which either effect is mostly active or inactive respectively. We then extend this investigation to the domain of the observables: astrometry and radial velocity. We also support and complement our theoretical investigations by an extensive mock-data analysis. Although we do consider radial velocity, our key insights regard astrometry. A complementary approach with a focus on the effects on gravitational redshift, radial velocity, and the timing of the pericentre passage is taken in the recent work of Takamori et al. (2020). The important topic of the detectability and separability of precession effects on timescales of a single orbit has also previously been investigated in earlier work through different lenses and by different techniques (e.g. Angélil & Saha 2014; Parsa et al. 2017).

We start in Sect. 2 with a description of our MBH + star + extended mass model. In Sect. 3 we apply this model in order to discuss the different effects of the relativistic correction and of the extended mass onto the stellar orbit. We also give prospects on the ability to improve on current upper bounds or to detect a dark mass if present. After this, we proceed with a mock-data analysis to support these claims in Sect. 4. Finally, we conclude with a discussion and prospects in Sect. 5.

2. Model for an MBH + star + extended mass

Below we introduce our theoretical model as a perturbed Kepler problem (Sect. 2.1) and give some background about the osculating equations as a formalism tailored to the quasi-Keplerian nature of the problem (Sect. 2.2). In this framework, we then give the explicit perturbative accelerations for the modelled effects (Sects. 2.3 and 2.4). Some further general remarks about the model can be found in Appendix B.

2.1. Perturbed Kepler problem

With its high mass, the MBH dominates the gravitational field in the S-star cluster. Consequently, the motion of each star in its vicinity is approximately governed by the Newtonian two-body problem MBH + star, such that the orbits are quasi-Keplerian (Merritt 2013, Sect. 4.2; Poisson & Will 2014, Sect. 3.3). The deviation from Keplerian motion is thereby due to relativistic corrections of the two-body problem on the one hand, and due to the perturbative gravitational forces of other bodies on the other hand. The latter may include other compact objects, but also continuous extended mass distributions, which is what we focus on here.

For the analysis to follow, we work with a perturbed Kepler orbit model for a star of mass m_s in the gravitational field created by an MBH of mass M_• together with a spherically symmetric continuous extended mass distribution centred at the location of the MBH. We assume that m_s is much lower than M_• and the extended mass. Hence we treat the star as a test particle such that formally, m_s = 0. This also ensures consistency with our assumption that the extended mass stays centred at the MBH. Placing the origin of the coordinate system at the MBH (Fig. 1) and denoting the stellar position by r, its equations of motion then take the form

(1)

Fig. 1. Definition of the Euler angles (Ω, ι, ω) giving the orbital orientation and the true anomaly f that gives the current position within the orbit. The static fundamental frame (X, Y, Z) (blue) is adapted to the observables, and the co-rotated Gaussian frame (n, λ, ez) (orange) is adapted to the symmetry of the (perturbed) two-body problem. See Fig. C.1 for an alternative convention.

with a two-component perturbative acceleration a_p = a_1PN + a_XM. The first term on the right-hand side of Eq. (1) denotes the dominant Newtonian two-body acceleration, which is responsible for the main Keplerian component of the stellar motion, and we denote r = |r|,n = r/r. G is the gravitational constant. The two accelerations a_1PN, a_XM are responsible for perturbations from Keplerian motion. a_1PN thereby denotes the relativistic correction to the Newtonian two-body equations of motion to first post-Newtonian (1PN) order (Merritt 2013; Poisson & Will 2014). a_XM, on the other hand, denotes the perturbative acceleration due to the presence of an extended mass distribution. We return to the explicit form of these terms in Sects. 2.3 and 2.4. Before this, we switch to a formalism that is tailored to the quasi-Keplerian nature of the problem.

2.2. Osculating equations

While the equations of motion in the form of Newton’s second law (Eq. (1)) are well suited in order to demonstrate the underlying physics of the problem, other formalisms are better tailored to its quasi-Keplerian nature. One such is the post-Keplerian formalism of Damour & Deruelle (1985), which is well established in pulsar-timing research, for instance. Another is the formalism of osculating orbits, which is widely used to study celestial mechanics in the solar system, for example (Merritt 2013, Sect. 4.2; Poisson & Will 2014, Sect. 3.3.2). The correspondence between the parameters of both formalisms is given in Klioner & Kopeikin (1994). We use here the latter framework, whose base equations (i.e. the equivalent to Eq. (1)) form a system of first-order evolution equations for a set of six Keplerian orbital elements in time. These are called the osculating equations. Different sets of orbital elements are equivalent. We follow the formulation of Poisson & Will (2014) and pick (p, e, ι, Ω, ω, f). In order: semi-latus rectum, eccentricity, inclination, argument of the ascending node, argument of pericentre, and true anomaly. The last four elements are position and orientation angles (Fig. 1). As auxiliary orbital element, we also consider the semi-major axis,

(2)

Furthermore, we decompose the perturbative acceleration along its co-rotated Gaussian frame components,

(3)

where e_z is the unit vector in the direction of orbital angular momentum and λ completes the right-handed orthonormal basis (Fig. 1). With this, the osculating equations take the form

(4)

and we quote the explicit form of the right-hand side functions F^a from Poisson & Will (2014) in Appendix A. In general, ℛ, 𝒮, 𝒲 are functions of the μ^a themselves. Importantly, the functions F^a are linear in ℛ, 𝒮, 𝒲, and in particular,

such that the unperturbed case reduces to the Kepler problem, in which the first five orbital elements are constants of motion, and the true anomaly f strictly monotonically increases at the above rate. In the general perturbed case, all orbital elements may vary with time, describing a precessing or otherwise evolving orbit.

From the orbital elements, the components of the position r and velocity v = ṙ of the secondary body wit respect to the primary can be obtained by recognising that an Euler rotation around the angles Ω, ι, and ω + f rotates the fundamental frame (X, Y, Z) into the Gaussian frame (n, λ, e_z) (Fig. 1). A detailed derivation of the resulting transformations can be found in Poisson & Will (2014, Sect. 3.2.5), from where we quote the results in Eqs. (A.7) and (A.8).

2.3. Explicit form of a_1PN

With the star as test particle, the relativistic correction to the Newtonian equations of motion to 1PN order is given by

(5)

where v = |v| (Merritt 2013, Eq. (4.166); Poisson & Will 2014, Eq. (10.1)). Projecting this onto (n, λ, e_z), we obtain

(6)

(7)

and 𝒲_1PN = 0 for the corresponding Gaussian components (Merritt 2013, Eq. (4.207); Poisson & Will 2014, Sect. 10.1.3). From the vanishing of 𝒲_1PN together with Eqs. (A.3) and (A.4), it is immediately clear that a_1PN does not cause an out-of-plane orbital precession: the Schwarzschild precession is in-plane.

2.4. Explicit form of a_XM

We proceed with calculating a_XM for a Newtonian spherically symmetric extended mass distribution ρ(r). The easiest way to do so is via Gauss’ law in integral form, which in general reads

(8)

The integrals are taken over a volume V and its boundary surface ∂V. dV and dσ denote the respective volume and surface elements. a_⊥ is the outward-pointing normal component to ∂V of the acceleration created by the mass distribution ρ. Exploiting the symmetry and choosing V to be a sphere of radius r, we have a_⊥ = ℛ_XM(r) such that the left-hand side of Eq. (8) gives 4πr²ℛ_XM(r). On the right-hand side, we obtain . Equating both results, we obtain

(9)

for arbitrary spherically symmetric distributions ρ(r). Furthermore, due to spherical symmetry, we have 𝒮_XM = 𝒲_XM = 0. As in the 1PN case, from the vanishing of 𝒲_XM, it immediately follows that a spherically symmetric extended mass does not cause an out-of-plane precession. However, in contrast to the 1PN case, the vanishing of 𝒮_XM and Eq. (A.1) show that here, p is also a constant of motion.

In the following, we restrict our considerations to Plummer and power-law cusp profiles

(10)

with density and scale parameters ρ₀, r₀ (Fig. 2). (Note the different interpretations of ρ₀ for both profiles: maximum (i.e. central) density (Plummer) versus density at r₀ (cusp)). We limit the cusp exponent γ to the range (0,3) for an outward falloff and integrability. A prominent example is a Bahcall-Wolf cusp, which corresponds to

(11)

Fig. 2. Density profiles (blue, left horizontal axis) and enclosed masses (orange, right horizontal axis) for the Plummer (solid) and cusp (dashed) profiles (Eq. (10)). The profiles have a scale parameter r0 given by Eq. (14) and density parameters ρ0 given by Eq. (15). The cusp exponent γ is given by Eq. (11). rp and ra denote the pericentre and apocentre distance of S2, respectively. The profiles are chosen such that 0.1% of the galactic centre MBH mass reside within ra, i.e. they correspond to the current upper bound of GRAVITY Collaboration (2020).

(Bahcall & Wolf 1976; Merritt 2013, Sect. 5.5.2). This is the distribution to which a spherically symmetric nuclear star cluster should relax after a sufficiently long time. Concerning the galactic centre, the cusp profile has some theoretical and observational backing as well, although with slightly smaller exponents (Schödel et al. 2018; Baumgardt et al. 2018; Gallego-Cano et al. 2018). However, we have to keep in mind that these results are concerned with the distribution of the cluster as a whole³, while here we are concerned with scales of about the orbit of S2, that is, with the innermost part of the cluster in the immediate vicinity of the MBH. For these scales, a profile that (unlike a cusp) plateaus at a finite density at its centre may appear more plausible, which is why we set a slightly stronger focus on the Plummer profile in our analyses (Plummer 1911). Originally, it has been found by modelling the distributions of globular clusters.

For Eq. (10) we can perform the radial integral of Eq. (9) analytically and obtain

(12)

Finally, to yield Eq. (4) in closed form, we substitute r = p/(1 + e cos f).

Summarising the results of this section and Sect. 2.2, we have a net perturbative acceleration a_p = a_1PN + a_XM with Gaussian components

(13)

with the constituents given by Eqs. (6), (7), and (12). Hence we now have our osculating equations (Eq. (4)) of consideration in explicit and closed form (see also Eqs. (A.1)–(A.6)). This concludes the summary of our MBH + star + extended mass model. Some further general remarks are given in Appendix B.

3. Effects of the extended mass on the stellar orbit

In the following, we investigate the impact of an extended continuous mass distribution on the osculating orbital elements (Sect. 3.2) and on the observables (Sect. 3.3) of one full orbit of S2. In Sect. 3.1 we present the setup to this end, for which we restrict ourselves to specific mass distributions motivated by the recent upper dark mass bounds of GRAVITY Collaboration (2020). Section 3.4 is then concerned with varying the mass profile parameters in order to estimate detection thresholds. We then conclude this section with some comments about possible reservations concerning our theoretical analysis and the conclusions drawn from it in Sect. 3.5.

3.1. Setup

We integrate the osculating equations (Eq. (4)) (see also Eqs. (A.1)–(A.6)) over a full orbit f ∈ [0, 2π] for the following perturbed Kepler models (Eq. (1)):

(I)

(IIa)

(IIb)

(IIIa)

(IIIb)

The Gaussian components of these perturbative accelerations are given by the respective combinations of Eq. (13). In picking the profile parameters for the extended masses, we orient ourselves at the recent upper bounds obtained in GRAVITY Collaboration (2020) and choose a scale factor

(14)

together with a density parameter ρ₀ such that 0.1% of M_• resides within the apocentre distance of S2 (Fig. 2). (Here we use the apocentre distance of the initial osculating orbit, that is, p₀/(1 − e₀) with p₀ and e₀ given by Eq. (17)). Consequently,

(15)

In all cases, we choose the ∼2018 pericentre of S2 as initial point and time t₀ at which we prescribe the osculating orbital elements. For these, as well as for the MBH mass M_• and the distance to earth R₀, we pick the recent best-fit values of GRAVITY Collaboration (2020, Table E.1). Rounded to two digits in units of convenience, we thus have

(16)

(17)

(18)

where the angles are given in radians and where . The negative sign of the inclination angle arises because GRAVITY Collaboration (2020) used the opposite convention for the sense of rotation of the inclination (Appendix C).

3.2. Impact on the osculating orbital elements

Because it is the quantity that directly encodes the pericentre shift, we start our investigation with the evolution of the argument of pericentre ω. Figure 3 shows plots of Δω(f) = ω(f)−ω₀ together with ω′(f) for the models (I), (IIa), and (IIIa). The plots for the cusp models (IIb) and (IIIb) are very similar in terms of qualitative features to their Plummer counterparts. We thus do not show them separately and demonstrate our arguments at the example of the Plummer model.

Fig. 3. Δω(f) = ω(f)−ω0 together with ω′(f) for the models (I), (IIa), and (IIIa). The shaded blue region f ∈ (5π/6, 7π/6) marks the orbital section at which the extended mass drives the orbital change, while the 1PN correction is nearly ineffective. The opposite holds for the rest of the orbit. The vertical line within the shaded blue region marks the true anomaly in 2021.96, which is approximately the time of publication (GRAVITY Collaboration 2022).

Firstly, comparing Figs. 3a and b, we see that both the Schwarzschild and the mass precession have a secular effect, that is, they yield an accumulated pericentre shift over one orbit. Moreover, we see that the Schwarzschild precession is prograde (ω′> 0) while the mass precession is retrograde (ω′< 0), as is well known (Jiang & Lin 1985; Rubilar & Eckart 2001; Merritt 2013). Secondly, we emphasise as the key observation from these plots that while the Schwarzschild precession is nearly ineffective for f ∈ (5π/6, 7π/6) and predominantly active elsewhere on the orbit, the exact opposite is true for the mass precession; see the shaded blue region in the plots. In other words, the different perturbations are driving the orbital change over the course of distinct orbital sections, and we identify these sections spatially in Sect. 3.3 below. Consequently, Fig. 3c shows that despite its comparatively low mass, the mass precession even dominates the Schwarzschild precession in the shaded blue region, causing a temporary retrograde precession in the combined model (IIIa)⁴. On the one hand, it does not come as a surprise that the Schwarzschild precession is strongest around pericentre and that the mass precession is strongest around apocentre. After all, relativistic effects are stronger at shorter distances, and the Newtonian pull from the spherical mass distribution is greater the more of the distribution is between the star and the MBH. On the other hand, it is not trivial that these orbital sections are as clear cut and disjoint as is evident from the plots.

To further underline the latter point, we compare Δω(f) for models (I), (IIIa), and (IIIb) in one plot (Fig. 4). (We recall that all three models share the data of Eqs. (16)–(18)). From the plot it is evident that the discrepancy in the pericentre shift per orbit between models (IIIa)–(IIIb) with an extended mass and model (I) without one almost exclusively builds up in the shaded blue region. The curves for all three models lie almost on top of each other before this orbital section and stay at approximately constant in the separation after. This strongly emphasises the clear-cut spatial separation of both precession effects. It further demonstrates that when the observational data are limited to the f ∈ ( − 5π/6, 5π/6) section of the orbit, then the sensitivity to a dark mass will be very weak, and will substantially improve by additionally sampling the complementary orbital section f ∈ (5π/6, 7π/6). In the latter orbital section, data are directly sensitive do a dark mass. The plot also shows that the curves for models (IIIa) and (IIIb) do not differ significantly, which we can understand from Fig. 2: Despite the different density profiles, the enclosed masses within a certain distance do not differ significantly, and the enclosed mass is responsible for the Newtonian acceleration.

Fig. 4. Δω(f) = ω(f)−ω0 for models (I), (IIIa), and (IIIb). The discrepancy in pericentre shift per orbit between models (IIIa)–(IIIb) with extended mass and model (I) without it predominantly builds up in the shaded blue region. The vertical line within this region marks the true anomaly of model (I) in 2021.96, which is approximately the time of publication (GRAVITY Collaboration 2022).

We now turn to the other orbital elements. As discussed in Sects. 2.3 and 2.4, there is no out-of-plane precession in our cases, such that Ω and ι are constants of motion. Additionally, we established that p is unaffected by a spherically symmetric extended mass, such that we would not learn about the interference of the Schwarzschild and mass precessions from investigating it. We could investigate e, but because p is not altered by the extended mass, the semi-major axis a carries the same information, as shown by Eq. (2). We thus choose to plot a together with a′(f) in Fig. 5 as it is related to astrometry in a more direct manner than eccentricity. Again, the corresponding cusp plots show the same qualitative features, and we do not plot them separately. In contrast to the pericentre precession, there is no secular effect on a such that the semi-major axis returns to its initial value after a full orbit. Similarly as with the pericentre precession, the intermediate changes in a mostly go in the opposite direction for both effects, and they again predominantly occur in distinct orbital sections. The extended mass has again the strongest effect for f ∈ (5π/6, 7π/6), precisely where the 1PN perturbation is nearly ineffective, and it causes a temporary contraction of a with the minimum at apocentre f = π. Consequently, in the plot for the combined model (IIIa) (Fig. 5c), the effect of the extended mass can be identified as the little dip in the plateau around apocentre.

Fig. 5. Δa(f) = a(f)−a0 together with a′(f), analogous to Fig. 3.

Even though the osculating orbital elements are not observables, the plots of Figs. 3–5 already allow us to give a first estimate of the extended mass impact on astrometry. In Figs. 3b and c we identify a temporary retrograde pericentre shift of ∼ − 2.7 arcmin in the shaded blue orbital section. By trigonometry, this accounts for an astrometric shift near apocentre in the direction parallel to the minor axis of about the apocentre distance a(π)(1 + e(π)) times this angle. Dividing by R₀, this translates into ∼191 μas angular distance modulo the projection onto the sky. On the other hand, we estimate directly from Figs. 5b and c an astrometric impact along the major axis of ∼ − 0.8 au times (1 + e(π)), which translates into ∼ − 183 μas angular distance modulo the projection onto the sky. In the following section, we compare these estimates to those obtained directly from the astrometry plots.

3.3. Impact on astrometry and radial velocity

Motivated by the above results, we transform by Eqs. (A.7) and (A.8) from the osculating orbital elements to the position r and velocity v of the star in order to shed direct light on the observables RA, Dec, and RV. These are right ascension, declination, and radial velocity, where the last is defined as minus the Z-component of v.

Figure 6 shows a plot of the orbit projected onto the plane of the sky for model (IIIa) and f ∈ [0, 2π]. (The respective plots for models (I) and (IIIb) are indistinguishable from Fig. 6 by naked eye, up to the time labels Fig. 3). Although the gap at pericentre is too small to be visible by naked eye, the orbit is not closed due to the pericentre shift. Labels inside the orbit indicate the true anomaly f of the respective points in the trajectory, while labels outside the orbit give the corresponding time in years. S2 is drawn at its position in 2021.96 ca. the time of publication. It is apparent now that the above emphasised orbital section f ∈ (5π/6, 7π/6) marks about the apocentre half of the orbit (blue section), while the rest of the orbit marks the pericentre half (orange section). Comparing this to GRAVITY Collaboration (2020, Sect. 2 and Fig. 1) we see that of the dataset used to derive the current upper dark mass bound (0.1% of M_• enclosed by S2) the VLTI/GRAVITY astrometric data is limited to the pericentre half f ∈ ( − 5π/6, 5π/6). The apocentre half is well covered by both VLT/NACO astrometric and VLT/SINFONI spectroscopic data, the former however being less accurate than GRAVITY⁵. From the insights we gathered in Sect. 3.2 we thus expect that a tracking of S2 with GRAVITY (and ELT/MICADO; Davies et al. 2018, 2021) over the current apocentre half ending in ca. 2033 will allow to substantially improve the current upper dark mass bounds, or even to detect one if present and not much smaller than those considered here. Note that as evident from the figures, at the time of publication (ca. 2021.96) we already entered this domain (GRAVITY Collaboration 2022).

Fig. 6. Orbit projected onto the plane of the sky for model (IIIa) and f ∈ [0, 2π]. The position of S2 is marked in 2021.96, which is approximately the time of publication.

Having identified the orbital sections f ∈ ( − 5π/6, 5π/6) and f ∈ (5π/6, 7π/6) with the pericentre and apocentre halves respectively, we can now also gain a more intuitive understanding of the point already emphasised in Sect. 3.2, that the spatial separation between the orbital sections in which the extended mass is ineffective and effective is so clear cut. As established in Sect. 2.4, the extended mass acceleration is given by a_XM = ℛ_XM(r)n (Eqs. (3) and (12)). For it to divert the course of the star from its otherwise Keplerian orbit, thus causing a precession, two things need to hold: Firstly, its magnitude R_XM(r) needs to be sufficiently high, and secondly, its direction n needs to have a sufficiently large normal component to the star’s velocity vector. The former is given when there is a sufficient amount of extended mass between the star and the MBH, i.e when r is sufficiently large. The latter is given in particular around both pericentre and apocentre. However from Figs. 1 and 6 we see that both requirements are only well satisfied simultaneously in the apocentre half of the orbit.

In Fig. 7 we examine the functional dependence of the astrometric quantities on f. Figure 7a shows the difference in RA(f) between models (I) and (IIIa). The plot is overlaid with that of −ω′(f) for model (IIa) in order to emphasise the connection with the estimate which we extracted from Figs. 3b and c in Sect. 3.2⁶. The orbit in the sky is oriented such that the projected minor-axis (major-axis) is roughly parallel to the RA (Dec) axis and the extended mass causes a counterclockwise pericentre shift when the star is in the apocentre half (Figs. 1 and 6). This causes a shift of points on the orbit near apocentre in roughly RA-direction. The ∼191 μas estimate of Sect. 3.2 hence concerns the right ascension, and as a crude estimate is in agreement with the peak of ∼150 μas of RA for the Plummer model close to apocentre in Fig. 7a. Furthermore, the plot shows that indeed the ΔRA curve mimics that of −ω′(f). This is expected over the part of the curves where ΔRA increases because the mass precession drives this deviation. It appears non-trivial, however, that this deviation decreases again after the peak influence of the extended mass. This demonstrates that not only is the apocentre half of the orbit the section in which the effect of the extended mass on astrometry starts to become discernible, but also that the most pronounced deviations are limited to this orbital section. This shows through yet another facet that it is indeed very important to resolve the region f ∈ [5π/6, 7π/6] well in order to capture the dark mass signature (see also Sect. 4 for a support of this claim by means of a mock-data analysis).

Fig. 7. Dark mass impact on astrometry. (a) ΔRAI, IIIa(f) = RAIIIa(f)−RAI(f) together with . (b) Same for Dec together with . (c) together with the analogous quantity for model (IIIb); note the different scales on the left and right vertical axes. Subscripts indicate the models to which the respective quantities correspond (Eqs. (I)–(IIIb)). The vertical line within the shaded blue region marks the true anomaly of model (I) in 2021.96, which is approximately the time of publication.

Figure 7b shows the same as Fig. 7 for the declination. Here the plot for the difference in Dec(f) is overlaid with that of a′(f) as the orientation of the orbit in the sky allows us to draw a direct connection between semi-major axis and declination. Again, we observe the mimicking of the two curves and a maximum deviation of ∼125 μas, in agreement with our crude estimate of ∼183 μas in Sect. 3.2. The deviation here also begins at the beginning of the apocentre half, then peaks during it, and declines again afterwards. Again, the cusp plots corresponding to Figs. 7a and b show the same qualitative features and we do not show them separately. In terms of magnitude, the astrometric impact of the cusp is greater, however, as we show in the following.

In Fig. 7c we consider the absolute value of the astrometric deviation in the plane of the sky,

(19)

between models (I) and (IIIa). To demonstrate the similar impacts of the Plummer and cusp models, we overlay the plot with the respective quantity for models (I) and (IIIb). As expected from Figs. 7a and b, the effect for the Plummer model peaks with ∼164 μas at about apocentre, and it is predominantly limited to the apocentre half of the orbit. In particular, the effect only kicks in once the star enters the apocentre half. The same holds for our Bahcall-Wolf cusp with 0.1% of M_• within the apocentre of S2 (Eqs. (10), (11), (14), (15)), for which the peak is stronger, however, with ∼420 μas.

The observable left to analyse is radial velocity. In Fig. 8 we plot the difference in RV(f) between models (I) and (IIIa) together with the respective quantity for model (IIIb). Although in contrast to astrometry, the deviation does not decline after the main influence of the extended mass, the effect still begins in this orbital section. Consequently, the lesson of this plot is similar to that learned from Fig. 4 in Sect. 3.2: Given data limited to the pericentre half, the sampling of the apocentre half of the orbit will clearly improve the detection prospects of a dark extended mass.

Fig. 8. ΔRVI, IIIa(f) = RVIIIa(f)−RVI(f) (blue) together with the analogous quantity for model (IIIb) (orange). The subscripts indicate the correspondence to the models (I) and (IIIa)–(IIIb). The vertical line within the blue shaded region marks the true anomaly of model (I) in 2021.96—ca. the time of publication.

3.4. Detection thresholds

In Sect. 3.3 we assessed the impact of an extended mass on the observables at the basis of the residuals between orbits with and without an extended mass component. These residuals are a priori not observables because they involve two orbits, while in an actual observation, nature only confronts us with one orbit, that is, with the actual one. The question we address now is how we can relate these residuals to observables, and how we can properly read them as detection thresholds for an extended mass based on the accuracy of the available data. For this analysis, we restrict ourselves to the Plummer model and to astrometry. It is then also part of our investigation of Sect. 4 to examine the thresholds we obtained here theoretically from the perspective of a mock-data analysis.

Suppose that we are observing one full orbit in a situation with an extended mass. First we recall a key observation from the previous sections: The impact of the extended mass is negligible in the pericentre half. This insight allows us to identify this orbital section of the observed orbit with the corresponding orbital section of a hypothetical and thus unobservable orbit without extended mass. Consequently, this identification lifts residuals such as that for absolute astrometry ΔA into the domain of observables (see Eq. (19) and Fig. 7).

When we now restrict this to absolute astrometry, the question left to answer is the accuracy required for data in which parts of the orbit in order to resolve a certain astrometric residual ΔA. We tackle this question the other way around: Suppose that we have one orbit of (astrometrically) evenly distributed data given with astrometric accuracies of σ_p in the pericentre half f ∈ [ − 5π/6, 5π/6] and σ_a in the apocentre half f ∈ [5π/6, 7π/6]. In a simple ansatz, we view the residual ΔA as given by the difference between two uncorrelated measurements with uncertainties σ_p and σ_a, respectively, such that by error propagation, we have

(20)

for the uncertainty of ΔA. The proportionality constant κ we introduced by hand in order to a posteriori account for details that we a priori neglected in our ansatz, such as that the uncertainty with which the overall orbits are determined decreases with the number of data points. We gauge κ for the cases of our interest below. Importantly, once we have κ, then with Eq. (20) we arrive at a relation between the astrometric accuracies of our data points and the resulting accuracy with which the astrometric residual ΔA in Fig. 7c can be resolved. In particular, we have

(21)

and because ΔA is directly related to the extended mass parameters, we thus found a way to estimate which extended masses we should be able to detect with 1σ given our data.

We demonstrate this with the six cases of Fig. 9, each having the same number of (astrometrically) evenly distributed data points on one full orbit f ∈ [ − 5π/6, 7π/6], ensuring a balanced sampling of the pericentre and apocentre halves, and consequently, of the signatures of the Schwarzschild and mass precessions (Fig. 6). The cases differ, however, in the accuracy of these data in different orbital sections. The case of Fig. 9a thereby mimics the data from which the current 1σ upper bound of 0.1% of M_• (4261 M_⊙) enclosed by S2 has been deduced, with the 50 μas errors on the pericentre half representing the GRAVITY data, and the 350 μas errors on the apocentre half representing the NACO data (Gillessen et al. 2017; GRAVITY Collaboration 2020, Sect. 2 and Fig. 1). We can thus use this case to roughly gauge κ against observation by evaluating Eq. (20) at the detection threshold (Eq. (21)) given by the peak value max(ΔA) = 164 μas, which we found for model (IIIa) in Fig. 7c. This yields κ ≈ 1/2, and adopting this value for all other cases, we can now calculate their σ_ΔA from Eq. (20) and again identify them with the respective threshold peaks max(ΔA) by Eq. (21). The resulting values are given in the second column of Table 1.

Fig. 9. Schematics of six cases with data over one full orbit, differing by the astrometric accuracy. The bottom semi-circles represent the pericentre half of the orbit, and the top semi-circles show the apocentre half. Outside numbers mark the true anomaly, and inside numbers indicate the astrometric accuracies of the data in the orbital sections.

With this we can now for instance examine how low a Plummer mass of the same spatial scale of Eq. (14) we would still be able to detect with ∼1σ for each case. To do this, we vary the density parameter ρ₀ of model (IIIa) until the peak value max(ΔA) in Fig. 7c meets the thresholds we just calculated. The resulting values are given in Table 1 together with the corresponding number of solar masses within the apocentre of S2. As more data are taken with GRAVITY or MICADO in the coming years, the accuracy of the observational data will gradually improve following the analogous sequence of cases from Figs. 9a–c. Hence, from our estimates of Table 1, we can await to be able to improve the current upper bound by about a factor 0.7 until the next apocentre passage in approximately 2026, and by about a factor 0.2 until a full orbit is sampled with GRAVITY (and MICADO) at its current accuracy performance of 50 μas. During the publication process of the present work, a new upper bound has been obtained with the data gathered until August 2021, which in our scheme falls between the cases of Figs. 9a and b (see the 2021.96 indicators in the figures of this section). The result agrees excellently with our predictions (GRAVITY Collaboration 2022). The cases of Figs. 9d–f would become relevant if the accuracy of GRAVITY could be improved to its performance goal of 10 μas. However, based on Table 1 we would only await a further substantial improvement of the dark mass sensitivity due to this improvement once a full orbit were be sampled with that accuracy.

3.5. Possible reservations

We conclude this section with some comments about possible reservations for our above theoretical analysis and the conclusions drawn from it.

Firstly, all model orbits of our preceding investigations have the same initial osculating orbit at pericentre f = 0 (Eqs. (17) and (18)), and the question may be whether a bias results from this particular choice that might endanger the generality of our results. If we had chosen apocentre f = π as the common initial location instead, then taking the example of Fig. 4, all three curves would indeed coincide in the middle of the apocentre half by construction, and curve I would be farthest apart from curves IIIa, b in the pericentre half. However, this does not change the fact that then also, the slopes of the curves with and without extended mass disagree in the shaded blue region, which cause the deviation in the first place. In this sense, the slopes are more important here in the interpretation than the absolute values, which is also why we overlaid most curves with respective to first derivatives (see for example Figs. 3b and 5b, in which the influence of the extended mass is isolated). Hence also for f₀ = π, as well as for any other f₀, we would have arrived at the same conclusions, that is, that the extended mass drives the orbital change almost exclusively in the apocentre half. We refer to Sect. 4 below to rout out any remaining doubts. There we double-check all our theoretically obtained results in a mock-data analysis, which is not subject to any bias such as the above.

Secondly, in this section, we investigated the functional form of different quantities with true anomaly f as opposed to time t. Because the times at which the star reaches a certain f do not exactly coincide between our models with and without extended mass, this means for the residuals of Figs. 7 and 8 that we compared the observables of the stars at slightly different times. The utility of our conclusions drawn from these plots for observation (in particular the detection threshold estimates of Table 1) might therefore be doubtful because observations are made in the time domain. We first note in this respect that the respective time differences are small (see the time labels in Figs. 3 and 5). Furthermore, the mock-data analysis of Sect. 4, which operates in the time domain, serves as an independent verification in this question as well.

Finally, we comment on how we translated astrometric measurement accuracies into obtainable 1σ detection thresholds in Sect. 3.4. This approach is both ad hoc and heuristic, and we gauged κ against an observational case that is only approximately represented by Fig. 9a. In particular, we neglected here the technical but important point of aligning the reference frames of the different instruments (Plewa et al. 2015; GRAVITY Collaboration 2020). Furthermore, we only took one apocentre half of NACO mock data in order to make this case directly comparable to the others, and the real data violate our assumption of having the same number of evenly distributed points in both orbital halves. An assessment of the accuracy of our theoretical detection thresholds of Table 1 is thus also part of the mock-data analysis in Sect. 4.

4. Mock data analysis

We now complement the theoretical results of Sect. 3 with a mock-data analysis. In Sect. 4.1 we lay out the common setup for the following investigations: Firstly, we assess in Sect. 4.2 the dark mass sensitivity of increasingly good astrometric datasets for one full orbit, and we assess our theoretically estimated detection thresholds of Table 1 (Sect. 3.4). Secondly, we investigate in Sect. 4.3 the dark mass sensitivity of data restricted to one half or three quarters of an orbit. Here we support in particular our theoretical claim that data limited to the pericentre half are not sensitive to a dark mass. In Sect. 4.4 we extend the latter analysis by two further cases with again a full orbit of data, but with flipped accuracies compared to the corresponding cases of Sect. 4.2. Finally, we discuss the parameter correlations in Sect. 4.5.

4.1. Setup

Throughout Sect. 4, we consider the 1PN + Plummer mass model (IIIa) with a scale parameter r₀ of 0.3 μas (Eqs. (10) and (14)) and we use our Python-based code osculating orbits in General Relativistic environments (OOGRE), which we developed while working on this paper. The code is private for now, but we have plans to make it public. OOGRE is a perturbed Kepler-orbit model code optimised for stars in the galactic centre. It is based on the formalism laid out in Sect. 2 with the following adjustments and additions:

Firstly, we now also took into account the Rømer delay, the motion of the Solar System, the special relativistic transverse Doppler shift, and the gravitational redshift into the model (Appendix D; Kopeikin & Ozernoy 1999; Grould et al. 2017). The component of the relative motion of the Solar System and the MBH parallel to the line of sight thereby introduces a Z-velocity parameter in the model. These additional effects would have unnecessarily complicated the discussion for the preceding theoretical investigations, but they are important in order to adequately simulate an observational setting.

Secondly, from the set of six initial osculating orbital elements (Eqs. (17) and (18)), we swapped the semi-latus rectum p₀ and the true anomaly f₀ for the orbital period P₀ and the time of the (approximately 2018) pericentre passage T₀⁷. We also switched to the opposite inclination angle convention (Appendix C) such that in total we have (P₀, T₀, e₀, i₀, Ω₀, ω₀) as our new set of initial osculating orbital elements. This also represents a complete set of Keplerian orbital elements. No information is lost. We can for instance recover p₀ via P₀ and e₀ using (2) together with Kepler’s third law of planetary motion. We made this change of parameters only to adapt OOGRE to our fitting code infrastructure, which interfaces in the same way to our (fully) general relativistic orbit model based on the GYOTO code (Vincent et al. 2011; Grould et al. 2016, 2017; Paumard et al. 2019). The equations integrated at the core of OOGRE remain the same as before (Eqs. (A.1)–(A.6)).

Finally, also for the sake of consistent interfacing, we considered these elements now as corresponding to the orbit that osculates at the apocentre passage at time T₀ − P₀/2 (approximately 2010). Based on this, we have a model (IIIa) that, in this section, is specified by the 10 undetermined parameters of the 11-parameter set

(22)

To create mock data, we chose the following specific model: For the initial osculating orbital elements, M_•, R₀, and we again selected the best-fit values of GRAVITY Collaboration (2020, Table 1). Except when noted otherwise, we used a density parameter ρ₀ of the value given in Eq. (15) such that the profile corresponds to the current 0.1% of M_• upper bound of GRAVITY Collaboration (2020). The exceptions are those cases of Sect. 4.2 by means of which we examine the detection thresholds of Table 1, as well as some cases of Sect. 4.3, in which we increase the extended mass density by a factor 100. Hence our model for creating the mock data is given by Eq. (22) with a choice for the extended mass density parameter ρ₀, and rounded to two digits,

(23)

(24)

(25)

For each case, we chose 114 mock observation dates, and for each date, we took one astrometric and one spectroscopic data point (RA, Dec, RV). The number 114 was chosen because in the dataset that was used to derive the current upper dark mass bound (0.1% of M_• enclosed by S2), there are 57 GRAVITY data points of S2 that are distributed over the pericentre half of the orbit (GRAVITY Collaboration 2020, Sect. 2, Fig. 1)⁸. An equally dense sampling of the apocentre half would yield 114 data points over one full orbit. For a chosen ρ₀, the model orbit was then sampled at the observation dates, and the mock data then follow by offsetting the sample points by errors drawn from normal distributions with standard deviations of 10 km s⁻¹ in RV and case-dependent standard deviations in astrometry. Finally, for each case, we considered an ensemble of 250 such mock datasets.

We then performed a Markov chain Monte Carlo (MCMC) fit of the ten-parameter model (22) (r₀ being fixed) to each mock dataset⁹. From Fig. 1, the model curves (RA(t), Dec(t), RV(t)) are given by the Y and X components of r(t) and the −Z component of v(t), respectively. They were calculated from Eqs. (A.7) and (A.8) with the osculating orbital elements resulting from the integration of Eqs. (A.1)–(A.6) for the models (22) with parameters out of a rectangular domain around Eqs. (23)–(25). The results are then discussed in particular by means of the statistics of the best-fit values and error estimates (as given by the median values and 1/2 times the difference between the 84% and 16% percentiles of the ensemble of 250 adequately burned-in MCMC chains) of the density parameter ρ₀ for each case. We also discuss parameter correlations at the hand of corner plot representations of the posterior distributions in the parameter space. In all cases in which the MCMC fits converge, the chains hit autocorrelation after 500 to 1500 steps with 100 walkers, after which we burned-in and continued to sample the central mode for another 2500 steps. This means that our statistics for each case is based on fits to 250 mock datasets, and parameter chains for each fit with 2500 × 100 mixing samples, that is, with the burn-in samples already discarded.

4.2. Full orbit of data with increasing sensitivity

In our first investigation, we fixed ρ₀ to the current upper bound value of Eq. (15) and determined how the sensitivity to this particular dark mass increased with the astrometric accuracy along the sequence of cases from Figs. 9a–f. The 114 data points are thereby roughly (astrometrically) evenly distributed over the orbit. Figure 10 shows the density parameter statistics for the fits to the 250 mock datasets of the case of Fig. 9a (as indicated in the top left corner of the first plot).

Fig. 10. Density parameter statistics of the fits to 250 mock datasets of the case of Fig. 9a with histograms of the best-fit parameters (a) and of the corresponding error estimates (b). The blue normal distributions have means and variances according to the data of the histograms. The orange normal distribution in (c) has a mean of Eq. (15) (i.e. the value the mock data are based on) and a variance equal to the mean of the distribution in (b). The latter distribution has marks for 1 and 2 one-sided standard deviations. The number of sigmas with which zero is rejected is shown in the top right.

Panel a shows a histogram of the best-fit values for ρ₀. Panel b shows the corresponding histogram of the MCMC error estimates. We note that the error distribution is much narrower (a difference of about two orders of magnitude in standard deviation) than the value distribution. Both histograms are overlaid with a normal distribution plotted in blue, whose means and variances are taken from the data of the respective histograms. Additionally, Fig. 10a also shows a normal distribution plotted in orange. Its mean equals Eq. (15), that is, the value on which the mock data are based, and its variance equals the mean of the error distribution (Fig. 10b). For large numbers of mock datasets, the two distributions should converge. In particular, the mean of the blue distribution should converge to the mean of the orange distribution. With our statistics from 250 mock datasets, these two means are already merely ∼9.02 × 10⁻¹² kg m³ or ∼0.078σ apart. We therefore chose to interpret our results at the hand of the orange distribution that was constructed in this way, although it would not make a difference for the accuracy required for our analysis if we chose to use the blue distribution. In this sense, we conclude that in the case of Fig. 9a with a ρ₀ given by Eq. (15), the null hypothesis ρ₀ = 0 of having no extended mass is rejected with about 1.5σ, as indicated in the top right corner of Fig. 10a. This is roughly the expected sensitivity because this case was chosen in coarse analogy to the data from which the current 1σ upper bound was determined by GRAVITY Collaboration (2020).

Repeating the above analysis for the remaining cases of Fig. 9 results in Fig. 11, where again the sketches in the top left corner of the plots indicate the respective cases, and the number in the top right corner represents the rejection of the null hypothesis ρ₀ = 0. The sensitivity of the data to the dark mass of the current upper bound increases with the cases from a–f, that is, with decreasing astrometric errors. Comparing Figs. 11a and b, we estimate that a continued monitoring of S2 with GRAVITY at its current performance until its next apocentre passage in 2026 will improve the sensitivity to this particular dark mass by 1σ. The gain achievable with data gathered until the time of publication of the present work (approximately 2021.96) should consequently lie somewhere between zero and 1σ. From Fig. 11c we infer that another 1.5σ are gained by monitoring until a full orbit of GRAVITY data is taken in 2033, even without the use of NACO data. At this point, we could think about limiting the dataset to GRAVITY in order to fully take advantage of the instrument’s superior systematics (see the discussion in Sect. 5.3). The remaining Figs. 11d and e give an impression of how the dark mass sensitivity would improve if GRAVITY would be able to lower its systematic uncertainties to the level of the statistical uncertainties. Notably, from comparing Figs. 11e and f, we see that the level of improvement would be particularly high if a full orbit were monitored with the increased accuracy.

Fig. 11. Density parameter statistics for all cases of Fig. 9 with the mock data created using ρ0 = 1.69 × 10−10 kg m3, corresponding to 0.1% of M• within the apocentre of S2 and to the current upper bound of GRAVITY Collaboration (2020).

Next, we examined the detection thresholds estimated in Sect. 3.4. For this we again used the cases of Fig. 9, but now with the density parameters of Table 1. Following the same procedure as above results in Fig. 12, from which we see that the estimated 1σ detection thresholds of Table 1 are accurate.

Fig. 12. Density parameter statistics for all cases of Fig. 9 with the mock data created using the density parameters corresponding to the detection thresholds of Table 1, which also lists the corresponding masses enclosed by S2.

4.3. Less than one full orbit of data

We now explore the dark mass sensitivity of data that are limited to half or three quarters of an orbit. The different cases are indicated in the sketches in the top left corner of the plots of Fig. 13, which are to be interpreted analogously to the cases of Fig. 9, just that now the data do not cover one full orbit and double lines indicate orbital sections in which the density of data points is doubled. Due to the latter, we have 114 data points in total in theses cases as well, which allows a fair comparison of the results to the result obtained with a full orbit of data. Furthermore, it allows us to investigate whether missing data from one orbital section can be made up for by adding more data points in the respective opposite orbital section. The key point we are investigating here, however, is the importance of data in different orbital sections for constraining an extended mass.

Fig. 13. Density parameter statistics for the cases with data limited to less than one full orbit of data, as indicated in the sketches in the top left corner. For data limited to the pericentre half, the MCMC chain for ρ0 does not converge, but spreads out in a random walk. This means that these data are not sensitive to a dark mass, while in contrast, data limited to the apocentre half are; (c), (d). Panels a and b: the density is increased by a factor 100 such that 10% of M• are enclosed by S2. Panels e and f: can be compared to Fig. 11c, showing that a doubled sampling of the respective orbital sections cannot make up for the data gap to one full orbit.

We start with the case of 114 mock observations in the pericentre half of the orbit, with 50 μas and 10 km s⁻¹ precision. Figure 13a shows that for our 0.1% of M_• extended mass of Eq. (15), the MCMC chain fails to converge (i.e. to hit autocorrelation) even after 20 000 steps with 100 walkers. As noted above, in all other cases, autocorrelation is reached already after only 500 to 1500 steps. The failure of the MCMC fit to converge matches our theoretical observation of Sect. 3 that data limited to the pericentre half are not sensitive to a dark mass. We elaborate further on the plausibility of this: Fig. 13a shows that the MCMC chain for the density parameter spreads out with an envelope of

(26)

where ρ₀ denotes the value used to create the mock data and n is the MCMC step. This is reminiscent of random walk behaviour, indicating that the walkers are sampling a flat posterior distribution, or (e.g.) the peak of a very broad Gaussian¹⁰. We strengthened our case further by increasing the extended mass by a factor 100, such that the mock data were created with ρ₀ = 1.69 × 10⁻⁸ kg m³ (i.e. 10% of M_• or ∼43 000 M_⊙ within the apocentre), and repeating the analysis. With this density, the retrograde mass precession over one full orbit could easily be observed by naked eye. Despite this, Fig. 13b reveals that the chain fails to converge here as well and instead spreads out with precisely the same envelope (Eq. (26)) as for the dark mass that is orders of magnitude lower. The fact that this picture did not change even after blowing up the extended mass by a factor 100 strongly indicates that the posterior is indeed flat in both cases. We conclude that data limited to a pericentre half indeed are not able to constrain a present dark mass, even if it is relatively high. For the analysis of GRAVITY Collaboration (2020), this means that although its accuracy is orders of magnitude better, the GRAVITY dataset could not have constrained the dark mass on its own, but only in combination with the NACO and/or SINFONI data, simply because by itself at the time, it was limited to the pericentre half (see also the discussion in Sect. 5.3).

We proceed to analyse the case of data limited to the apocentre half. The fits converge and yield the distributions of Figs. 13c and d for 0.1% and 10% of M_• extended masses, respectively. The data are strongly sensitive to the higher extended mass and reject ρ₀ = 0 with 30.6σ. The lower extended mass is also constrained, but compared to Fig. 11c, the constraints are much weaker than in the case of a full orbit of data with the same number of observations, that is, a null rejection of 0.2σ as opposed to 4.0σ. We conclude that while data limited to the apocentre half are indeed sensitive to a dark mass, they are not strongly sensitive, and the loss of information from the pericentre half cannot be made up for by a denser sampling. We can understand this in the light of our discussion in Sect. 3.4: If we do not have samples in the pericentre half, we lack information about the part of the trajectory that plays the role of the hypothetical orbit without an extended mass, and without this, we cannot measure deviations such as Eq. (19) via which we infer the dark mass.

Next, we investigate the dark mass sensitivity of data limited to the orbital half from pericentre to apocentre, which are created with a ρ₀ given by Eq. (15). By intuition, we expect these data to be slightly more sensitive to a dark mass than the previous data because they catche both the peak and a tail of the curve in Fig. 7c. Indeed, a comparison of Figs. 13c and e shows a null rejection of 0.3σ versus 0.2σ. Although only slightly more sensitive, this difference is significant because the standard deviations of the error statistics are sufficiently small (see the analogous situation in Fig. 10). By comparison to Fig. 11c, it is clear that the missing information from the opposite orbital half cannot be made up for by a denser sampling here either.

Finally, we study the case of three quarters of an orbit of data ending in the next apocentre, where the sampling is doubled in the last quarter. The ρ₀ at the base of the mock data is again given by (15). For this setting, we find a 1.5σ rejection of ρ₀ = 0 (Fig. 13f). This is a strong improvement in comparison to the preceding three cases, which only had half an orbit of data. A comparison to Fig. 11b shows, however, that completing the orbit with data points still constrains the dark mass more strongly than a doubled sampling in the opposite quarter, even when the astrometric accuracy in the completing quarter is much lower. For the current observations of S2 with GRAVITY, this suggests that at least a fraction of NACO data will have to be added to the pool until a full orbit has been sampled with GRAVITY (and MICADO) until about 2033 in order to improve the dark mass constraints. Post 2033, GRAVITY data may be used alone to fully take advantage of the instrument’s superior systematics (see the discussion in Sect. 5.3).

4.4. Flipped accuracies

In the previous subsection, we confirmed amongst other results that data limited to a pericentre half are not sensitive to a dark mass, while data limited to an apocentre half are. However, we also showed that the latter sensitivity is weak, such that a full orbit of data is still much more constraining than an apocentre half with the same number of observations. We understood this phenomenon in the light of our discussion of Sect. 3.4 by the lack of information to measure discrepancies such as Eq. (19) via which we detect the extended mass.

We now extend this investigation and study two further cases with ρ₀ given by Eq. (15), which are related to those of Figs. 9a and e by flipping the astrometric accuracies between the orbital halves. The results are shown in Fig. 14. Comparing Figs. 11a and 14a as well as Figs. 11e and 14b, we see that if we have a full orbit of data, with different astrometric accuracies in the pericentre and apocentre halves (the accuracies in RV remain 10 km s⁻¹ everywhere), then the cases with better data in the pericentre half constrain the dark mass more strongly than their mirror cases. The difference is larger between the cases with 10 μas and 50 μas errors, but it is present in both comparisons. This result may seem counter-intuitive at first glance because in Sect. 4.3 we clearly confirmed that the pericentre half by itself is not sensitive to a dark mass. Naively, we might therefore expect that the apocentre half is also the more constraining half when a full orbit of data is available. However, from our discussion in Sect. 3.4, we understand that data in the pericentre half are necessary to constrain the hypothetical orbit without an extended mass in order to measure discrepancies such as Eq. (19), by which we determine the dark mass. Although the apocentre half is the half that is directly sensitive to the dark mass density, the pericentre half therefore appears to place the more stringent constraints on the initial osculating orbital elements, which indirectly seems to be the dominant factor in the overall constraining. In Table 2 we list the means of the error distributions (i.e. of the distributions that are analogous to that of Fig. 10b) for the initial osculating orbital elements for the cases we just discussed. Indeed, the errors of the elements are greater by a factor 2 to 3 in the cases with larger errors in the apocentre half than in their mirror cases.

Fig. 14. Density parameter statistics for the two cases indicated by the sketches in the top left corner. A comparison of the resulting null rejections to those of the mirror cases Figs. 11a and e shows that given one full orbit of data, the data in the pericentre half play the stronger role in constraining the dark mass than the data in the apocentre half.

4.5. Correlations

We conclude our analysis with an investigation of parameter correlations at the hand of corner plot representations of the posterior distributions¹¹. Our posterior distributions show already well-known correlations such as those between M_• and R₀, ω₀ and Ω₀, and ω₀ and e₀ (e.g. Gillessen et al. 2017; GRAVITY Collaboration 2020). We do not discuss these here, but instead focus on the correlations of the extended mass density ρ₀ with other parameters, as well as on the case-dependent overall qualitative features of the corner plots.

Figure 15 shows five correlations of ρ₀ with other parameters at the hand of the respective two-dimensional projections of the posterior distribution. These five correlations are present in all our cases. The diagrams plotted here correspond to a fit to a mock dataset of the case of Fig. 11c. The corresponding full corner plot is shown in Fig. E.1. Evidently, ρ₀ very strongly correlates with the time of pericentre passage of the initial osculating orbit T₀. This is clear from the perspective that in order to fit the same data, a later initial time can to a certain extent be counteracted by a stronger retrograde in-plane precession, and thus by a higher dark mass density. In this context, it is also interesting to compare this correlation to that between T₀ and the parameter f_SP in GRAVITY Collaboration (2020, Appendix E). f_SP controls the strength of the Schwarzschild precession, which is prograde, and hence these parameters are strongly anti-correlated. The correlation between ρ₀ and P₀ can also be understood via this link between the precession effect and a time shift. Similarly, we can also interpret the correlation between ρ₀ and the argument of pericentre of the initial osculating orbit ω₀: A larger ω₀ shifts the pericentre in the prograde sense, and thus can be countered by a higher retrograde precession, and consequently by a larger ρ₀. The weaker correlations of ρ₀ with Ω₀ and ρ₀ with e₀ are carried over from the already well-known correlations of the respective latter parameters with ω₀.

Fig. 15. Five correlations between the dark mass density ρ0 and other parameters that are observed in all cases. The particular plots shown here correspond to a fit of the case of Fig. 11c as indicated by the sketch (see Fig. E.1 for the full corner plot).

In terms of the overall qualitative features of the corner plots, we do not see variations within each case, but there are variations between the cases. Within the cases with a full orbit of data, these variations are small, but there is a clear contrast to the corner plots for cases with data that are limited to orbital sections, for which almost all parameters correlate mutually strongly. In Appendix E we show three corner plots to demonstrate this.

5. Discussion

In the following, we conclude with a short summary of the scope of our work (Sect. 5.1) followed by a collection of our general results (Sect. 5.2) and the implications and findings of direct relevance to present and future observation (Sect. 5.3). Finally, we give a short outlook about the utility of our analysis in the wider scope of galactic centre science (Sect. 5.4).

5.1. Summary

We explored how different precession effects of a star orbiting an MBH can clearly be separated from each other despite their net interference on secular timescales. The key for this was to determine their distinct signatures (in particular with respect to their locations), which they inscribe in the orbit on timescales of a single period. We focused on separating the secular interference of the (prograde) Schwarzschild and (retrograde) mass precessions, the former being the lowest-order relativistic deviation from Keplerian motion, and the latter being caused by a (dark) continuous extended mass distribution around the black hole.

We performed three layers of analysis: Firstly, we examined the impact of both precessions on the osculating orbital elements (in particular, on the argument of pericentre ω) (Sect. 3.2). Secondly, we investigated the resulting impact on the observables, that is, on astrometry and radial velocity (Sect. 3.3). Thirdly, we fitted model orbits to mock datasets of a series of cases, which differ by their astrometric accuracies in different orbital sections (Sect. 4).

5.2. General results

In summary, based on the above analyses, we list our key findings below.

1. While the mass precession almost exclusively impacts the orbit in the apocentre half, the Schwarzschild precession almost exclusively impacts it in the pericentre half, allowing for a clear separation of the effects (Sects. 3 and 4).

2. Data that are limited to the pericentre half are not sensitive to a dark mass, while data limited to the apocentre half are sensitive to it, but only to a limited extent (Sects. 3 and 4).

3. A full orbit of data is required to substantially constrain a dark mass. We traced this back in part to the circumstance that the pericentre half data are required to properly constrain the base Keplerian osculating orbit elements (Sect. 4).

4. Based on a full orbit of astrometric and spectroscopic data, the astrometric component in the pericentre half plays the stronger role in constraining the dark mass than the astrometric component in the apocentre half (Sect. 4).

An important general implication of these results is that despite their secular interference, the Schwarzschild and mass precessions can clearly be separated and measured independently of each other, for example, in a fit of a model comprising both components to a full orbit dataset of sufficient precision. Moreover, while these results have been obtained based on the orbit of S2 and with a variety of specific mock datasets, we have no reason to assume that they would not also hold in greater generality.

5.3. Implications for observation

Beyond these general statements, we also gathered results that are directly relevant for current and future observations of S2 with VLTI/GRAVITY and with upcoming instruments such as GRAVITY+ (Eisenhauer 2019) and ELT/MICADO. We repeat and expand on these points below.

Firstly, we interpret the current observational dark mass 1σ upper bound of 0.1% of M_• (∼4000 M_⊙) within the S2 orbit of GRAVITY Collaboration (2020) in the light of our above results. Result 2 implies that despite its superior astrometric precision, the corresponding GRAVITY dataset could not have constrained the dark mass unilaterally because it was limited to the pericentre half of the orbit. Based on this and result 3, the combination with NACO and/or SINFONI data was thus crucial. Furthermore, in the context of this combined dataset, we know from result 4 that the GRAVITY data have played the dominant role with regard to astrometry in constraining the dark mass.

Secondly, from these conclusions together with results 1 and 2, it is clear that the dark mass sensitivity of this dataset will increase substantially as also the apocentre half becomes populated with GRAVITY data points in the coming years. At the time of publication of this work (approximately 2021.96), we have entered this domain (see the figures of Sect. 3 and GRAVITY Collaboration 2022). We quantified this predicted gain in sensitivity by estimating the 1σ detection thresholds that can be achieved with one full orbit of data of different astrometric accuracies (Table 1 and Figs. 9 and 12). However, our theoretical cases and mock datasets have primarily been designed to extract the general statements of results 1–4 and only secondarily to mimic past, present, and future observational datasets and systematics (see also the discussion in Sect. 3.5). With this caveat, our estimates are as follows: With GRAVITY data collected until the next apocentre passage (approximately in 2026), at the current performance of the instrument of 50 μas, it will be possible to push the 1σ upper bound down to ∼3000 M_⊙ enclosed by S2. The bound in reach until the time of publication consequently lies somewhere between this value and the current bound. Once a full orbit of GRAVITY data is gathered (approximately in 2033), the 1σ upper bound can be pushed down to ∼1000 M_⊙ within the S2 orbit, without adding any NACO data to the pool.

Despite the importance of a compound of GRAVITY plus NACO and/or SINFONI data until 2033 (see result 3), a post-2033 restriction to GRAVITY as the sole astrometric component would have certain advantages because we would then fully benefit from the edge that interferometry has in systematics compared to classical imaging with 10-m class single-dish telescopes in the task of the astrometric tracking of S2. Firstly, in particular in the crowded environment of the galactic centre, astrometry with classical imaging is more strongly confusion limited (e.g. Trippe et al. 2010). Secondly, classical imaging suffers from static and variable image distortions (Plewa et al. 2015; Service et al. 2016). Thirdly, GRAVITY measures the astrometric position of S2 as a direct offset from Sgr A*, and thus from the MBH, while classical imaging with a 10-m class telescope does not, and has to bootstrap through the radio reference frame (Plewa et al. 2015; Sakai et al. 2019). Overall, while interferometry is certainly not free from systematics, astrometry from classical imaging involves more steps, each of which adds complexity and systematic uncertainties, which need to be modelled with additional free parameters (GRAVITY Collaboration 2022, Sect. 2).

Finally, we also estimated the degree to which the dark mass sensitivity could be further increased if the accuracy of GRAVITY were improved to ±10 μas (see Table 1 and Figs. 11 and 12). We found that this would have a particularly large impact once again one full orbit of data is gathered with this accuracy because by then, also the (in combination) more important pericentre half (see result 4) would also be sampled with maximum precision. Then, the 1σ upper bound would be pushed down to merely ∼200 M_⊙. At this low order of magnitude, it is likely that the continuous extended mass model would have to be refined to reflect the possible graininess of a distribution that could be made up of faint stars, stellar remnants, black holes, and/or dark matter. New observational constraints on the latter could come from the GRAVITY+ upgrade as well as from ELT/MICADO starting in approximately 2026, as both instruments will be able to observe the galactic centre to significantly fainter magnitudes (see also GRAVITY Collaboration 2021).

5.4. Outlook

For the larger picture, we hope that our work will serve as blueprint for analogous studies of the interference and separability of other effects. This should be useful in particular considering the challenging task of detecting the MBH spin via measuring the Lense-Thirring precession, which is the next great milestone in the field. Due to the smallness of the Lense-Thirring precession, it is of vital importance to examine how it interferes with other effects, and to explore how it can be separated from them (Wex & Kopeikin 1999; Merritt 2013; Zhang et al. 2015; Yu et al. 2016; Grould et al. 2017; Waisberg et al. 2018; Qi et al. 2021).

Acknowledgments

We acknowledge the support of CNRS [PNCG, PNGRAM] and Paris Observatory [CS, PhyFOG].

References

Appendix A: Osculating equations

The equations below are given in full generality, such that m denotes the total mass of the two bodies. We recall that in the extreme mass ratio case of an MBH + star, the star can be treated as a test particle, such that formally, m = M_∙. The osculating equations (Eq. (4)) in the form of Poisson & Will (2014, Eqs. (3.64)–(3.66)) read

(A.1)

(A.2)

(A.3)

(A.4)

(A.5)

(A.6)

where a dot denotes a derivative with respect to time t, and where ℛ, 𝒮, 𝒲 denote the Gaussian components of the perturbative acceleration (see Eq. (3) and Fig. 1). The latter are in general functions of the orbital elements themselves.

As mentioned at the end of Sect. 2.2 and shown in detail in Poisson & Will (2014, Sect. 3.2.5), the transformation law from the osculating orbital elements to r and v results from an Euler rotation around Ω,ι and ω + f,

(A.7)

(A.8)

Appendix B: Remarks about the model

Sect. 2 presented our MBH + star + extended mass model in the form of a perturbed Kepler problem with a perturbative acceleration a_p = a_1PN + a_XM. We thus correct the Newtonian two-body problem relativistically to 1PN order while treating the extended mass in a purely Newtonian manner. We comment on the domain of applicability of such a model below.

Firstly, the truncation at 1PN order for the two-body interaction is legitimate as long as the star orbits the MBH at distances much larger than its Schwarzschild radius. In the case of S2, the pericentre is ∼ 1400 Schwarzschild radii away from the black hole, such that post-Newtonian theory is well applicable.

Secondly, a purely Newtonian treatment of the extended mass distribution is justified as long as the error made by neglecting relativistic corrections for the extended mass is small compared to the modelled effects at play. This is the case as long as the extended mass is much lower than the mass of the MBH because then its 1PN correction will be negligible compared to the 1PN correction of the MBH, which is the only other modelled effect at play. As mentioned in the introduction, current upper bounds lie at the 0.1% level of the MBH mass, and we can conclude that our model is adequate. (We argue here in analogy to Poisson & Will 2014, p. 481), who laid out the negligibility of PN effects stemming from third-body perturbations in the solar system.

Finally, we compared our model to that of Rubilar & Eckart (2001) and Takamori et al. (2020), who included the extended mass directly into the Newtonian + 1PN two-body terms by replacing the MBH mass by an enclosed mass function of r. This is an elegant way to automatically include 1PN corrections for the extended mass in the model. However, it is ad hoc in the sense that it violates a criterion under which the 1PN two-body term has been derived in the first place: namely that the masses are isolated and well separated (Merritt 2013, p. 132; Poisson & Will 2014, Sect. 9.3).

Appendix C: Orbital angle conventions

When we introduced the angles (Ω,ι,ω) that give the orientation of an orbit as shown in Fig. 1 we followed the convention of modern theory based standard textbooks such as Merritt (2013) or Poisson & Will (2014). In particular, in the observation-based literature, the alternative convention shown in Fig. C.1 is also frequently encountered (see also the literal description in Paumard et al. 2006, Appendix A). Both conventions agree with respect to Ω and ω, but their inclination angles are taken in opposite directions: In Fig. 1, ι is taken in the sense of a right-handed rotation around the line of node axes pointing towards the ascending node, while in Fig. C.1, i is taken in the left-handed sense. Consequently, both conventions yield the same orbits when ι = −i. Furthermore, in Fig. 1 (Fig. C.1) a star transitioning through the ascending node is moving towards (away from) Earth.

Fig. C.1. Alternative orbital angle convention. L denotes the orbital angular momentum.

Appendix D: Approximated Rømer delay

Fig. 1 shows that modulo the time R_0/c that the light travels from the MBH to Earth, the difference between the time of emission and the time of observation on Earth t_o is given by

(D.1)

This is the Rømer equation, which has to be solved implicitly for the with calculated by the chosen orbit model if the full effect is taken into account. In practice, this can increase the computational cost and might become a problem in particular in tasks such as MCMC fitting, in which model orbits are created more frequently. However, at least for S2, we obtain an excellent approximation by substituting by its truncated first-order Taylor series around t_o, resulting in the approximated Rømer equation

(D.2)

In contrast to Eq. (D.1), Eq. (D.2) can be solved explicitly for t_e, while the error is small. For S2, the maximum deviation in astrometry of one orbit calculated with the full versus the approximated Rømer delay accounts for merely 1.2 × 10⁻² μas, which is orders of magnitude smaller than the effects we consider in this paper.

This approximation, as we used it here, has already been used in previous work (GRAVITY Collaboration 2018b, Sect. A.7). Takamori et al. (2020) dropped the first-order term in Eq. (D.2), by which the maximum deviation in astrometry within one orbit of S2 accounts for ∼ 2 μas for S2. However, the strongly increased accuracy by including the first-order term comes at negligible additional computational cost because Ż is in any case calculated by the model via Eq. (A.8). Moreover, the accuracy might be increased even further by adding the second-order Taylor series term, which comes at negligible additional computational cost as well because the acceleration on the star is given a priori.

Appendix E: Corner plots

In Figs. E.1–E.3 we present the corner plot representations of the posterior distributions for the cases of Figs. 11c, 13c, and 13e. The plot for the cases of Fig. 13d with a density higher by 100 times looks qualitatively indistinguishable from its lower-mass counterpart. We note that in these plots, the inclination i follows the convention of Fig. C.1.

Fig. E.1. Corner plot representation of the posterior distribution of a fit to a mock dataset corresponding to the case of Fig. 11c, as indicated by the sketch in the top right corner.

Fig. E.2. Corner plot representation of the posterior distribution of a fit to a mock dataset corresponding to the case of Fig. 13c, as indicated by the sketch in the top right corner. The plot for the case of Fig. 13d (i.e. with a density higher by 100 times) looks qualitatively indistinguishable.

Fig. E.3. Corner plot representation of the posterior distribution of a fit to a mock dataset corresponding to the case of Fig. 13e, as indicated by the sketch in the top right corner.

All Tables

All Figures

	Fig. 1. Definition of the Euler angles (Ω, ι, ω) giving the orbital orientation and the true anomaly f that gives the current position within the orbit. The static fundamental frame (X, Y, Z) (blue) is adapted to the observables, and the co-rotated Gaussian frame (n, λ, ez) (orange) is adapted to the symmetry of the (perturbed) two-body problem. See Fig. C.1 for an alternative convention.
In the text

Fig. 2. Density profiles (blue, left horizontal axis) and enclosed masses (orange, right horizontal axis) for the Plummer (solid) and cusp (dashed) profiles (Eq. (10)). The profiles have a scale parameter r0 given by Eq. (14) and density parameters ρ0 given by Eq. (15). The cusp exponent γ is given by Eq. (11). rp and ra denote the pericentre and apocentre distance of S2, respectively. The profiles are chosen such that 0.1% of the galactic centre MBH mass reside within ra, i.e. they correspond to the current upper bound of GRAVITY Collaboration (2020).

In the text

Fig. 3. Δω(f) = ω(f)−ω0 together with ω′(f) for the models (I), (IIa), and (IIIa). The shaded blue region f ∈ (5π/6, 7π/6) marks the orbital section at which the extended mass drives the orbital change, while the 1PN correction is nearly ineffective. The opposite holds for the rest of the orbit. The vertical line within the shaded blue region marks the true anomaly in 2021.96, which is approximately the time of publication (GRAVITY Collaboration 2022).

In the text

	Fig. 4. Δω(f) = ω(f)−ω0 for models (I), (IIIa), and (IIIb). The discrepancy in pericentre shift per orbit between models (IIIa)–(IIIb) with extended mass and model (I) without it predominantly builds up in the shaded blue region. The vertical line within this region marks the true anomaly of model (I) in 2021.96, which is approximately the time of publication (GRAVITY Collaboration 2022).
In the text

	Fig. 5. Δa(f) = a(f)−a0 together with a′(f), analogous to Fig. 3.
In the text

	Fig. 6. Orbit projected onto the plane of the sky for model (IIIa) and f ∈ [0, 2π]. The position of S2 is marked in 2021.96, which is approximately the time of publication.
In the text

Fig. 7. Dark mass impact on astrometry. (a) ΔRAI, IIIa(f) = RAIIIa(f)−RAI(f) together with . (b) Same for Dec together with . (c) together with the analogous quantity for model (IIIb); note the different scales on the left and right vertical axes. Subscripts indicate the models to which the respective quantities correspond (Eqs. (I)–(IIIb)). The vertical line within the shaded blue region marks the true anomaly of model (I) in 2021.96, which is approximately the time of publication.

In the text

	Fig. 8. ΔRVI, IIIa(f) = RVIIIa(f)−RVI(f) (blue) together with the analogous quantity for model (IIIb) (orange). The subscripts indicate the correspondence to the models (I) and (IIIa)–(IIIb). The vertical line within the blue shaded region marks the true anomaly of model (I) in 2021.96—ca. the time of publication.
In the text

	Fig. 9. Schematics of six cases with data over one full orbit, differing by the astrometric accuracy. The bottom semi-circles represent the pericentre half of the orbit, and the top semi-circles show the apocentre half. Outside numbers mark the true anomaly, and inside numbers indicate the astrometric accuracies of the data in the orbital sections.
In the text

Fig. 10. Density parameter statistics of the fits to 250 mock datasets of the case of Fig. 9a with histograms of the best-fit parameters (a) and of the corresponding error estimates (b). The blue normal distributions have means and variances according to the data of the histograms. The orange normal distribution in (c) has a mean of Eq. (15) (i.e. the value the mock data are based on) and a variance equal to the mean of the distribution in (b). The latter distribution has marks for 1 and 2 one-sided standard deviations. The number of sigmas with which zero is rejected is shown in the top right.

In the text

	Fig. 11. Density parameter statistics for all cases of Fig. 9 with the mock data created using ρ0 = 1.69 × 10−10 kg m3, corresponding to 0.1% of M• within the apocentre of S2 and to the current upper bound of GRAVITY Collaboration (2020).
In the text

	Fig. 12. Density parameter statistics for all cases of Fig. 9 with the mock data created using the density parameters corresponding to the detection thresholds of Table 1, which also lists the corresponding masses enclosed by S2.
In the text

Fig. 13. Density parameter statistics for the cases with data limited to less than one full orbit of data, as indicated in the sketches in the top left corner. For data limited to the pericentre half, the MCMC chain for ρ0 does not converge, but spreads out in a random walk. This means that these data are not sensitive to a dark mass, while in contrast, data limited to the apocentre half are; (c), (d). Panels a and b: the density is increased by a factor 100 such that 10% of M• are enclosed by S2. Panels e and f: can be compared to Fig. 11c, showing that a doubled sampling of the respective orbital sections cannot make up for the data gap to one full orbit.

In the text

	Fig. 14. Density parameter statistics for the two cases indicated by the sketches in the top left corner. A comparison of the resulting null rejections to those of the mirror cases Figs. 11a and e shows that given one full orbit of data, the data in the pericentre half play the stronger role in constraining the dark mass than the data in the apocentre half.
In the text

	Fig. 15. Five correlations between the dark mass density ρ0 and other parameters that are observed in all cases. The particular plots shown here correspond to a fit of the case of Fig. 11c as indicated by the sketch (see Fig. E.1 for the full corner plot).
In the text

	Fig. C.1. Alternative orbital angle convention. L denotes the orbital angular momentum.
In the text

	Fig. E.1. Corner plot representation of the posterior distribution of a fit to a mock dataset corresponding to the case of Fig. 11c, as indicated by the sketch in the top right corner.
In the text

	Fig. E.2. Corner plot representation of the posterior distribution of a fit to a mock dataset corresponding to the case of Fig. 13c, as indicated by the sketch in the top right corner. The plot for the case of Fig. 13d (i.e. with a density higher by 100 times) looks qualitatively indistinguishable.
In the text

	Fig. E.3. Corner plot representation of the posterior distribution of a fit to a mock dataset corresponding to the case of Fig. 13e, as indicated by the sketch in the top right corner.
In the text