# Constraining a companion of the galactic center black hole, Sgr A

Clifford M Will, Smadar Naoz, Aurélien Hees, Alexandria Tucker, Eric Zhang, Tuan Do, Andrea Ghez

► To cite this version:

Clifford M Will, Smadar Naoz, Aurélien Hees, Alexandria Tucker, Eric Zhang, et al.. Constraining a companion of the galactic center black hole, Sgr A. The Astrophysical Journal, 2023, 959 (1), pp.58. 10.3847/1538-4357/ad09b3. hal-04186952

HAL Id: hal-04186952

https://hal.science/hal-04186952v1

Submitted on 5 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
OPEN ACCESS

<figure><img src="image_2.png" /></figure>

CrossMark

# Constraining a Companion of the Galactic Center Black Hole Sgr A*

Clifford M. Will¹,² ID, Smadar Naoz³,⁴ ID, Aurélien Hees⁵ ID, Alexandria Tucker¹,⁶,⁷ ID, Eric Zhang⁸ ID, Tuan Do³ ID, and Andrea Ghez³ ID

¹Department of Physics, University of Florida, Gainesville, FL 32611, USA

²Institut d'Astrophysique, Sorbonne Université, F-75014 Paris, France

³Department of Physics and Astronomy, University of California, Los Angeles, CA 90095, USA

⁴Mani L. Bhaumik Institute for Theoretical Physics, Department of Physics and Astronomy, UCLA, Los Angeles, CA 90095, USA

⁵SYRTE, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, LNE, 61 avenue de l'Observatoire, F-75014 Paris, France

⁶Department of Physics, University of Illinois, Urbana, IL 61801, USA

⁷Illinois Center for Advanced Studies of the Universe, Department of Physics, University of Illinois, Urbana, IL 61801, USA

⁸Department of Physics and Astronomy, University of California, Riverside, CA 92507, USA

Received 2023 August 22; revised 2023 November 1; accepted 2023 November 3; published 2023 December 6

## Abstract

We use 23 yr of astrometric and radial velocity data on the orbit of the star S0-2 to constrain a hypothetical intermediate-mass black hole orbiting the massive black hole Sgr A* at the Galactic center. The data place upper limits on variations of the orientation of the stellar orbit at levels between 0°02 and 0°07 per year. We use a combination of analytic estimates and full numerical integrations of the orbit of S0-2 in the presence of a black hole binary. For a companion intermediate-mass black hole outside the orbit of S0-2 (1020 au), we find that a companion black hole with mass $m_c$ between $10^3$ and $10^5 M_\odot$ is excluded, with a boundary behaving as $a_c \sim m_c^{1/3}$. For a companion with $a_c < 1020$ au, a black hole with mass between $10^3$ and $10^5 M_\odot$ is excluded, with $a_c \sim m_c^{-1/2}$. These bounds arise from quadrupolar perturbations of the orbit of S0-2. Significantly stronger bounds on an inner companion arise from the fact that the location of S0-2 is measured relative to the bright emission of Sgr A* and that separation is perturbed by the “wobble” of Sgr A* about the center of mass between it and the companion. The result is a set of bounds as small as $400 M_\odot$ at 200 au; the numerical simulations suggest a bound from these effects varying as $a_c \sim m_c^{-1}$. We compare and contrast our results with those from a recent analysis by the GRAVITY collaboration.

*Unified Astronomy Thesaurus concepts: Supermassive black holes (1663); Astrophysical black holes (98); Galactic center (565); Gravitation (661); Three-body problem (1695); Milky Way dynamics (1051)*

## 1. Introduction

Sagittarius A* (Sgr A*) is a compact, bright radio source at the center of the Milky Way. Recent technological advances, such as the advent of adaptive optics, have made it possible to observe stars orbiting this source. The results imply that this is the likely location of a supermassive black hole (SMBH) of about 4 million solar masses (e.g., Ghez et al. 2000, 2008; Gillessen et al. 2009), surrounded by a cluster of stars (e.g., Ghez et al. 2003; Gillessen et al. 2009; Lu et al. 2013). Combined infrared (e.g., Keck observations; Witzel et al. 2018), radio, and X-ray observations (e.g., JVLA and Chandra observations; Dibi et al. 2016; Capellupo et al. 2017) have revealed hot emission from gas near the event horizon of Sgr A*. Observations by the Event Horizon Telescope Collaboration have provided evidence for the “shadow” of the black hole (Event Horizon Telescope Collaboration et al. 2022). Thus, the proximity of the Milky Way’s Galactic center provides a unique laboratory for addressing issues in the fundamental physics of SMBHs, their impact on the central regions of galaxies, and their role in galaxy formation and evolution.

The hierarchical nature of the galaxy formation paradigm suggests that galaxy mergers may result in the formation of binaries of SMBHs (e.g., Di Matteo et al. 2005; Hopkins et al. 2006; Robertson et al. 2006; Callegari et al. 2009). While

observations of SMBH binaries are challenging, there exist several confirmed binary candidates with subparsec to hundreds of parsec separations (e.g., Sillanpaa et al. 1988; Rodriguez et al. 2006; Komossa et al. 2008; Bogdanović et al. 2009; Boroson & Lauer 2009; Dotti et al. 2009; Batcheldor et al. 2010; Deane et al. 2014; Liu et al. 2014, 2016; Li et al. 2016; Bansal et al. 2017; Kharb et al. 2017; Runnoe et al. 2017; Pesce et al. 2018). Additionally, observations of dual active galactic nuclei with kiloparsec-scale separations have been suggested as SMBH binary candidates (e.g., Komossa et al. 2003; Bianchi et al. 2008; Comerford et al. 2009; Green et al. 2010; Liu et al. 2010; Smith et al. 2010; Comerford et al. 2018; Stemo et al. 2020).

If Sgr A* is a member of a binary, could its companion be an intermediate-mass black hole (IMBH), that is, a black hole with a mass in the range of hundreds to thousands of solar masses? Recent observations by the LIGO/Virgo/Kagra collaboration have now confirmed the existence of $100 M_\odot$ black holes (e.g., GW190521; Abbott et al. 2020a, 2020b). Our Galactic center may harbor IMBHs as a result of a possible minor merger with a low-mass or dwarf galaxy or even with a globular cluster. Such a scenario was considered by Rashkov & Madau (2013), who suggested that if IMBHs serve as the seeds of SMBHs in the center of galaxies, hierarchical galaxy evolution could yield many IMBHs in our galaxy. Additionally, a combination of theoretical and observational arguments have led to speculation that IMBHs may exist in the central parsec of the Galaxy (e.g., Hansen & Milosavljević 2003; Maillard et al. 2004; Gürkan & Rasio 2005; Gualandris & Merritt 2009; Chen & Liu 2013; Fragione et al. 2020; Generozov & Madigan 2020;

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** A hierarchical three-body system consisting of the Sgr A*-IMBH binary and the star S0-2. The binary can be either inside or outside the orbit of the star.</figcaption></figure>

Naoz et al. 2020; Zheng et al. 2020; Rose et al. 2022, 2023; Zhang et al. 2023).

