© The Authors 2023

1. Introduction

For over a quarter century, near-infrared (NIR) astronomy has been successfully used to observe the nearest galaxy center available to astronomers, that of the Milky Way. The Galactic Center (GC) was discovered to be an incredibly rich stellar environment (Genzel et al. 1994, 1997, 1996; Eckart & Genzel 1996; Ghez et al. 1998, 2003; Schödel et al. 2002), which allows it to serve as a laboratory for studying the fundamental astrophysical properties of environments governed by the potential of a central supermassive black hole (SMBH). One ongoing experiment involves obtaining measurements for the positions of the stars within the inner 2 arcsec surrounding Sgr A* using a combination of adaptive optics (AO) and Speckle data sets spanning from 1992 to 2019 paired with the enhanced resolution of the interferometric data collected using the GRAVITY instrument on the Very Large Telescope Interferometer (VLTI) from 2017 to 2022 (Schödel et al. 2002, 2003; Ghez et al. 2003, 2008; Eisenhauer et al. 2005; Gillessen et al. 2009, 2017; GRAVITY Collaboration 2018a, 2022). In total, orbital determinations have been performed for ∼50 stars in the GC, all of which contribute to the understanding of the GC environment (Schödel et al. 2002, 2003; Ghez et al. 2003, 2008; Eisenhauer et al. 2005; Gillessen et al. 2009, 2017; Meyer et al. 2012; Boehle et al. 2016). This work continues the monitoring study of stars in the GC by expanding the projected outer search radius to the 7 arcsec surrounding Sgr A*, a regime occupied by the “clockwise disk” (Paumard et al. 2006; Bartko et al. 2010; Yelda et al. 2014; von Fellenberg et al. 2022).

In particular, we utilize data collected from 2002 to 2019 using the combined Nasmyth Adaptive Optics System (NAOS)- Near-Infrared Imager and Spectrograph (CONICA) instrument better known as NACO and the Spectrograph for INtegral Field Observations in the Near Infrared (SINFONI) mounted on the European Southern Observatory’s Very Large Telescope (ESO VLT) to track stellar motions for ∼2000 stars. Section 2 provides an overview of the NACO AO data utilized to detect accelerations as well as the SINFONI datacubes that contain the stellar spectra. Section 3 details the source finding algorithm developed to track the motions of the stars within the field. Section 4 discusses the methods employed to detect the stellar accelerations. Section 5 describes two methods for deriving orbital motions for the stars studied. The first pertains to stars that have measurable radial velocities from the SINFONI datacubes, while the second introduces a method that uses Monte Carlo sampling over a flat radial velocity prior for stars that do not have measured radial velocities. Section 6 relates the results from this study to previous works that focused on the structure of the nuclear star cluster environment represented by the GC. Finally, a summary of results as well as suggestions for future observations that would further this study are provided in Sect. 7.

2. Data

This work employs data from two instruments, namely NACO and SINFONI, that have since been decommissioned from ESO VLT in Paranal, Chile. The NACO images used in this work utilize the 13.3 mas pixel scale. It is typical for a data set to contain approximately 2 hours of data, with a single-detector integration time of ∼15 s (Gillessen et al. 2009). After integrating a few times on one field of view (FOV) it would be changed to a new pointing. There were typically four pointings used in total, such that the central 4 arcsec appear in all frames (Gillessen et al. 2009). The reduction of the data had already been performed (see, e.g., Schödel et al. 2002; Bartko et al. 2009; Gillessen et al. 2009, 2017); however, the process generally followed, which includes sky subtraction and flat fielding, is outlined in Gillessen et al. (2009). The reduction also includes processes that separate sources, namely moderate deconvolution in the central 5 arcsec with the Lucy-Richardson algorithm (Gillessen et al. 2009). In total, 85 epochs were selected from the archive of NACO data based on a visual inspection looking for the best images from each year between 2002 to 2019, a list of which is included in Tables A.1 and A.2.

Also key to this work, SINFONI is a NIR (1.1–2.45 micron) spectrograph that operated from 2004 until 2019 (Eisenhauer et al. 2003; Yang et al. 2019). SINFONI is composed of two instruments, namely MACAO, which is an AO module, and the Spectrometer for Infrared Faint Field Imaging (SPIFFI; Eisenhauer et al. 2003; Bonnet et al. 2003). Again, data cubes utilized in this work have been previously analyzed (e.g., Eisenhauer et al. 2005; Paumard et al. 2006; Bartko et al. 2009; Pfuhl et al. 2011). The majority of measurements were taken from a combined cube featuring observations from 18–19 August 2004. In this case, the observations have undergone a correction to remove features from the Earth’s atmosphere. Additionally, this cube is useful because it has a FOV of 8 arcsec × 8 arcsec, meaning most of the stars studied in this work are contained in this single cube. For stars that were not included in the FOV covered by this cube, or if the spectra are too noisy, a tool developed by von Fellenberg et al. (2022) was used to identify cubes with smaller FOVs featuring the relevant source. A table containing all of the stacked exposures is provided in Table B.1.

3. The source finder

In order to create a source finder, we selected a base image that we used to define the list of sources that were tracked over the 85 NACO images. The image 20110330 was selected due to the fact that it is the best image closest to the mean epoch of the sample, which will reduce correlations in the free parameters later on. We also considered how to best mitigate the deconvolution effects that occur around bright sources and can be misidentified as sources by the finder. We chose to perform image multiplication using three of the best NACO images including the base epoch and one from each of the 2010 and 2012 observing runs. By multiplying these images together we minimize the number of “fake” sources incorrectly identified by the finder. The images multiplied together were selected because they are the closest in time to the base epoch, thus minimizing the movement of the sources between the multiplied epochs. After image multiplication was performed, we used the Python package DAOStarFinder (Stetson 1987) to compile a list of sources to search for in all of the epochs. In total we detect 1839 sources in the multiplied image outside the central 2 arcsec. The next step toward creating a source finder algorithm that is successful for multiple epochs is to establish an image stacking procedure. If the images are stacked on top of each other, accounting for slightly different pointings and telescope rotation, the pixel search grid can be limited to less than 2 pixels. This is a generous limit, since stars with separations greater than 2 arcsec move significantly slower than the innermost stars. At this separation, stars generally have proper motions corresponding to less than one pixel between the epochs studied.