In an earlier paper (Naoz et al. 2020), we constrained the allowable parameter space of an IMBH at the center of our Galaxy using the 20+ yr of observations of the star S0-2, which orbits the SMBH Sgr A* with an orbital period of 16 yr and an eccentricity of about 0.88. The recent closest approach of this star to Sgr A* (pericenter) has been used to test and confirm the prediction of general relativity (GR) for the relativistic redshift (e.g., GRAVITY Collaboration et al. 2018; Do et al. 2019) and the advance of the pericenter (GRAVITY Collaboration et al. 2020). The star S0-2 has been observed for more than two decades, and its orbit is sufficiently regular that, if there is a companion to Sgr A*, it is either quite close to the main black hole or well outside the orbit of S0-2, making this three-body system somewhat hierarchical in nature (Figure 1).

In this paper, we expand on Naoz et al. (2020) by developing improved analytical limits on a possible companion whose orbit is internal to that of S0-2 and limits on a companion external to the orbit of S0-2 and by obtaining bounds on a companion using direct numerical fits using publicly available data on the orbit of S0-2 from the UCLA Galactic Center Group (Do et al. 2019). In Section 2, we use the equations of motion for hierarchical triple systems expanded to quadrupole order together with a suitable averaging procedure to produce analytic estimates of the bounds on an IMBH. In Section 3, we describe the numerical methods by which we obtain bounds using the full array of data on S0-2. In Section 4, we review and update other bounds on a hypothetical companion. Section 5 makes concluding remarks.

## 2. Bounds on a Companion IMBH: Analytic Estimates

We consider a hierarchical triple system consisting of the massive black hole Sgr A* and a lighter black hole companion, with masses $m_*$ and $m_c$, respectively, and a star of mass $m_*$ such as S0-2 (see Figure 1). The two black holes may orbit each other within the orbit of the star, or the star may orbit the massive black hole in the presence of the lighter black hole orbiting outside the pair. We will denote these cases as the inner and outer companion cases, respectively. In the inner companion case, we assume that the ratio $a_c/a$ of the inner and outer semimajor axes is small or that the inner orbital period is short compared with the stellar orbital period. We treat the outer body, the star, as a massless test particle. It has no effect on the inner binary, but its orbit is perturbed by the varying

multipole moments of the inner binary's gravitational field. These have been denoted "inverse eccentric Kozai-Lidov" perturbations (e.g., Naoz et al. 2017; Zanardi et al. 2017). For the case of an outer companion, the ratio $a/a_c$ and the ratio of the stellar orbit period to the period of the companion black hole are assumed to be small. The resulting dynamical evolution from this evolution is known as the "eccentric Kozai-Lidov" (e.g., Kozai 1962; Lidov 1962; Naoz 2016).

We describe the three-body system in a coordinate system whose z-axis is parallel to the system's total angular momentum (see Figure 1); since the star has negligible mass, this implies that the z-axis is perpendicular to the plane of the black hole binary, which becomes the x-y reference plane (also called the "invariable plane"). The stellar orbit is inclined by an angle $\iota$ to the reference plane, intersecting it along the "line of nodes" at an angle $\Omega$ relative to the reference x-axis. The equivalent line of nodes for the companion orbit makes an angle $\Omega_c = \Omega - \pi$ relative to the x-axis. The pericenter angle of the companion orbit is $\omega_c$ from the line of nodes, or $\Omega_c + \omega_c$ from the x-axis. The star's pericenter is an angle $\omega$ from the ascending node of the stellar orbit. Each orbit is characterized by semimajor axes $a_c$ and $a$ and eccentricities $e_c$ and $e$ as usual. For simplicity, we assume that the two black holes have zero spin.

The Newtonian equations of motion for the three-body system are given by

$$ \frac{d^2\mathbf{X}}{dt^2} = -\frac{G(m_* + m_*)}{R^2}\mathbf{N} - Gm_c\left(\frac{\mathbf{X} - \mathbf{x}}{|\mathbf{X} - \mathbf{x}|^3} + \frac{\mathbf{x}}{r^3}\right), \quad (1a) $$

$$ \frac{d^2\mathbf{x}}{dt^2} = -\frac{G(m_* + m_c)}{r^2}\mathbf{n} + Gm_*\left(\frac{\mathbf{X} - \mathbf{x}}{|\mathbf{X} - \mathbf{x}|^3} - \frac{\mathbf{X}}{R^3}\right), \quad (1b) $$

where $\mathbf{X}$ and $\mathbf{x}$ are the positions of S0-2 and the companion, respectively, relative to the SMBH, with $R = |\mathbf{X}|$, $r = |\mathbf{x}|$, $\mathbf{n} = \mathbf{x}/r$, and $N = X/R$. For the purpose of our analytic estimates, we will henceforth set $m_* = 0$; its mass of about 14 $M_\odot$ will be included in the numerical integrations to be discussed in Section 3.

However, here we must account for two observational subtleties. The first is that the observed astrometric position of S0-2 on the sky is defined using a reference frame attached to the bright emission at the location of Sgr A*; put differently, the observations are differential measurements between the two images. The overall orientation of this reference frame is tied to galactic masers (Sakai et al. 2019). Now, in the absence of a companion black hole, this reference frame is the standard inertial frame attached to the center of mass of the Sgr A* star system. But if we are to consider a hypothetical companion, we must allow for the ensuing "wobble" of Sgr A* relative to the new center of mass located somewhere between the two black holes. The most straightforward way to do this is to describe the orbit of S0-2 by a variable that aligns with the observations, namely,

$$ \mathbf{X} \equiv \mathbf{x}_* - \mathbf{x}_c. \quad (2) $$

This is, in fact, the variable defined in Equation (1a). The second is that the data on S0-2 also include radial velocity (RV) measurements, which are defined relative to a local standard of rest, that is, an inertial frame. However, in the numerical analysis of the data, to be discussed in Section 3, these measurements are consistently referred to the astrometric frame
tied to Sgr A* by the transformation $V \rightarrow V - (m_c/M)v$, where $M = m_\odot + m_c$. Thus we will work in the astrometric reference frame throughout.

Depending on whether the star is inside or outside the orbit of the companion, we will expand Equation (1a) in powers of the small ratio of the two semimajor axes, treat the unperturbed orbital solution in terms of osculating orbit elements as described above, and insert the perturbing accelerations into the so-called Lagrange planetary equations, which govern the variation of the orbit elements with time.

We wish to obtain the observable change of the orbit elements that characterize the orientation of the orbit of S0-2, namely, $i, \Omega$, and $\omega$, that would be induced by the companion black hole. Upper limits on these variations have been placed using data from the two observational groups; the eccentricity $e$ of S0-2's orbit has also been well measured, leading to a rough upper bound on its rate of variation over the observation period, but we will not incorporate this in our study.

## 2.1. Orbit Elements in the Reference Frame and on the Sky

The star's angular orbit elements $i, \Omega$, and $\omega$ discussed above are defined in the reference system of Figure 1 and are not the same as the observed inclination, ascending node, and pericenter angles ($i_{sky}, \Omega_{sky}, \omega_{sky}$), defined with respect to the line of sight ($e_{Z,sky}$) and a basis on the plane of the sky ($e_{X,sky}, e_{Y,sky}$). The transformation of the orbit elements between these two bases is quite complicated, particularly since we do not know the orientation of the orbital plane of the hypothetical companion black hole a priori. However, there exists a set of invariant quantities that give a direct link between the two types of orbit elements. These arise from $j$ and $r$, respectively, the angular momentum and Runge-Lenz unit vectors for the stellar orbit, defined in the system reference frame by (see, e.g., Poisson & Will 2014)

$$ 
\begin{aligned}
j &\equiv \sin i(\sin\Omega e_X - \cos\Omega e_Y) + \cos i e_Z, \\
r &\equiv (\cos\Omega \cos\omega - \cos i \sin\Omega \sin\omega) e_X \\
  &\quad + (\sin\Omega \cos\omega + \cos i \cos\Omega \sin\omega) e_X + \sin i \sin\omega e_Z.
\end{aligned}
\quad (3) 
$$

As these unit vectors are mutually orthogonal ($j \cdot r = 0$), the resulting 3 degrees of freedom uniquely define $i, \Omega$, and $\omega$, and together with the unit vector $l \equiv j \times r$ they define a basis of vectors for the stellar orbit that can be related by a suitable rotation to either the sky basis or the reference basis defined by the companion black hole orbit. However, the following scalar quantities, constructed from $dj/dt, dr/dt, dl/dt$, and the three basis vectors, are invariant under rotations (note that $j \cdot dj/dt = r \cdot dr/dt = l \cdot dl/dt = 0$):

$$ 
\begin{aligned}
\mathcal{I}_1 &= \frac{dr}{dt} \cdot j = -\frac{dj}{dt} \cdot r = \sin \omega \frac{d\iota}{dt} - \sin i \cos \omega \frac{d\Omega}{dt}, \\
\mathcal{I}_2 &= \frac{dj}{dt} \cdot l = -\frac{dl}{dt} \cdot j = -\cos \omega \frac{d\iota}{dt} - \sin i \sin \omega \frac{d\Omega}{dt}, \\
\mathcal{I}_3 &= \frac{dr}{dt} \cdot l = -\frac{dl}{dt} \cdot r = \frac{d\varpi}{dt},
\end{aligned}
\quad (4) 
$$

where the variable $\varpi$ is defined by the relation $d\varpi/dt = d\omega/dt + \cos i d\Omega/dt$. Thus, a measurement on the sky of the variables that appear in Equations (4) yields, via construction

of the three invariants, values of those invariants expressed in the reference system defined by the companion's orbit.

However, at present, apart from $d\omega/dt$ for the pericenter of S0-2 (GRAVITY Collaboration et al. 2020), which is dominated by the GR precession, no measurements of $d\Omega/dt$ or $d\iota/dt$ exist, only rms errors on their values, consistent with zero. So in order to bound a hypothetical companion black hole, our strategy will be to convert the uncertainties in the rates of change $d\Omega/dt$ and $d\iota/dt$ into bounds, along with the uncertainty in $d\omega/dt$ after subtracting the GR effect. For simplicity, we will use two invariants constructed from $\mathcal{I}_1, \mathcal{I}_2$, and $\mathcal{I}_3$, namely,

$$ 
\begin{aligned}
\left| \frac{dj}{dt} \right|^2 &= \mathcal{I}_1^2 + \mathcal{I}_2^2 = \left( \frac{d\iota}{dt} \right)^2 + \sin^2 i \left( \frac{d\Omega}{dt} \right)^2, \\
\left| \frac{d\varpi}{dt} \right|^2 &= \mathcal{I}_3^2.
\end{aligned}
\quad (5) 
$$

Because we have no a priori knowledge of the orientation of the companion black hole's orbit, it makes sense to marginalize or average these invariants over those orientations. This implies averaging them over the two-sphere parameterized by $i$ and $\omega_c$ and the circle parameterized by $\omega$.

A similar approach was done in our previous work (Naoz et al. 2020) for an inner companion, averaging over the orbits. As we highlight below, here we relax the double averaging approach and expand this to the outer companion as well.

## 2.2. Outer Companion

We first consider a companion black hole outside the orbit of S0-2 so that S0-2 and Sgr A* form the inner orbit of the hierarchical triple and the companion is the outer perturber. With $m_\odot = 0$ in Equation (1b), $x$ evolves as a Keplerian orbit. Expanding Equation (1a) in powers of $R/r$ through quadrupole order, we find the equation of motion for S0-2,

$$ \frac{d^2\mathbf{X}}{dt^2} = -\frac{Gm_*N}{R^2} + \frac{Gm_cR}{r^3}(3(\mathbf{N} \cdot \mathbf{n})\mathbf{n} - \mathbf{N}). \quad (6) $$

The first term in Equation (6) is the Keplerian acceleration of S0-2 in the potential of the SMBH. We define its osculating orbit elements in the usual manner, with $X = RN, R = a(1 - e^2)/(1 + e \cos F)$, and the unit vector $N$ described in the $x-y-z$ basis using orbit elements $i, \Omega$, and $\omega$, along with sines and cosines of the true anomaly $F$, which satisfies the equation $dF/dt = \sqrt{Gm_*a(1 - e^2)}/R^2$ (see Poisson & Will 2014). The second term is the conventional quadrupole perturbation due to the outer black hole. To obtain the evolution of the orbit elements, we drop the first Keplerian term in Equation (6) and treat the remaining term as a perturbation of the Keplerian orbit; we find the components of the perturbations along the radial unit vector $N$, perpendicular to that vector but in the orbital plane and perpendicular to the orbital plane, and plug them into the Lagrange planetary equations.

The conventional approach would be to carry out the double time average of the equations over the inner and outer orbits (also called the “secular approximation”). However, we are in a regime where the inner and outer orbital periods are not necessarily very different, making the secular approximation suspect. In addition, we are less interested in detailed equations for the long-term evolution of the orbit elements of S0-2 than in estimates for changes in its orientation over the orbit and a half corresponding to the actual observations.
Accordingly, our method will be as follows. Holding the phase of the companion fixed, we first integrate the planetary equations for S0-2's orbit from $-\pi$ to $2\pi$; this roughly corresponds to the ~24 yr of observation of S0-2 from the apocenter through two pericenters. We divide by 1.5 orbital periods to get a per-orbit rate of change. We then construct the two invariants $|d\mathbf{j}/dt|^2$ and $|d\mathbf{\varpi}/dt|^2$ and marginalize them over the unknown orbital phase of the companion and over the unknown inclination $i$ and the two unknown pericenter angles $\omega$ and $\omega_c$. The results are