To stack the images and account for offsets in pointings and telescope rotations between epochs, we utilized third-order polynomial transformations to translate between the different NACO images and project the previously determined location for each source into the new epoch. As such, each epoch requires its own unique set of polynomial coefficients to perform this transformation. In order to calculate these coefficients, we identified a set of stars that appear in all epochs. In this case, 91 reference stars are used, which are the same as those defined in Gillessen et al. (2009) and Plewa et al. (2015). The position of each reference star is determined by hand for the base image then searched for specifically at all other epochs using a simplified version of the finder. Since these are particularly bright stars in the field, the finder performs very well which meant we could use simpler stacking techniques that do not require the set of transformations being determined. In this case, simply ensuring the positions of two bright sources are shifted such that they stack on top of the locations determined for the sources in the previous epoch was sufficient to obtain the remaining reference star locations.

Once these pixel positions were known for each epoch, we could solve for polynomial coefficients that transform the pixel positions from one observation to another. In this case, two sets of third-order polynomials were generated for each epoch. The first set advances the positions by one epoch, while the other advances the positions by the number of epochs since the 20110330 epoch. This is done in order to minimize the number of “lost” stars between epochs, or the number of stars for which no position can be determined by the finder in the new epoch. Essentially, if the star was detected in the previous epoch, the transformation coefficients from that epoch were used to project the pixel position to the pixel coordinates of the following epoch. These then served as an initial guess for the finder to determine the stellar position in the new epoch. If the star was not detected in the previous epoch, however, then the position of the star in the 20110330 image was transformed with the second set of coefficients to the new pixel coordinate system and again served as the initial guess for the finder. Stars can be lost in interim epochs either due to being a faint source in poor observing conditions or due to a confusion event. It is not possible to simply relate all epochs only using the polynomials that transform from the 20110330 image, as this set of coefficients is less successful at providing a good initial guess for the finder as the time between epochs increases. This is due to the motions of the sources that cannot be accounted for a priori¹. As such, fewer stars from the original list of 1839 are lost for the epochs closest to 20110330 as this secondary set of transformations is often successful at finding the star if it is detectable. Within the nine epochs selected from 2011, on average the loss of stars between each epoch is 0.6%. Whereas, within the six epochs used from 2002–2003, the average loss of stars between each epoch is approximately 2%. This increase is due to the fact that as the number of years increases from this mean epoch, the transformation often does not provide a sufficiently good guess to locate lost stars between epochs due to their stellar motion over time. This highlights another reason to select the mean epoch of the sample as the base image. We also note that although stars that appear to be moving fast may be preferentially removed due to the design of the source finder, this will not bias the sample of detected physical accelerations. Since at these separations the average motion between epochs should be less than a pixel, the search box utilized by the source finder is sufficient to detect stars moving with physical motions. In this case, losses of stars with physical motions that occur are mostly because they are faint stars lacking strong contrast with the background of the NACO observation. Common losses also include the loss of the stars located on the edge of the FOV that may be lost between epochs due to slightly different pointings.

Once there was an initial guess for the pixel position of the star, this was fed into the finder. The finder requires the definition of a fitting box to perform a Gaussian fit to the data contained in the box. To account for the fact that the sources do not appear uniform in size due to photon spilling in the CCD detector, we used multiple side lengths for this fitting box in iteration. These side lengths typically ranged from 8 to 16 pixels to get a good fit. Additionally, the finder performs a fit at each pixel moving outward from the initial guess position by 0 to 2 pixels in every direction. A simplified schematic of the process is shown in Fig. 1.

Fig. 1. Illustration of the source finding algorithm. First, the position from Epoch(N) is projected to Epoch(N+1) using the polynomial transformations. The transformed position is then used as an initial guess for the finder. Each grid square in the zoomed-in panel, which shows a NACO source profile, corresponds to a Gaussian search box of some side length. The fitter tries fitting a Gaussian profile to the pixels in each of the surrounding grid squares. The Gaussian fit with the lowest χ2 (white circle) is then selected as the position for this star in the new epoch. For clarity, the search boxes appear here as a grid; however, the actual algorithm utilizes overlapping search boxes.

The fit was performed by LevMarLSQFitter, which utilizes the Levenberg-Marquardt algorithm and least squares statistic to return a best-fit result (The Astropy Developers 2021). The fit corresponding to the lowest χ² was then selected. We verified the new position of the stars by visual inspection of the resulting pixel positions plotted onto the NACO image of the new epoch. We find that the source finder is capable of consistently determining the location of the source in a new epoch as required.

4. Calculating accelerations

With the stellar positions determined over all epochs, we now require a transformation that relates the pixel positions to on-sky coordinates. As previously mentioned, we used a set of 91 reference stars to define the GC coordinate system to a level of 0.7 mas accuracy (Gillessen et al. 2009). Once pixel positions for all the reference stars were obtained, we could determine a polynomial best fit that could convert pixel positions to right ascension (RA) and declination (Dec.) for each reference star at the relevant epoch. This was done using a least-squares fit and defining cubic transformations for both RA and Dec. individually. Applying the transformations to the stars detected by the source finder yields their RA and Dec. coordinates as a function of time. We then fit the positions for the star over all epochs it was observed using a second-order polynomial to determine the proper motions of the stars in RA and Dec.

In order to define the significance and direction of the accelerations detected for these stars, it is beneficial to rotate the coordinate system such that the motion of the stars is decomposed into its radial and tangential components. Similarly to the previous polynomial fit, we fit a second-order polynomial using numpy.polyfit (The NumPy Developers 2008) to the motion as a function of time in this rotated coordinate system such that the results have the form

$$\begin{array}{cc}& \frac{1}{2}(a_{\mathrm{radial}}\pm \mathrm{\Delta }{a}_{\mathrm{radial}}){(t-t_{\mathrm{ref}})}^2\hfill \\ \hfill & +(v_{\mathrm{radial}}\pm \mathrm{\Delta }{v}_{\mathrm{radial}})(t-t_{\mathrm{ref}})+(p_{\mathrm{radial}}\pm \mathrm{\Delta }{p}_{\mathrm{radial}}),\hfill \end{array}$$(1)

where a_radial is the acceleration in the radial component, v_radial is the velocity in radial component, p_radial is the radial component position of the star and the relative Δ terms refer to their respective errors, which are determined from the diagonal elements of the covariance matrix. The use of the (t − t_ref) term is due to the fact that the epochs are measured in years since 2010 (which is the approximate mid-point of the epoch list). The results are of the same form for the motion in the tangential direction. We used the polynomial fits to calculate the tangential and radial accelerations for each star and determined the significance of these accelerations. The significance of the radial acceleration, as well as its size and sign, indicates whether the motion of the star is dictated by the central SMBH. With this, we could quickly determine which stars show significant accelerations toward Sgr A* as it is simply a case of examining the second-order term and its error. This acceleration term must be significant in the radial direction and consistent with zero in the tangential direction. The accelerations of stars measured to be significant in the tangential direction are an indicator of the confusion for these measurements, as there is no known body in the GC that would govern the stars to move with this type of significant tangential motion (GRAVITY Collaboration 2022).

Next, we defined a condition to determine when the radial acceleration is considered sufficiently significant. We did this by calculating the significance in both the radial and tangential directions for each source. First the tangential significance was calculated for each star from the following:

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{tangential}}=\frac{1}{\sqrt{{\chi }_r^2}}\frac{a_{\mathrm{tangential}}}{\mathrm{\Delta }{a}_{\mathrm{tangential}}},\end{array}$$(2)

where ${\chi }_r^2$ is the reduced chi-square of the second-order polynomial fit to the tangential component of the data, a_tangential is the coefficient from the second-order fit corresponding to the acceleration and Δ a_tangential is the error. A histogram of the tangential acceleration significance for all sources was fit with two Gaussian profiles as shown in Fig. 2. The solid red line shows the Gaussian that describes the data in full. However, we find that erroneous identifications and nonphysical accelerations dominate the wings of the histogram. For this reason, the second fit shown by the dashed red line excludes the wings and provides a more accurate Gaussian fit for the peak of the distribution that contains the relevant data. The standard deviation of this Gaussian distribution is measured to be ∼1.212 and corresponds to the width of the dashed fit. This standard deviation was used to normalize the radial significance, namely,

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{radial}}=\frac{1}{\sqrt{{\chi }_r^2}{\sigma }_{\mathrm{tangential}}}\frac{a_{\mathrm{radial}}}{\mathrm{\Delta }{a}_{\mathrm{radial}}},\end{array}$$(3)

Fig. 2. Histogram illustrating the distribution of tangential acceleration significance within the sample of stars examined. The histogram is normalized such that the area under the curve is equal to 1. The tangential significance for each star is calculated from Eq. (2). Two Gaussian fits to the histogram are plotted. The first includes the wings of the histogram when finding a best fit (solid curve) that is dominated by sources with nonphysical accelerations, while the second excludes the data located in the wings of the histogram (dashed curve). Therefore, the standard deviation for the curve, which excludes the wings, was calculated and used to normalize the radial acceleration significance. In this case, the standard deviation is found to be ∼1.212.

where σ_tangential now refers to the standard deviation of the distribution shown in Fig. 2. The result is an adjusted significance that we used to define a list of stars that have a radial acceleration due to Sgr A*.

Additionally, we excluded sources that exhibit an apparently nonphysical acceleration. We did this by calculating the maximally possible acceleration that is reached assuming z = 0, where z is the line-of-sight distance to the star:

$$\begin{array}{c}\hfill a_{(z=0)}=-\frac{G{M}_{\bullet }}{r_{2\mathrm{D}}^2}.\end{array}$$(4)

Here, G is the gravitational constant, M_• is the mass of Sgr A* (here taken to be 4.3 × 10⁶ M_⊙GRAVITY Collaboration 2022) and r_2D is the projected distance of the source from Sgr A*. If the second-order fit produced a coefficient that corresponds to a greater acceleration than this maximum, the acceleration is nonphysical and this star was rejected from the list. In this case, the number of stars exhibiting nonphysical accelerations increases as we look in rings of increasing separations from Sgr A*. This demonstrates that the reasons for these nonphysical accelerations are largely due to the source finder misidentifying a star after a confusion event. Additionally, as the projected distance from Sgr A* increases there are increased gaps between epochs where stars are detected due to the source locations not being recovered by the source finder. The stars further out are often fainter and depending on the pointing of the telescope, may be excluded from the FOV for some epochs. These are challenging conditions for the source finder that can lead to misidentifications that result in accelerations that exceed physical limits, sometimes by full orders of magnitude. In total, of the 677 stars that have negative accelerations, only 168 stars (approximately 25%) are not excluded due to nonphysical accelerations. Stars that exhibit physical accelerations are not likely to have their accelerations falsely increased by a similar effect, as inspection of the source positions in RA and Dec. would clearly reveal if source misidentifications have occurred.

Extending the use of the Lucy-Richardson algorithm to be applied to faint stars beyond separations of 5 arcsec could improve the contrast between the sources and the background, thus improving the positions determined by the source finder, and in turn, reducing the number of stars featuring nonphysical accelerations.

We must now define our final list of stars with significant and physical accelerations due to Sgr A*. For this to be the case, both the acceleration and the significance calculated in the radial direction must be negative. For this reason, we set the condition such that σ_radial ≤ −3. For stars with a physical acceleration that does not meet this condition, continued monitoring may help improve the significance of their radial acceleration in the future.

With all of these conditions considered, in total, 54 stars within 2–7 arcsec from Sgr A* are found to exhibit significant radial accelerations. These stars are key to understanding the dynamics resulting from the potential of Sgr A* and are, therefore, the focus of the remainder of this work. These stars are summarized in Table 1.

With this final list of stars that experience a significant radial acceleration, we can perform further analysis to determine more information about their motion within the GC. Hereafter, “radial acceleration” will be referred to as “astrometric acceleration” to avoid confusion with the traditional “radial velocity” term referring to the v_z component of the motion. The polynomial fits in both RA and Dec. for all 54 of these stars are provided in Tables C.1 and C.2.

5. Determining stellar orbits

All 54 of the stars featuring significant astrometric accelerations in Table 1 were then classified as “5D constrained” stars. This means that, of the six phase-space coordinates required to determine the six orbital elements, five were known. The known elements included the on-sky position of the star (x, y) and the proper motion of the star (v_x, v_y). Additionally, by determining the acceleration from fitting the on-sky orbit trace in Sect. 4, we could determine the absolute value of the z coordinate for each star, given by

$$\begin{array}{c}\hfill |z|=\sqrt{|{\left(\frac{-G{M}_{\bullet }{r}_{2\mathrm{D}}}{a_{2\mathrm{D}}}\right)}^{2/3}-r_{2\mathrm{D}}^2|},\end{array}$$(5)

where a_2D is the 2D acceleration term, $a_{2\mathrm{D}}=\sqrt{a_x^2+a_{\mathit{y}}^2}$, and similarly, r_2D is the 2D position vector.

The missing phase-space coordinate was then the velocity in the z coordinate, or the radial velocity of the star. We used two methods to extract radial velocity information for the 54 stars demonstrating significant acceleration toward Sgr A* and determine their orbital solutions. Section 5.1 details the first method involving the determination of radial velocities from spectroscopic observations. For stars where we could not obtain a reliable spectrum, we applied Monte Carlo sampling in radial velocity to determine a range of physically viable orbital solutions. This method is discussed in Sect. 5.2.

5.1. Measured radial velocities from SINFONI data