$$ 
\begin{aligned}
\left| \frac{d\mathbf{j}}{d\tau} \right|^2 &= \zeta^2 \mathcal{A}_{\text{out}} \alpha^{-6}, \\
\left| \frac{d\mathbf{\varpi}}{d\tau} \right|^2 &= \zeta^2 \mathcal{B}_{\text{out}} \alpha^{-6},
\end{aligned}
\quad (7) $$

where $\zeta \equiv m_c/m_\bullet$, $\alpha \equiv a_c/a$, which in this case is greater than 1; $\tau$ is time measured in units of the stellar orbital period, $P = 2\pi(a^3/m_\bullet)^{1/2}$; and

$$ \mathcal{A}_{\text{out}} = \frac{(8 + 24e_c^2 + 3e_c^4)(81\pi^2(2 + 6e^2 + 17e^4) + 512e^2(1 - e^2))}{1080(1 - e_c^2)^{9/2}(1 - e^2)} $$

$$ \mathcal{B}_{\text{out}} = \frac{(8 + 24e_c^2 + 3e_c^4)(64(1 - 2e^2)^2 + 567\pi^2e^2(1 - e^2))}{270(1 - e_c^2)^{9/2}e^2} \quad (8) $$

where $e_c$ is the eccentricity of the companion black hole's orbit.

## 2.3. Inner Companion

We next consider a companion whose orbit is inside that of S0-2, so that $a_c < a$. Expanding Equation (1a) (with $m_\star = 0$) in powers of $r/R$ we obtain the equation of motion

$$ \begin{aligned} \frac{d^2\mathbf{X}}{dt^2} = & -\frac{GMN}{R^2} - \frac{Gm_c\mathbf{n}}{r^2} - \frac{Gm_c r}{R^3}(3(\mathbf{N} \cdot \mathbf{n})\mathbf{N} - \mathbf{n}) \\ & -\frac{3Gm_c r^2}{2R^4}(5(\mathbf{N} \cdot \mathbf{n})^2\mathbf{N} - 2(\mathbf{N} \cdot \mathbf{n})\mathbf{n} - \mathbf{N}). \end{aligned} \quad (9) $$

The first term in Equation (9) is the Keplerian acceleration of the star in the asymptotic field of the binary of mass $M$. The second term is the “fictitious” acceleration of S0-2 caused by the acceleration of the astrometric reference frame attached to Sgr A*, which “wobbles” around the center of mass of the inner binary, while the third term is a “dipole” term caused by the fact that the vector $\mathbf{X}$ no longer points toward the center of mass of the black hole binary. The final term is the conventional quadrupolar perturbation, except for the fact that the mass factor in Equation (9) is $m_c$ instead of the usual reduced-mass factor $\mu = m_c m_\bullet / M$.

Carrying out the same procedure as for the outer companion, we obtain

$$ 
\begin{aligned}
\left| \frac{d\mathbf{j}}{d\tau} \right|^2 &= \frac{\zeta^2}{(1 + \zeta)^2} (\mathcal{A}_{\text{in}}\alpha^4 + \mathcal{B}_{\text{in}}\alpha^2 + \mathcal{C}_{\text{in}}\alpha^{-1} + \mathcal{D}_{\text{in}}\alpha^{-4}), \\
\left| \frac{d\mathbf{\varpi}}{d\tau} \right|^2 &= \frac{\zeta^2}{(1 + \zeta)^2} (\mathcal{E}_{\text{in}}\alpha^4 + \mathcal{B}_{\text{in}}\alpha^2 - \mathcal{C}_{\text{in}}\alpha^{-1} + \mathcal{F}_{\text{in}}\alpha^{-4}),
\end{aligned}
\quad (10) $$

<table><caption>**Table 1**
Orbit Elements on the Sky of Sgr A*</caption><thead><tr><th>Orbit Element</th><th>Value</th><th>σ<sub>Stat</sub></th><th>σ<sub>Syst</sub></th><th>Bound on Variation</th></tr></thead><tbody><tr><td>e</td><td>0.8858</td><td>0.0004</td><td>2.8 × 10<sup>-5</sup></td><td>2.9 × 10<sup>-4</sup> yr<sup>-1</sup></td></tr><tr><td>i (deg)</td><td>133.82</td><td>0.18</td><td>0.13</td><td>0.02 deg yr<sup>-1</sup></td></tr><tr><td>ω (deg)</td><td>66.11</td><td>0.24</td><td>0.077</td><td>0.07 deg yr<sup>-1</sup></td></tr><tr><td>Ω (deg)</td><td>227.49</td><td>0.29</td><td>0.11</td><td>0.07 deg yr<sup>-1</sup></td></tr></tbody></table>

**Notes.** Orbits elements from Table 1 of Do et al. (2019) and bounds on their variations (from Hees et al. (2017) and Naoz et al. (2020). The quantities $\sigma_{\text{Stat}}$ and $\sigma_{\text{Syst}}$ denote the statistical and systematic errors on the elements.

where $\alpha \equiv a_c/a$, and

$$ \mathcal{A}_{\text{in}} \equiv \frac{(8 + 40e_c^2 + 15e_c^4)(81\pi^2 + 16e^2)}{540(1 - e^2)^4}, $$

$$ \mathcal{B}_{\text{in}} \equiv \frac{8(2 + 3e_c^2)}{27(1 - e^2)^2}, $$

$$ \mathcal{C}_{\text{in}} \equiv \frac{32}{27(1 - e^2)}, $$

$$ \mathcal{D}_{\text{in}} \equiv \frac{(2 + e_c^2)[16(1 - e^2) + 81\pi^2e^2]}{54(1 - e_c^2)^{5/2}(1 - e^2)}, $$

$$ \mathcal{E}_{\text{in}} \equiv \frac{(8 + 40e_c^2 + 15e_c^4)[16(1 + 2e^2)^2 + 243\pi^2e^2]}{540(1 - e^2)^4}, $$

$$ \mathcal{F}_{\text{in}} \equiv \frac{(2 + e_c^2)[81\pi^2(1 - e^2) + 16e^2]}{54(1 - e_c^2)^{5/2}e^2}. \quad (11) $$

The effects of the terms in the equation of motion (9) can be seen in these expressions, with the quadrupole term ($\sim r^2$) appearing (squared) in the $\alpha^4$ term, the dipole term ($\sim r$) appearing in the $\alpha^2$ term, the wobble effect ($\sim r^{-2}$) appearing in the $\alpha^{-4}$ term, and the $\alpha^{-1}$ term representing a cross term between the wobble and dipole effects.

However, the step where we integrated over the orbit of S0-2 while holding the companion fixed in its orbit, while perhaps not unreasonable when the orbital periods are comparable ($a_c \sim a$), becomes problematic when the period of the companion is much shorter than that of the star. In this situation, the wobble and dipole effects on the orbit elements actually average to zero over one orbit of the companion. We will discuss this in detail when we compare these estimates with the data in the next subsection.

## 2.4. Observational Constraints

With over 23 yr of astrometric measurements and 19 yr of RV measurements, the Galactic Center Group has measured the orbit elements of S0-2 with reasonable accuracy (Do et al. 2019) (see Table 1). The GRAVITY collaboration has measured the orbit elements of the star (denoted by them as S-2) with similar uncertainties (GRAVITY Collaboration et al. 2018). Apart from a published measurement of $d\omega/dt$, which agrees with the GR prediction (GRAVITY Collaboration et al. 2020), no significant change in the other orbit elements has been detected to date.

In addition, the publicly available data from Do et al. (2019) have made it possible to estimate upper limits on a linear drift for each of S0-2's orbital elements. The orbital fit methodology is described in the Supplementary Materials of Do et al. (2019). The parameters included in the orbital fit were the mass of Sgr
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Bounds on the mass and semimajor axis of a companion IMBH in an orbit with $e_c = 0.6$. Red: bounds from $|dj/dt|$. Blue: bounds from $|d\varpi/dt|$. Dotted lines denote bounds on an inner companion purely from quadrupole perturbations of $|dj/dt|$ and $|d\varpi/dt|$. The regions labeled “a,” “b,” and “c” are discussed in Section 2.4. The tan area denotes a gravitational-wave damping timescale for the companion shorter than $10^7$ yr.</figcaption></figure>

$A^*$, its distance and line-of-sight velocity, and its position and velocity on the plane of the sky, the six standard orbital elements for S0-2. In addition, a linear drift for each orbital parameter was included. Statistical tests for model selection based on Bayesian evidence (see Do et al. 2019) showed that no significant deviations from zero were measured. An estimate of the 95% upper limit on a linear drift of S0-2's orbital elements was derived from the posterior probability distribution of the fit combined with an estimate of the systematic uncertainty derived from a “jackknife” analysis at the level of the reference frame construction (see Boehle et al. 2016; Do et al. 2019; Sakai et al. 2019). As a result, a 95% confidence upper limit on $|d\omega_{\text{sky}}/dt|$ was reported in Hees et al. (2017), and bounds on $|d\Omega_{\text{sky}}/dt|$ and $|d\dot{i}_{\text{sky}}/dt|$ were reported in Naoz et al. (2020). These are summarized in Table 1. The analyses described above did not include the possibility of an IMBH companion; that will be the subject of the next section.

Using the observed sky-basis orbit elements for S0-2, $a = 1020$ au, $e = 0.886$, $i = 134^\circ$, and the estimates shown in Table 1, we construct the observed bounds on the invariants $|dj/dt|$ and $|d\varpi/dt|$. Combining these with Equations (7) and (10), we obtain the bounds plotted in Figure 2. We have not utilized the bound on variations in the eccentricity in part because, for an inner companion, the eccentricity is constant to quadrupole order, while for an outer companion, the bound is not significantly different from the bounds obtained using the angular invariants, so we have not displayed that bound in Figure 2. The tan area in Figure 2 denotes the region where gravitational-wave emission would have caused a companion to merge with Sgr A* within the ~10 million year age of S0-2.

The red and blue lines in Figure 2 labeled “outer companion” show the $a_c \sim m_c^{1/3}$ trend characteristic of quadrupolar

perturbations of the orbit of S0-2. Companions below and to the right of those curves are excluded. The lines labeled “inner companion” begin with the $a_c \sim m_c^{-1/2}$ dependence expected for quadrupole perturbations and are extended as dotted lines labeled “inner, quadrupole only.” The bounds indicated by the dotted lines are largely consistent with the bounds shown in our earlier paper (Naoz et al. 2020). The region labeled “c” between the solid and dotted lines is excluded because of quadrupolar perturbations. However, as $a_c$ decreases, the effect of the wobble of Sgr A* begins to dominate, and the curves bend over to display an $a_c \sim m_c^{1/2}$ dependence. One might be tempted to conclude that a large set of companions below the lines labeled “inner companion” are excluded. However, we must recall that these curves were obtained by holding the companion fixed in its orbit while integrating the perturbations over the 1.5 orbits of S0-2. In fact, the companion’s period is shorter than that of the star, and therefore we should have integrated over the companion’s orbit first, as called for in the secular approximation. In that case, the wobble and dipole contributions from Equation (9) integrate precisely to zero. In reality, over the 1.5 orbits of S0-2, the perturbations due to the wobble and dipole terms will not average precisely to zero but will be suppressed relative to what is implied by Equation (10), depending on the specific relation between the orbital periods. In other words, the bound labeled “inner companion” could be porous, such that, for smaller $a_c$, some companions might be allowed because their wobble and dipole perturbations are sufficiently suppressed. This corresponds to the region labeled “a.” On the other hand, for a given $a_c$, the wobble and dipole effects grow linearly with $m_c$, and thus, for higher-mass companions, the wobble and dipole perturbations will be larger and more likely to lead to exclusion, corresponding to the region labeled “b.”

This long discussion illustrates the difficulty of drawing firm analytic conclusions in regimes where the hierarchical assumption is only marginally true, particularly when the frequency of the perturbation is shorter than that of the orbit being investigated. In such a case, one must turn to full numerical integrations of the equations of motion in hopes of obtaining a truer picture. Those will be the subject of the next section.

### 3. Bounds on a Companion IMBH: Inference Using Galactic Center Data and Integration of the Three-Body Equations of Motion

In this section we carry out a full numerical integration of the equations of motion of S0-2 in the presence of a hypothetical IMBH companion and perform a Bayesian exploration of the parameter space using the publicly available data from UCLA’s Galactic Center Group. The model fitted to the data is new and will be described in detail in this section, while the methodology and data used in our analysis are identical to the ones from Do et al. (2019) and will be briefly summarized.

In this analysis, we consider the Newtonian motion of the three bodies (the SMBH, its possible companion, and the star S0-2), but we also include $1/c$ relativistic correction in the expression of the RVs (relativistic redshift). In total, the model depends on 20 free parameters:

1. The masses of the SMBH $m$ and its companion $m_c$.

2. Six parameters describing the initial conditions of the companion with respect to the SMBH. We use osculating
elements at the reference epoch J2000: semimajor axis $a_c$, eccentricity $e_c$, inclination $i_c$, argument of periastron $\omega_c$, longitude of ascending node $\Omega_c$, and mean anomaly at J2000 $m_{0,c}$.

3. Six parameters corresponding to the initial conditions of S0-2. We use osculating elements at the reference epoch J2000: orbital period $P$, eccentricity $e$, inclination $i$, argument of periastron $\omega$, longitude of ascending node $\Omega$, and time of closest approach $t_0$.

4. The distance $R_0$ between the solar system and the Galactic center.

5. Four parameters to parameterize a possible drift of the reference frame ($x_0, y_0, v_{x0}, v_{y0}$, see Do et al. 2019).

6. One offset $v_{z0}$ for the RVs.

In the code, we fix the mass of the star S0-2 $m_*$ to 13.6 $M_\odot$, the nominal value estimated by Habibi et al. (2017). The osculating elements from both the black hole companion and S0-2 are transformed into Cartesian positions and velocities at J2000 using regular Keplerian transformations. Note that for S0-2, we distinguish two cases: (1) if $a_c > 1000$ au, we use the mass of the SMBH to convert the osculating elements into Cartesian coordinates since in that case, S0-2 is orbiting around the SMBH and it is perturbed by the outer body; and (2) if $a_c < 1000$ au, the transformation between osculating elements and Cartesian coordinates is performed using the total mass of the binary system since in this case, S0-2 is orbiting around the center of mass of the binary system (see Equations (6) and (9)).

From the Cartesian coordinates and velocities at J2000, we integrate the Newtonian Equations (1) for the three-body system. We safely neglect the first post-Newtonian corrections to these equations of motion, considering that the data set used in this analysis is not sensitive to these; see the discussion in Do et al. (2019). We integrated these equations of motion from J2000 forward and backward in time in order to cover the full observational time span, that is, from 1995 to 2018.

From the results of the numerical integration, we compute both the astrometric and the RV observable. We take into account the Römer time delay, which is due to the fact that the speed of light is finite, and thus the signal from the star takes a certain amount of time to propagate through S0-2's orbit in the z-direction. To first order in $1/c$, this delay can be approximated by Do et al. (2019)

$$t_{\text{em}} = t_{\text{obs}} - \frac{Z(t_{\text{obs}})}{c}, \qquad (12)$$

where $t_{\text{obs}}$ is the epoch of observation, $t_{\text{em}}$ the epoch of emission of the light, and $Z(t)$ the third component of $X(t) = (X, Y, Z)$.

The astrometric observations are the relative sky position of S0-2 with respect to the SMBH:

$$\bar{X}(t_{\text{obs}}) = \frac{X(t_{\text{em}})}{R_0} + x_0 + v_{x0}(t_{\text{obs}} - t_{\text{J2000}}), \qquad (13a)$$

$$\bar{Y}(t_{\text{obs}}) = \frac{Y(t_{\text{em}})}{R_0} + y_0 + v_{y0}(t_{\text{obs}} - t_{\text{J2000}}), \qquad (13b)$$

where $(\bar{X}, \bar{Y})$ are the astrometric observables and $x_0, y_0, v_{x0}$, and $v_{y0}$ model a 2D offset and linear drift of the reference frame.

On the other hand, the RV is not defined with respect to the SMBH but rather to the center of mass of the system. This is due to the fact that the RV is defined with respect to the local standard of rest (Dehnen & Binney 1998; Bland-Hawthorn &

Gerhard 2016). More precisely, the RV observable is computed as

$$RV(t_{\text{obs}}) = V_z(t_{\text{em}}) - \frac{m_c}{m_* + m_c} v_z(t_{\text{em}}) + v_{z0} + v_{GR}, \qquad (14)$$

where $V_z = dZ/dt$ (the line-of-sight component of the velocity of S0-2 with respect to the SMBH), $v_z = dz/dt$ (the line-of-sight component of the velocity of the companion with respect to the SMBH), and $v_{z0}$ is a constant velocity offset that accounts for possible systematic effects in the RV measurement or in the correction to the velocity of the local standard of rest. The first term in this equation is the standard Newtonian velocity projected along the line of sight. The second term is a correction to take into account the fact that the measured RVs are expressed in the local standard of rest (i.e., the origin of the RV reference frame is the dynamical center of mass of the binary black hole system, and this term corrects for the motion of the SMBH with respect to the center of mass of the system). This contribution becomes nonnegligible only for heavy companions orbiting close to the SMBH (i.e., small $a_c$ and large $m_c$). Finally, the last term encompasses the first relativistic corrections, which include the gravitational redshift from both the SMBH and its companion and the transverse Doppler predicted by special relativity (see Do et al. 2019 and GRAVITY Collaboration et al. 2018).

The data set used in this analysis is the same as that used in Do et al. (2019) and in Naoz et al. (2020). It consists of 45 astrometric positional measurements (spanning 24 yr) and 115 RVs spanning 18 yr. We used astrometric measurements obtained from the W. M. Keck Observatory by using speckle imaging (a technique to overcome blurring from the atmosphere by taking very short exposures and combining the images with software) from 1995–2005 and adaptive optics imaging from 2005–2018. These measurements are expressed in the reference frame developed in Sakai et al. (2019) and Jia et al. (2019) and are publicly available in Do et al. (2019). In addition, we used RV obtained from six spectroscopic instruments: one from NIRSPEC (Near-Infrared Spectrograph) on Keck, six from NIRC2 (Near-Infrared Camera 2) on Keck, 54 from OSIRIS (OH-Suppressing Infra-Red Imaging Spectrograph) on Keck, nine from NIFS (Near-infrared Integral Field Spectrometer) on Gemini, four from IRCS (Infrared Camera and Spectrograph) on Subaru, and 41 from SINFONI (SINgle Faint Object Near-IR Investigation) on the Very Large Telescope. The Keck, Subaru, and Gemini data are presented in Do et al. (2019; see also Chu et al. 2018), while the Very Large Telescope data are reported in Gillessen et al. (2017).

In this analysis, we use Gaussian likelihoods for the RV and the astrometric measurements. The RVs are supposed to be independent and normally distributed. We use a Gaussian likelihood for the astrometric measurements, including correlations between the measurements. The covariance matrix for the astrometric measurement depends exponentially on the sky-projected distance between two measurements. It is parameterized by a correlation length $\lambda$ and a mixing parameter $p$, which are both fitted simultaneously with all other model parameters. A detailed discussion of the likelihood used can be found in Section 1.5.1 of the Supplementary Materials from Do et al. (2019). Finally, following the analysis from Do et al. (2019), we also fit for an offset for the NIRC2 RVs.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Bounds on a companion IMBH from numerical simulations. The gray scale represents the posterior probability density distribution for the mass of the SMBH companion ($m_c$) and its semimajor axis ($a_c$) marginalized over the other 21 fitted parameters and normalized to 1. The scale of shades at the right of the figure indicates the relative probability. White regions are excluded. The black dashed line corresponds to the 95% confidence area. The blue curves are the same as in Figure 2. The dotted green curve corresponds to a $a_c \sim 1/m_c$ behavior. The tan area is excluded (more precisely is not relevant) as a result of gravitational-wave damping for the companion.</figcaption></figure>

In total, in this analysis, we fit simultaneously for 23 parameters: 20 model parameters, one offset for the NIRC2 RVs, and two parameters to model the correlations between the astrometric measurements. We perform a Bayesian inference for model fitting, using nested sampling to estimate the posterior probability distribution via the MultiNest package (Feroz & Hobson 2008; Feroz et al. 2009). The resulting 2D posterior for the parameters ($a_c, m_c$) marginalized over all the other 21 parameters is presented in Figure 3. The overall shape of the confidence area is similar to the one found in Gualandris et al. (2010). We notice two regimes representing an outer companion and an inner companion, with a turning point corresponding to a semimajor axis of the same order of magnitude as that of S0-2. The upper part of the figure shows an exclusion region (in white) fully compatible with the analytic estimate of Section 2.2, as indicated by the blue curve in Figure 3, taken directly from Figure 2. It shows the $ac \sim m_c^{1/3}$ behavior expected from quadrupolar perturbations of the S0-2 orbit.

The lower part of the figure shows quite different behavior. In addition to excluding high-mass inner companions that would induce quadrupolar perturbations on the orbit of S0-2 (to the right of the line labeled “inner, quadrupole only”), the results also exclude companions well to the left of that line, corresponding to the region labeled “b” in Figure 2. Here quadrupolar perturbations of S0-2’s orbit are very small, and its orbit serves as a “fixed reference” for observing the “wobble” of Sgr A* induced by the companion, much as the wobble of stars relative to a fixed background served to discover the first exoplanets. The absence of such an effect in the data serves to exclude companions, for example, with masses as small as 400 $M_\odot$ at 200 au. For companions of lower mass, the wobble is too small to be detected, as depicted by solutions with viable companions in the lower left-hand corner of Figure 3.

Other than revealing the potential importance of the wobble effect, the analytic approach does not do a good job of characterizing the bounds in this region of parameter space in detail, largely because it involves integrations over the observation time, which tend to wash out the effect. By contrast, the numerical integrations incorporate the full time dependence of the wobble and dipole effects, including correlations with other effects. The dotted green line denotes an approximate $a_c \sim m_c^{-1}$ dependence of the bound in this regime, which would suggest the effect of the “dipole” term in Equation (9), but how robust this is remains to be seen.

## 4. Other Bounds on a Companion

In a recent paper, the GRAVITY collaboration used numerical simulations to constrain the possibility of a companion IMBH (GRAVITY Collaboration et al. 2023). They used 4 yr of astrometric data (2017–21) and 21 yr of spectroscopic data, while we used 23 and 18 yr worth, respectively. Both analyses used the relative separation between S0-2 and Sgr A* as the fundamental variable, and both used similar sets of fitted parameters. They plotted posterior density distributions for allowed companions, one plot for an inner companion and one for an outer companion. For ease of comparison, in Figure 4 we have reproduced the two main components of their Figure 1 with the left (right) panel corresponding to an inner (outer) companion (note their axes are in linear scale and 0"/125 corresponds to 1020 au). The shades of blue in the GRAVITY plots correspond to 39%, 86%, and 99% confidence (from dark to light).

In the left panel, we have overplotted the analytic bound (in red) from pure quadrupole perturbations from an inner companion and the stronger bound (in green) inferred from our numerical simulations, which we have suggested result from the effects of the wobble of Sgr A*. The gravitational-wave bound is also shown for reference. In a region where we find no candidate companions—above and to the right of both the green and red curves—the GRAVITY analysis seems to find significant numbers of solutions.

In the right panel of Figure 4, we have plotted in red the analytic bounds for both an outer companion and an inner companion with quadrupole perturbations. Between these two curves, we argue that no companions should exist (and our numerical integrations find none), while the GRAVITY analysis shows a number of candidate companions. Above the “outer companion” curve, both analyses agree on the presence of companions compatible with the observations. The GRAVITY analysis actually goes beyond our study, showing that such outer companions may have a destabilizing effect on the S-star cluster, thus providing additional potential constraints on their existence.

Additional bounds on a hypothetical companion have resulted from limits on the wobble of Sgr A* relative to the distant quasars and from studies of the effect of a companion on the distribution of inclinations of the S-star cluster. These primarily exclude high-mass companions ($>2000 M_\odot$) exterior to the orbit of S0-2. Figure 13 of Gualandris & Merritt (2009) presents a summary of those bounds. Zhang et al. (2023) placed bounds using a stability criterion for nonhierarchical triple systems, arguing that certain companions could induce changes in the semimajor axis of S0-2 of order unity within the lifetime of the star. The resulting bound is consistent with the excluded region in Figure 3, although somewhat weaker (i.e., to the
<figure><img src="image_7.png" /></figure>

**Figure 4.** Comparison between the results of this paper and an analysis by the GRAVITY Collaboration et al. (2023). The left panel corresponds to an inner companion, and the right panel to an outer companion. The shaded blue regions denote companion solutions with 39%, 86%, and 99% confidence, from dark to light, respectively. White regions are excluded. The dashed red curves show our analytic bounds (taken from Figure 2) for purely quadrupole perturbations. The dashed green curve (taken from Figure 3) is a rough fit to our numerical results for low-mass inner companions. Above and to the right of both curves in the left panel, we predict no companions. Between the two dashed red curves in the right panel, we also predict no companions. Background figure reproduced from GRAVITY Collaboration et al. (2023) via Creative Commons CC-BY 4.0 license.

right). Broderick et al. (2011) discussed bounds that could be achieved using millimeter very long baseline interferometry.

## 5. Conclusions

We have used astrometric and RV data on the orbit of the star S0-2 to constrain a hypothetical IMBH orbiting the massive black hole Sgr A* at the Galactic center. We employed a combination of analytic estimates and full numerical integrations of the orbit of S0-2 in the presence of a black hole binary. For companions with masses above $10^3 M_\odot$, we found a wedge-shaped region in the space of $a_c$ versus $m_c$ (in log scale) where companions are excluded because their quadrupolar perturbations would induce changes in the orientation of S0-2 larger than the observations allow. Our analytic estimates and numerical simulations were in agreement in this regime.

For lower-mass companions inside the orbit of S0-2, analytic estimates suggested that the wobble of Sgr A* about the center of mass of its orbit with the companion would be the main observable effect but did not give reliable exclusion curves, probably because the averaging methods obscured important short-timescale effects. However, the numerical simulations verified the importance of the wobble effect and excluded a significant region of the $a_c - m_c$ parameter space, down to masses as small as 400 $M_\odot$ at 200 au.

## Acknowledgments

C.M.W. is grateful for the hospitality of the Institut d'Astrophysique de Paris, where parts of this work were carried out, and for partial support from NSF PHY 19-09247 and 22-07681. S.N. acknowledges partial support from NASA ATP 80NSSC20K0505 and the NSF-AST 2206428 grant and

thanks Howard and Astrid Preston for their generous support. A.M.G. is supported by NSF AST-1909554, the Gordon and Betty Moore Foundation, and the Arthur Levine and Lauren Leichtman Chair in Astrophysics at UCLA. A.T. is supported by NSF PHY-2212983. We are grateful to Odele Straub for useful comments and for permission to use Figure 1 from the GRAVITY Collaboration et al. (2023).

## ORCID iDs

Clifford M. Will https://orcid.org/0000-0001-8209-0393

Smadar Naoz https://orcid.org/0000-0002-9802-9279

Aurélien Hees https://orcid.org/0000-0002-2186-644X

Alexandria Tucker https://orcid.org/0000-0002-6656-2933

Eric Zhang https://orcid.org/0000-0002-7611-8377

Tuan Do https://orcid.org/0000-0001-9554-6062

Andrea Ghez https://orcid.org/0000-0003-3230-5055

## References

Abbott, R., Abbott, T. D., Abraham, S., et al. 2020a, PhRvL, 125, 101102

Abbott, R., Abbott, T. D., Abraham, S., et al. 2020b, ApJL, 900, L13

Bansal, K., Taylor, G. B., Peck, A. B., Zavala, R. T., & Romani, R. W. 2017, ApJ, 843, 14

Batchelder, D., Robinson, A., Axon, D. J., Perlman, E. S., & Merritt, D. 2010, ApJL, 717, L6

Bianchi, S., Chiaberge, M., Piconcelli, E., Guainazzi, M., & Matt, G. 2008, MNRAS, 386, 105

Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529

Boehle, A., Ghez, A. M., Schödel, R., et al. 2016, ApJ, 830, 17

Bogdanović, T., Eracleous, M., & Sigurdsson, S. 2009, ApJ, 697, 288

Boroson, T. A., & Lauer, T. R. 2009, Natur, 458, 53

Broderick, A. E., Loeb, A., & Reid, M. J. 2011, ApJ, 735, 57

Callegari, S., Mayer, L., Kazantzidis, S., et al. 2009, ApJL, 696, L89

Capellupo, D. M., Haggard, D., Choux, N., et al. 2017, ApJ, 845, 35

Chen, X., & Liu, F. K. 2013, ApJ, 762, 95

Chu, D. S., et al. 2018, ApJ, 854, 12
Comerford, J. M., Griffith, R. L., Gerke, B. F., et al. 2009, ApJL, 702, L82

Comerford, J. M., Nevin, R., Stemo, A., et al. 2018, ApJ, 867, 66

Deane, R. P., Paragi, Z., Jarvis, M. J., et al. 2014, Natur, 511, 57

Dehnen, W., & Binney, J. J. 1998, MNRAS, 298, 387

Di Matteo, T., Springel, V., & Hernquist, L. 2005, Natur, 433, 604

Dibi, S., Markoff, S., Belmont, R., et al. 2016, MNRAS, 461, 552

Do, T., Hees, H., Ghez, A., et al. 2019, Sci, 365, 664

Dotti, M., Montuori, C., Decarli, R., et al. 2009, MNRAS, 398, L73

Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2022, ApJL, 930, L12

Feroz, F., & Hobson, M. P. 2008, MNRAS, 384, 449

Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601

Fragione, G., Loeb, A., Kremer, K., & Rasio, F. A. 2020, ApJ, 897, 46

Generozov, A., & Madigan, A.-M. 2020, ApJ, 896, 137

Ghez, A. M., Morris, M., Becklin, E. E., Tanner, A., & Kremenek, T. 2000, Natur, 407, 349

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

Ghez, A. M., Salim, S., Weinberg, N. N., et al. 2008, ApJ, 689, 1044

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

Gillessen, S., Plewa, P. M., Eisenhauer, F., et al. 2017, ApJ, 837, 30

GRAVITY Collaboration, Abuter, R., & Amorim, A. 2018, A&A, 615, L15

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

GRAVITY Collaboration, Straub, O., & Bauböck, M. 2023, A&A, 672, A63

Green, P. J., Myers, A. D., Barkhouse, W. A., et al. 2010, ApJ, 710, 1578

Gualandris, A., Gillessen, S., & Merritt, D. 2010, MNRAS, 409, 1146

Gualandris, A., & Merritt, D. 2009, ApJ, 705, 361

Gürkan, M. A., & Rasio, F. A. 2005, ApJ, 628, 236

Habibi, M., et al. 2017, ApJ, 847, 120

Hansen, B. M. S., & Milosavljević, M. 2003, ApJL, 593, L77

Hees, A., Do, T., Ghez, A. M., et al. 2017, PhRvL, 118, 211101

Hopkins, P. F., Hernquist, L., Cox, T. J., et al. 2006, ApJS, 163, 1

Jia, S., Lu, J. R., Sakai, S., et al. 2019, ApJ, 873, 9

Kharb, P., Lal, D. V., & Merritt, D. 2017, NatAs, 1, 727

Komossa, S., Burwitz, V., Hasinger, G., et al. 2003, ApJL, 582, L15

Komossa, S., Zhou, H., & Lu, H. 2008, ApJL, 678, L81

Kozai, Y. 1962, AJ, 67, 591

Li, Y.-R., Wang, J.-M., Ho, L. C., et al. 2016, ApJ, 822, 4

Lidov, M. L. 1962, P&SS, 9, 719

Liu, J., Eracleous, M., & Halpern, J. P. 2016, ApJ, 817, 42

Liu, X., Greene, J. E., Shen, Y., & Strauss, M. A. 2010, ApJL, 715, L30

Liu, X., Shen, Y., Bian, F., Loeb, A., & Tremaine, S. 2014, ApJ, 789, 140

Lu, J. R., Do, T., Ghez, A. M., et al. 2013, ApJ, 764, 155

Maillard, J. P., Paumard, T., Stolovy, S. R., & Rigaut, F. 2004, A&A, 423, 155

Naoz, S. 2016, ARA&A, 54, 441

Naoz, S., Li, G., Zanardi, M., de Elía, G. C., & Di Sisto, R. P. 2017, AJ, 154, 18

Naoz, S., Will, C. M., Ramirez-Ruiz, E., et al. 2020, ApJL, 888, L8

Pesce, D. W., Braatz, J. A., Condon, J. J., & Greene, J. E. 2018, ApJ, 863, 149

Poisson, E., & Will, C. M. 2014, Gravity: Newtonian, Post-Newtonian, Relativistic (Cambridge: Cambridge Univ. Press)

Rashkov, V., & Madau, P. 2014, ApJ, 780, 187

Robertson, B., Bullock, J. S., Cox, T. J., et al. 2006, ApJ, 645, 986

Rodriguez, C., Taylor, G. B., Zavala, R. T., et al. 2006, ApJ, 646, 49

Rose, S. C., Naoz, S., Sari, R., & Linial, I. 2022, ApJL, 929, L22

Rose, S. C., Naoz, S., Sari, R., & Linial, I. 2023, ApJ, 955, 30

Runnoe, J. C., Eracleous, M., Pennell, A., et al. 2017, MNRAS, 468, 1683

Sakai, S., Lu, J. R., Ghez, A., et al. 2019, ApJ, 873, 65

Sillanpaa, A., Haarala, S., Valtonen, M. J., Sundelius, B., & Byrd, G. G. 1988, ApJ, 325, 628

Smith, K. L., Shields, G. A., Bonning, E. W., McMullen, C. C., Rosario, D. J., & Salviander, S. 2010, ApJ, 716, 866

Stemo, A., Comerford, J. M., Barrows, R. S., et al. 2021, ApJ, 923, 36

Witzel, G., Martinez, G., Hora, J., et al. 2018, ApJ, 863, 15

Zanardi, M., de Elía, G. C., Di Sisto, R. P., et al. 2017, A&A, 605, A64

Zhang, E., Naoz, S., & Will, C. M. 2023, ApJ, 952, 103

Zheng, X., Lin, D. N. C., & Mao, S. 2020, ApJ, 905, 169