To obtain spectra for the stars listed in Table 1, we used datacubes from the ESO’s decommissioned SINFONI instrument. Stars were identified in the IR image via cross-comparison to NACO images with the nearest epoch to the SINFONI datacube. Once identified, we extracted their spectra. From the stellar spectra, we could then classify the stars as either spectroscopically old or young depending on the characteristic features present in the spectrum of the star. This process is explained in detail by von Fellenberg et al. (2022) but can be summarized as follows: the Bracket-γ (Br-γ) line located at ∼2.166 μm or the, often less discernible, helium absorption line at ∼2.113 μm is considered as a criterion for young stars and the CO band heads that start at ∼2.3 μm for old stars.

In some cases the noise level in the spectrum was too high and we could not confidently extract a line from the spectrum. We also utilized line maps for stars that had faint absorption lines to confirm their origin was consistent with the location of the source and the expected full width at half maximum source diameter.

The radial velocities were then determined from the measurement of the midpoint of the relevant line and the Doppler equation. We then corrected these to the local standard of rest. The results for the stars that had detected spectra are contained in Tables D.1 and D.2. The errors for the radial velocities measured with SINFONI are large due to low signal-to-noise ratios and intrinsically large line widths (especially for faint sources). This means the resulting orbital elements have errors dominated by this measurement.

In total, we could determine radial velocities for 20 of the 54 stars, either by comparison to previously studied catalogs (e.g., Yelda et al. 2014; von Fellenberg et al. 2022) or by extracting the radial velocities from SINFONI spectra. Of the 20 stars, three stars already have orbital elements published in previous works, namely 3–74, 3–103, and 3–134, which are otherwise known as the reference stars R74, R30, and R34, respectively. These stars were left in the sample to serve as a control group. A total of ten of the stars for which radial velocities could be determined were found to be classified as old stars, while the remaining ten were determined to be young stars. Of the young stars, three had not been identified previously. For all stars that have measured radial velocities, we calculated the six orbital elements using all six phase-space coordinates up to the degeneracy induced by the unknown sign of the z coordinate.

It is often the convention to express orbital elements as the Keplerian parameters, which are: the semimajor axis, a, the eccentricity, e, the inclination, i, the argument of periapsis, ω, the longitude of the ascending node, Ω, and the epoch of the pericenter passage, t_p. The resulting orbital solutions for the 20 stars with radial velocity measurements are contained in Tables D.1 and D.2. We also visualized these orbits by plotting them together as shown in Fig. 3. With the addition of these stars at larger separations from Sgr A*, we have increased the number of known orbits within the GC by ∼40%.

Fig. 3. All stars with full orbital solutions determined using negative z coordinates (upper) and positive z coordinates (lower). The black plus sign denotes the location of Sgr A* at the (0, 0) coordinate, and NACO data points are marked by the filled circular points along each orbit. We also include the orbit of the S star S2 in red to provide a sense of scale.

If it were possible to observe a change in the radial velocity for these stars, we would be able to determine the correct sign for the z coordinate. Additionally, this would mean seven phase space coordinates are known while only six degrees of freedom apply to the calculation of the orbital elements for each star. With the additional phase-space coordinate, these stars could be used to constrain the potential of Sgr A* and, therefore, the SMBH mass. In this way, a prediction of when it would be possible to observe a significant change in radial velocity for each star is valuable. If it is assumed that a change of Δ v_rad = 50 km s⁻¹ would be robustly detectable, then it follows that the time that must elapse from now until such a change is observable is given by dividing the change in radial velocity by the radial acceleration, a_z, of the star. The choice of 50 km s⁻¹ is due to the fact that this change in radial velocity will be detectable using the Enhanced Resolution Imager and Spectrograph (ERIS), which has begun commissioning at the ESO VLT in 2022. The values obtained are provided for each star in Table 2.

5.2. Monte Carlo sampling over flat radial velocity prior

For stars that do not have measurable radial velocities, we used Monte Carlo sampling over a flat prior for the radial velocity to determine possible orbital solutions. The method follows closely that outlined by Bartko et al. (2009) and Lu et al. (2008). However, for the stars studied in Bartko et al. (2009) and Lu et al. (2008), the missing phase-space coordinate was the z coordinate, which has been determined in this work by detecting the curvature in the on-sky orbit trace.

In this case, the missing phase-space coordinate, here v_z, was subject to a series of assumptions. Firstly, we assumed that the star is gravitationally bound to Sgr A*. As such, a maximum velocity is readily calculated as the escape velocity of the star. If the maximum velocity vector is then ${\mathit{v}}_{\mathrm{esc}}={\mathit{v}}_{\max }=\sqrt{{\mathit{v}}_x^2+{\mathit{v}}_{\mathit{y}}^2+{\mathit{v}}_{z,\max }^2}$, the absolute value of v_z, max is then determined from

$$\begin{array}{c}\hfill |v_{z,\max }|=\sqrt{|\frac{2G{M}_{\bullet }}{r}-v_x^2-v_y^2|}.\end{array}$$(6)

In this case, we considered all solutions in the range v_z  ∈  [−v_esc, v_esc] by generating 20 000 velocities between these two values from a flat distribution. This was performed twice to accommodate the fact that it is again not possible to determine whether the z coordinate is positive or negative, which results in each star having 40 000 physically possible orbital solutions. We plotted these solutions onto a spherical Hammer-Aitoff projection map with axes corresponding to the longitude of the ascending node, Ω, and the inclination, i. The result is two curved traces representing possible orbital solutions for each star on the map.

An example of such a map of solutions is provided in Fig. 4. The coordinates used to produce Fig. 4 are those used in von Fellenberg et al. (2022). However, previous works (see, e.g., Paumard et al. 2006; Gillessen et al. 2009, 2017) use a slightly different definition, which is compared to the coordinates used in this work in Appendix E.

Fig. 4. Projected Ω vs. i traces determined via Monte Carlo sampling in radial velocity for the star 4–19. The ranges of the determined physically possible solutions are shown by the gray points, and the exact orbital solutions determined via radial velocity measurements from SINFONI data are indicated by the blue star markers for both the positive z (navy) and negative z (light blue) solutions. The errors in [x, y, z, vx, vy] are taken into account by sampling these values from a Gaussian distribution, but the resulting error bars are excluded from the figure. This consideration results in an uncertainty in Ω for the positive z solution of dΩ = 20.7° and dΩ = 15.8° for the negative z solution; for i it is given as di = 6.7° and di = 7.4° for positive and negative z solutions, respectively.

We could also determine the direction of rotation of the stellar orbit from these projections. More specifically, if a star has solutions located on the northern hemisphere of the map, then the star orbits clockwise, and if located in the southern hemisphere, the orbit is counterclockwise.

6. The Galactic Center environment

We can use the results discussed in the previous section to further understand the environment of the nuclear star cluster in the GC at these separations from Sgr A*. More specifically, in Sect. 6.1 we use the Hammer-Aitoff projection plots to compare previous works that looked for overdensities of stars in these (Ω, i) coordinates to indicate which stars could belong to disks or features in the GC environment (e.g., Bartko et al. 2009; Paumard et al. 2006; Yelda et al. 2014; von Fellenberg et al. 2022). Once disk membership is assumed for a star, we use this information to examine the mass enclosed within the orbits of these stars by removing the dependence of the black hole mass from the stellar z coordinate. We examine this in Sect. 6.2. Finally, by comparison to a simulated isotropic cluster, we investigate further structural information via trends in the orbits of the stars studied in Sect. 6.3.

6.1. Comparison to known disk structures

Figure 5 shows the Ω − i projection of spectroscopically young stars with known orbital elements corresponding to the solutions found in Sect. 5.1 that coincide with known overdensities in the GC environment. These observable disk features in the GC are those calculated in von Fellenberg et al. (2022; hereafter, vF22) and are visually summarized as ellipses in Fig. 5. The sample studied in vF22 contained stars with separations from Sgr A* up to 32 arcsec and, therefore, contains the regions studied in this work. As such, this map allows for a comparison between the stellar orbital solutions calculated in Sect. 5.1 and the known disk features in the GC. In this case, only stars that were identified as young stars using stellar spectra are considered since these features are typically populated only by young stars (see, e.g., Paumard et al. 2006; Bartko et al. 2010; Genzel et al. 2010; vF22). The features include the inner clockwise disk (CW1), the outer clockwise disk (CW2), the counterclockwise disk/filament (CCW/F1), as well as outer filament 2 (F2) and outer filament 3 (F3). In this case, only stellar orbital solutions that are consistent (within errors) with the values of i and Ω representing disk features are plotted. As such, Fig. 5 indicates that three of the stars in this work could belong to either CW2 or F2. We can use this information to predict the sign of the z coordinate by assuming the correct sign of z corresponds to the solution that results in the star belonging to a disk feature. This, of course, self-amplifies the disk features. However, it can be argued that this is a reasonable assumption since these features indicate overdensities in the stellar population in these coordinates. If a star has both positive and negative z solutions that are consistent with a feature, we cannot use this information to lift the orbital degeneracy. For this reason, only stars with a single z-solution consistent with a disk feature are listed in the summary of possible disk members provided in Table 3. The results are generally consistent with prior knowledge of the structure of the nuclear star cluster. It was determined by vF22 that particular features dominate at different separations from Sgr A*. From the stars studied here, it is clear that mostly stars within the first 2″ to 3″ are found on the outer clockwise disk. Unfortunately, Fig. 4 of vF22 also indicates that the slice with features at the lowest significance is from 4″ to 8″, meaning a smaller fraction of stars at these separations lie on features. Further, no single disk feature seems to dominate the density map at these separations. Since this largely represents the stars with larger separations studied in this work, it is not surprising that eight of the young stars with determined orbital parameters do not appear to be consistent with any features. Additionally, it could be argued that the source 3–103 is unlikely to belong to the Outer Filament 2, as this feature is most prominent at larger separations. Generally, the results of this study are consistent with the understanding of the nuclear star cluster presented in vF22.

Fig. 5. Projected Ω vs. i orbital solutions for the three young stars that have known orbital elements that correspond to a known GC disk feature (within error bars). All features are marked here by ellipses (see, e.g., Bartko et al. 2009; Yelda et al. 2014; von Fellenberg et al. 2022). Disk features are as calculated for Fig. 7 of von Fellenberg et al. (2022; hereafter, vF22).

We can use a similar analysis for stars that did not have measurable radial velocities to reveal which ranges of solutions appear to be consistent with the GC features. Figure 6 can be used to constrain the possible z values for stars that have solutions corresponding to disk features. Here, if a star without a measurable spectral radial velocity had orbital solutions determined from the Monte Carlo sampling of radial velocities that coincide with any of the known GC features, they were plotted onto Fig. 6. Again, we could use this information to predict the correct sign of the z coordinate for these stars. Table 4 provides a list of stars that only appear on disk features using one of the signs for the z coordinate. A total of eight stars have solutions that correspond to disk features if a negative z coordinate is used, while six with positive z coordinates result in solutions lying on a disk features. For stars that list multiple features, there may be one feature that is more suitable. For example, for stars that have solutions consistent with both the clockwise disk and the outer filaments, the solutions corresponding to the clockwise disk are more reasonable. This is because the CW disk is the dominant feature for stars at the separations from Sgr A* consistent with this sample, and the filaments are generally dominant for larger separations (vF22). For 20 stars, it was either the case that no solutions were consistent with disk features or that solutions corresponding to disk features were possible with both positive and negative z coordinates.

Fig. 6. Projected Ω vs. i traces determined via Monte Carlo sampling over a flat distribution in radial velocity for all stars that have a subset of solutions that corresponds to at least one of the known GC disk features using only a single sign for the z coordinate. The features are shown here by ellipses (see, e.g., Bartko et al. 2009; Yelda et al. 2014; vF22). Disk features are as calculated for Fig. 7 of vF22. The color of the trace indicates the 3D position vector of the source given by $r_{3\mathrm{D}}=\sqrt{x^2+{\mathit{y}}^2+z^2}$. This illustrates that the stars with the lowest r3D appear on the CW disk. This is consistent with the findings of vF22.

6.2. Using disk membership to estimate the mass distribution

The stars listed in Tables 3 and 4 along with the assumption of disk membership can be used to further understand the distribution of mass located within their orbits. We did this using the longitude of the ascending node and inclination of the disk features consistent with their original orbital solutions. Using these Ω and i values provides a novel method for calculating the z coordinate for each star that does not directly rely on the use of the accelerations calculated from the on-sky orbit trace. This means the final z coordinate will not depend on the mass of Sgr A*, as is the case for the absolute value of z calculated from Eq. (5). This was determined using

$$\begin{array}{c}\hfill z=\frac{x(sinisin\mathrm{\Omega })+y(sinicos\mathrm{\Omega })}{-cosi},\end{array}$$(7)

where x and y are the stellar positions in Dec. and RA, respectively, and the selected Ω and i coordinates are those that were determined in vF22 for the corresponding disk feature listed in Tables 3 and 4. We then calculated the enclosed mass within the orbit of the star by solving Eq. (5) for M_•. Figure 7 shows the median of the resulting masses enclosed within the orbits of each star as a function of 3D distance. Again, for stars with known spectroscopic ages, only young stars were included in the calculation. In this case, since we are only demonstrating this method as a possible way to estimate the mass distribution at these separations, a simplified approach is employed. This means, no marginalization over the priors was used to determine uncertainties for the mass enclosed in individual stellar orbits.

Fig. 7. Distribution of mass contained within the orbits of various sources. The median of the stars that have orbital solutions consistent with disk features in this work is shown by the red square marker. In this case, since the z values used are derived from the disk priors and we did not marginalize over priors, no formal error estimates for individual points can be provided. However, the uncertainty for the median of all the points in this work is estimated as the spread of the individual points. Also included are the values from the flare data (GRAVITY Collaboration 2018b), S2 from both Do et al. (2019) and GRAVITY Collaboration (2022), multiple S stars as measured with the VLT (GRAVITY Collaboration 2022), and stars that are members of the clockwise disk (GRAVITY Collaboration 2022). In this work, the enclosed mass is consistent within errors with the clockwise disk value, meaning no additional mass is required to account for the mass contained in the orbits of these stars. The blue line shows the sum of MSgr A* and the extended stellar mass distribution from previous works (e.g., Ghez et al. 2008; Gillessen et al. 2009; Genzel et al. 2010; Alexander 2017; Schödel et al. 2018; Baumgardt et al. 2018; GRAVITY Collaboration 2022).

The enclosed mass calculated for stars in this work is consistent within its error with the enclosed mass of the clockwise disk stars in GRAVITY Collaboration (2022). Since seven of the fourteen stars used in this calculation have solutions corresponding to either CW1 or CW2 and most of the stars are at separations consistent with the CW disk it is not surprising to see that the enclosed mass calculated here is consistent with this value. However, the points determined in this work do not further constrain the mass distribution since the required assumption of disk membership is heavily influenced by the uncertainties in radial velocity for stars with extracted spectral lines. For stars with unknown radial velocities, it is clear that the spread in possible orbital solutions also contributes to the error in the estimated mass enclosed within the potential orbits of these stars. This further highlights the importance of reducing the uncertainties in the stellar spectra used to calculate the radial velocities for these stars as well as needing better spectra for stars where radial velocities could not be determined. Future instrumentation with increased resolution promise to provide better estimates of disk membership and, thus, contribute to constraining the distribution of mass within the GC using stars with greater separations from Sgr A*.

6.3. Comparison to a simulated isotropic cluster

It is also useful to compare the studied sample of 54 stars with significant accelerations toward Sgr A* to a simulated isotropic cluster to determine the nature of the nuclear star cluster at these separations. The procedure to calculate the six phase-space coordinates for each cluster member is detailed in Genzel et al. (2003) and vF22, where first an isotropic distribution of orbits was determined, which we then converted to stellar positions and proper motions. The distributions representing the orbital elements for stars in the resulting isotropic cluster are shown in Fig. F.1.

In order to compare this simulated cluster to the data, we calculated the angular momentum statistic (or j statistic; Genzel et al. 2003) and the high-eccentricity statistic (or h statistic; Madigan et al. 2014) for stars in the sample as well as the simulated isotropic cluster members. Both statistics were developed as a means to determine estimates for the orbital elements of stars located far enough away from Sgr A* that full orbital solutions could not be determined with less than a decade of observation. In this work, the j and h statistics are used only to serve as a powerful visual comparison between the simulated isotropic cluster and the sample of stars found to have significant astrometric accelerations toward Sgr A*.

Figure 8 shows the comparison between the observed stars with significant accelerations and the isotropic cluster via histograms of calculated values of the j statistic. Here, the values of j = −1, 0 and 1 are shown by the vertical lines and correspond to CCW tangential, CW radial and CW tangential orbits, respectively. Interestingly, it appears that the observed stars exhibit a larger number of CCW tangential on-sky orbits than CW tangential orbits, whereas the isotropic cluster exhibits an equal fraction of both. We also see that a larger portion of the measured stars exhibit j values consistent with CW radial orbits as there is a surplus of stars that have j ∼ 0 when compared to the isotropic cluster. More quantitatively, when a Kolmogorov-Smirnov goodness-of-fit test (KS test) is performed using scipy.stats.ks_2samp (The SciPy Community 2008) to compare the calculated j values for both the observed stars and the simulated isotropic cluster, a KS score of ∼0.12 is returned. This corresponds to a p value of ∼9.4 × 10⁻⁹, meaning we can reject the null hypothesis in favor of the scenario where these values are drawn from different distributions.

Fig. 8. Number of stars binned as a function of the calculated j statistic for both the isotropic cluster (blue) and observed stars (magenta). Vertical lines indicate CW tangential (yellow), CCW tangential (orange), and CW radial (red) on-sky projected orbits at 1, −1, and 0, respectively. These results suggest an overdensity of stars with CCW tangential orbits compared to the isotropic cluster that displays equal amounts of CW and CCW tangential orbits. To best model observations, all members of the isotropic cluster with separations of < 2 arcsec from Sgr A* are discarded before plotting.

This inconsistency with the simulated isotropic cluster is also demonstrated in Fig. 9. The vertical lines in Fig. 9 indicate low |h| values (orange) corresponding to radial orbits, high |h| values (yellow) indicating tangential orbits and unbound orbits (red) at $|h|=\sqrt{2}$. None of the stars appear above the $|h|=\sqrt{2}$ line due to the condition applied when searching for significant accelerations, which excludes all nonphysical (i.e., unbound) accelerations and the fact that the cluster members are bound to Sgr A* by definition. Here it is clear that the observed stars have a different distribution of |h| values than the isotropic cluster, shifted toward lower |h| values. Comparison to the isotropic cluster also indicates that this sample may be biased to high eccentricity orbits since a higher proportion of measured stars feature |h|∼0. Again, using the KS test now to compare the |h| values for the simulated cluster and the observed stars, a KS score of ∼0.32 is returned corresponding to a p value of ∼2.0 × 10⁻¹⁵. This again indicates that we should reject the null hypothesis in this case and conclude the observed stars do not originate from the same isotropic distribution as the simulated cluster. Overall, the results suggest that the observed sample of stars is not well represented by an isotropic cluster.

Fig. 9. Number of stars binned as a function of the calculated |h| statistic for both the isotropic cluster (blue) and observed stars (magenta). The vertical lines in this case indicate high |h| values corresponding to tangential orbits (yellow), low |h| values corresponding to radial orbits (orange), and unbound orbits (red). This indicates that radial orbits (or low |h| values) are more common within the sample than orbits within the simulated isotropic cluster. Again, to best model observations, all members of the isotropic cluster with separations of < 2 arcsec from Sgr A* are discarded before plotting.

In summary, these comparisons to either known disk features in the GC environment using Hammer-Aitoff projection maps or the simulated isotropic cluster using the j and h statistics are consistent with results of previous studies that have found that the GC environment is round but not isotropic due to disks and features being present in the stellar population. Further, disk membership could be a useful tool to lift orbital degeneracies and constrain the mass enclosed within the orbits of stars with full orbital solutions. However, further work is required to increase the precision of the radial velocity measurements in order to produce meaningful constraints on the enclosed mass in the GC.

7. Summary

This work extends the long-term monitoring program of stars in the GC to those located at greater separations from Sgr A*. In order to accomplish this, we created a source finder that can track more than 1800 stars over 17 years of NACO data within the complex environment represented by the GC nuclear star cluster. We then calculated the proper motions for these stars using position information across these epochs in order to determine second-order fits for the motion of these stars in rotated radial and tangential components of motion. A subset that includes 54 of these stars presented significant astrometric accelerations in the direction of Sgr A*. For these stars, we used the measured 2D accelerations to determine the absolute value of the z coordinate for each source.

Further, we determined radial velocities from SINFONI data for 20 of the stars with significant astrometric accelerations. The radial velocities could then be used to transform the six phase-space coordinates into orbital elements that describe the motion of the stars around Sgr A*. This included three previously studied reference stars, as well as 17 stars that previously had no determined orbits. In total, we determined that ten of the stars are spectroscopically old and the remaining ten young. Three of these young stars had no previous spectroscopic classification.

We also considered the orbital solutions of stars without radial velocity measurements using Monte Carlo sampling. By plotting the sets of physically possible solutions on Hammer-Aitoff projection maps, we identified trends in the orbits. Almost all of the stars were found to include ranges of orbital solutions that could classify them as disk members. If only one sign for the z coordinate of these stars resulted in solutions that were consistent with a disk feature, the corresponding sign for the z position was selected and the resulting orbit was used to further understand the GC environment. This was the case for 20 5D constrained stars. Of the stars with orbital solutions, three had values of inclination and longitude of the ascending node that corresponded to one of three disk features, namely CW1, CW2, or F2 (although the last case is unlikely due to separation considerations). We then assumed the disk membership to be true in order to derive a z coordinate that depends only on the i and Ω of the disk to which the star belongs and the stellar positions in RA and Dec. We then derived the mass enclosed within the orbit of the star. We find that the median enclosed mass we calculated, when compared to previous works, agrees within uncertainties with the expected mass enclosed in the clockwise disk (GRAVITY Collaboration 2022). However, this measurement is not able to constrain the enclosed mass. This is likely due to the fact that the radial velocities for these stars are too uncertain. The GC environment was also probed using the j and h statistics of the stars with significant radial accelerations toward Sgr A* by comparing them to those calculated for members of a simulated isotropic cluster. We find that the sample is not well represented by the modeled isotropic cluster, which is consistent with the fact that disks and features are observed in the GC (Paumard et al. 2006; Bartko et al. 2010; Genzel et al. 2010; vF22). Future observations could be used to further this work, particularly in the field of IR spectroscopy. Within the next ten years, for example, it should be possible to measure a change in radial velocity greater than 50 km s⁻¹ for the stars contained in Table 2, namely: 3–19, 4–120, 3–20, 4–56, 4–36, 3–53, 3–91, 5–246, 3–130, 3–161, 3–158, and 3–134. This would make it possible to constrain the z coordinate for these stars, removing the degeneracy of their orbital solutions as well as allowing them to contribute to constraining the potential of Sgr A*.

One instrument that will be crucial to studies of the GC in the coming years is ERIS, the installation of which was completed at the UT4 of the ESO VLT in early 2022 and which will serve as the replacement for NACO and SINFONI (Davies et al. 2018). ERIS contains an upgraded version of SPIFFI, namely SPIFFIER, which should result in better throughput and higher spectral resolution. The NIR camera, NIX, is also expected to have a comparable performance to NACO. Specifically, ERIS is optimized for GC observing with the resolution and wavelength required for the continued monitoring of fainter and closer stars (Davies et al. 2018).

Another exciting future instrument for the monitoring of stars featured in this work is the Multi-AO Imaging Camera for Deep Observations (MICADO), which will be the first instrument installed on the ESO Extremely Large Telescope (ELT; Davies et al. 2021). MICADO will have a resolution of only three times less than the current resolution of GRAVITY while also having a FOV comparable to that of NACO (Davies et al. 2021). As such, MICADO will peer deeper into the outer regions of the GC than any previous instrument and will be essential to furthering this work. Beyond imaging, it is necessary to acquire spectra for stars that feature significant accelerations toward Sgr A*. Definite orbital solutions could still not be produced for 34 of the 5D constrained stars detected in this work, due to a lack of radial velocity measurements from stellar spectra. As such, spectroscopy from ERIS and, further into the future, MICADO and HARMONI on the ELT promises to provide improved resolution compared to SINFONI, which will accurately measure radial velocities for these stars. This increased precision will also contribute to improving the constraints on the enclosed mass within the orbits of stars at these separations.

In conclusion, this work not only increases the number of known orbits within the GC by ∼40%, it also confirms previous studies that have found evidence for multiple disk structures within the nuclear star cluster environment. In addition, this work provides strong scientific justification for follow-up observations using future instrumentation that will continue to provide valuable information about the stars at greater separations from Sgr A*.

Acknowledgments

This publication is based on observations collected at the European Southern Observatory and the authors thank ESO and the ESO/Paranal staff, without whom this work would not be possible. We also thank the referee for providing useful comments which helped to improve the paper. A.Y. acknowledges support from the European Space Agency’s Young Graduate Trainee program. J.S. was supported in part by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy - EXC-2094 - 390783311.

References

Appendix A: List of NACO data sets

Appendix B: List of SINFONI data sets

Appendix C: Polynomial fits for RA and Dec. motion of stars with significant astrometric acceleration

Appendix D: Orbital elements for stars with measured radial velocities

Appendix E: Coordinate system used for Hammer-Aitoff diagrams

The coordinate system used to plot the Hammer-Aitoff projection maps in this work follows the conventions used in von Fellenberg et al. (2022). It relates to previous works that use different conventions for the plotted coordinates such as Paumard et al. (2006), Gillessen et al. (2009), and Gillessen et al. (2017) by rotating Ω by 90 degrees relative to the coordinate system used in this work. To demonstrate this, Fig. E.1 shows the equivalent of Fig. 5 in the coordinate system used in previous works.

Fig. E.1. Demonstration of the coordinate system used in previous works for Hammer-Aitoff projection maps (see, e.g., Gillessen et al. 2009, 2017) by mapping Fig. 5 to this coordinate system. The transformation is most simply described as a 90 degree rotation in Ω. The projected (Ω, i) coordinate corresponding to the orbital solution for S2 from Gillessen et al. (2017) is also included for reference (green) in addition to the points representing the three stars also shown in Fig. 5 (navy).

Appendix F: Distributions used to define the simulated isotropic cluster

Fig. F.1. Distributions of orbital elements used in Sect. 6.3 to the generate isotropic cluster, including the semimajor axis, eccentricity, inclination, argument of periapsis, longitude of the ascending node, and true anomaly. By definition, all stars in the cluster are bound to the potential of Sgr A*.

All Tables

All Figures

Fig. 1. Illustration of the source finding algorithm. First, the position from Epoch(N) is projected to Epoch(N+1) using the polynomial transformations. The transformed position is then used as an initial guess for the finder. Each grid square in the zoomed-in panel, which shows a NACO source profile, corresponds to a Gaussian search box of some side length. The fitter tries fitting a Gaussian profile to the pixels in each of the surrounding grid squares. The Gaussian fit with the lowest χ2 (white circle) is then selected as the position for this star in the new epoch. For clarity, the search boxes appear here as a grid; however, the actual algorithm utilizes overlapping search boxes.

In the text

Fig. 2. Histogram illustrating the distribution of tangential acceleration significance within the sample of stars examined. The histogram is normalized such that the area under the curve is equal to 1. The tangential significance for each star is calculated from Eq. (2). Two Gaussian fits to the histogram are plotted. The first includes the wings of the histogram when finding a best fit (solid curve) that is dominated by sources with nonphysical accelerations, while the second excludes the data located in the wings of the histogram (dashed curve). Therefore, the standard deviation for the curve, which excludes the wings, was calculated and used to normalize the radial acceleration significance. In this case, the standard deviation is found to be ∼1.212.

In the text

	Fig. 3. All stars with full orbital solutions determined using negative z coordinates (upper) and positive z coordinates (lower). The black plus sign denotes the location of Sgr A* at the (0, 0) coordinate, and NACO data points are marked by the filled circular points along each orbit. We also include the orbit of the S star S2 in red to provide a sense of scale.
In the text

Fig. 4. Projected Ω vs. i traces determined via Monte Carlo sampling in radial velocity for the star 4–19. The ranges of the determined physically possible solutions are shown by the gray points, and the exact orbital solutions determined via radial velocity measurements from SINFONI data are indicated by the blue star markers for both the positive z (navy) and negative z (light blue) solutions. The errors in [x, y, z, vx, vy] are taken into account by sampling these values from a Gaussian distribution, but the resulting error bars are excluded from the figure. This consideration results in an uncertainty in Ω for the positive z solution of dΩ = 20.7° and dΩ = 15.8° for the negative z solution; for i it is given as di = 6.7° and di = 7.4° for positive and negative z solutions, respectively.

In the text

	Fig. 5. Projected Ω vs. i orbital solutions for the three young stars that have known orbital elements that correspond to a known GC disk feature (within error bars). All features are marked here by ellipses (see, e.g., Bartko et al. 2009; Yelda et al. 2014; von Fellenberg et al. 2022). Disk features are as calculated for Fig. 7 of von Fellenberg et al. (2022; hereafter, vF22).
In the text

Fig. 6. Projected Ω vs. i traces determined via Monte Carlo sampling over a flat distribution in radial velocity for all stars that have a subset of solutions that corresponds to at least one of the known GC disk features using only a single sign for the z coordinate. The features are shown here by ellipses (see, e.g., Bartko et al. 2009; Yelda et al. 2014; vF22). Disk features are as calculated for Fig. 7 of vF22. The color of the trace indicates the 3D position vector of the source given by $r_{3\mathrm{D}}=\sqrt{x^2+{\mathit{y}}^2+z^2}$. This illustrates that the stars with the lowest r3D appear on the CW disk. This is consistent with the findings of vF22.

In the text

Fig. 7. Distribution of mass contained within the orbits of various sources. The median of the stars that have orbital solutions consistent with disk features in this work is shown by the red square marker. In this case, since the z values used are derived from the disk priors and we did not marginalize over priors, no formal error estimates for individual points can be provided. However, the uncertainty for the median of all the points in this work is estimated as the spread of the individual points. Also included are the values from the flare data (GRAVITY Collaboration 2018b), S2 from both Do et al. (2019) and GRAVITY Collaboration (2022), multiple S stars as measured with the VLT (GRAVITY Collaboration 2022), and stars that are members of the clockwise disk (GRAVITY Collaboration 2022). In this work, the enclosed mass is consistent within errors with the clockwise disk value, meaning no additional mass is required to account for the mass contained in the orbits of these stars. The blue line shows the sum of MSgr A* and the extended stellar mass distribution from previous works (e.g., Ghez et al. 2008; Gillessen et al. 2009; Genzel et al. 2010; Alexander 2017; Schödel et al. 2018; Baumgardt et al. 2018; GRAVITY Collaboration 2022).

In the text

Fig. 8. Number of stars binned as a function of the calculated j statistic for both the isotropic cluster (blue) and observed stars (magenta). Vertical lines indicate CW tangential (yellow), CCW tangential (orange), and CW radial (red) on-sky projected orbits at 1, −1, and 0, respectively. These results suggest an overdensity of stars with CCW tangential orbits compared to the isotropic cluster that displays equal amounts of CW and CCW tangential orbits. To best model observations, all members of the isotropic cluster with separations of < 2 arcsec from Sgr A* are discarded before plotting.

In the text

Fig. 9. Number of stars binned as a function of the calculated |h| statistic for both the isotropic cluster (blue) and observed stars (magenta). The vertical lines in this case indicate high |h| values corresponding to tangential orbits (yellow), low |h| values corresponding to radial orbits (orange), and unbound orbits (red). This indicates that radial orbits (or low |h| values) are more common within the sample than orbits within the simulated isotropic cluster. Again, to best model observations, all members of the isotropic cluster with separations of < 2 arcsec from Sgr A* are discarded before plotting.

In the text

Fig. E.1. Demonstration of the coordinate system used in previous works for Hammer-Aitoff projection maps (see, e.g., Gillessen et al. 2009, 2017) by mapping Fig. 5 to this coordinate system. The transformation is most simply described as a 90 degree rotation in Ω. The projected (Ω, i) coordinate corresponding to the orbital solution for S2 from Gillessen et al. (2017) is also included for reference (green) in addition to the points representing the three stars also shown in Fig. 5 (navy).

In the text

	Fig. F.1. Distributions of orbital elements used in Sect. 6.3 to the generate isotropic cluster, including the semimajor axis, eccentricity, inclination, argument of periapsis, longitude of the ascending node, and true anomaly. By definition, all stars in the cluster are bound to the potential of Sgr A*.
In the text