The Orbital Structure and Selection Effects of the Galactic Center S-star Cluster

The galactic center hosts one of the most mysterious structures of the Milky Way, the S-star cluster (SSC; see Genzel et al. 2010 for a review). This spheroidal ensemble of at least 50 young, preferentially B-type stars with masses in the range of 8–20 M_⊙ and ages of less than 15 Myr (Gillessen et al. 2017; Habibi et al. 2017) is confined to a region of 0.04 pc from Sgr A* (Genzel et al. 2003; Ghez et al. 2003a; Gillessen et al. 2009). Their eccentric to nearly parabolic orbits have been reconstructed based on proper motion measurements (Genzel et al. 1997; Ghez et al. 1998, 2003b, 2008; Gillessen et al. 2009; Schödel et al. 2009; Genzel et al. 2010; Boehle et al. 2016). This kinematic data not only confirms the existence of a super-massive black hole (SMBH) with a mass M_• = 4.297 × 10⁶ M_⊙ at the Galactic Center (Genzel et al. 1997; Ghez et al. 1998, 2005; Schödel et al. 2002; Eisenhauer et al. 2005). It also provides valuable clues on the dynamical evolution of these recently formed stars in its proximity (Levin & Beloborodov 2003).

The origin of the SSC is puzzling as single stars typically condense out of molecular cloud cores with radii of order 0.1 pc (Hester et al. 1996; Mellema et al. 2006; Forbes et al. 2021). This cloud size is larger than the whole SSC. Therefore, the S-stars either formed under very exotic conditions (Goodman 2003; Thompson et al. 2005; Levin 2007; Nayakshin et al. 2007; Cantiello et al. 2021; Ali-Dib & Lin 2023) or were transported into this tiny, hostile environment (Syer et al. 1991; Artymowicz et al. 1993; Davies & Lin 2020) of the central SMBH. If the S-stars formed in situ, some violent star formation event must have happened very close to the SMBH a few megayears ago (Morris 1993; Ghez et al. 2003a). Currently, the S-stars are embedded in a 10⁶ K hot, diffuse, gas bubble (Xu et al. 2006; Murchikova et al. 2019) which is not a favorite condition for star formation. One should, however, note that not all of the gas is hot. A cluster of cold G-clouds with temperatures of order 10⁴ K has been discovered in the same region as the SSC (Gillessen et al. 2012, 2013; Pfuhl et al. 2015; Plewa et al. 2017). These gas clumps have very small masses, of order several Earth masses. Their nature and origin are still a matter of debate (Burkert et al. 2012; Murray-Clay & Loeb 2012; Ballone et al. 2013; Scoville & Burkert 2013; Owen & Lin 2023) and at the moment it is not clear whether they are related to the origin of the S-stars.

Today, nearly 200 massive stars have been detected in the Sgr A* proximity. In addition to the SSC, ∼20% of the known stars belong to a subgroup of O/WR stars with top-heavy luminosity function and estimated age ∼6 ± 2 Myr (Ghez et al. 2003a). They reside in a possibly warped disk with clockwise rotation, commonly referred to as the clockwise disk (CWSs, Levin & Beloborodov 2003; Lu et al. 2006; Paumard et al. 2006; Bartko et al. 2009; Löckmann & Baumgardt 2009; Lu et al. 2009; Bartko et al. 2010; Yelda et al. 2014). The disk stars have semimajor axis in the range of 0.04 pc ≲ a_⋆ ≲ 0.5 pc with eccentricities of order e_⋆ ∼ 0.2–0.4. Interestingly, the smallest values of a_⋆ = 0.04 pc correspond to the sharp outer edge in the semimajor axis distribution of the SSC. This coexistence and the very similar ages provide evidence that the SSC and its surrounding disk population might have a common origin. The remaining fraction of the young stars are the "off-disk" stars (ODSs) with much less well-determined, probably large eccentricities and inclinations, as well as a non-isotropic distribution, surrounding the disk-star population. Several over-dense concentrations in their orbital angular-momentum orientation distribution have been identified, including the controversial identification of a second disk, inclined by ∼110° with respect to the clockwise disk (Paumard et al. 2006; Lu et al. 2009; Yelda et al. 2014; Ali et al. 2020). The luminosity function of the "off-disk" stars is similar to that of the S-stars and less top-heavy than the disk stars (von Fellenberg et al. 2022). There are suggestions that two disks of massive stars formed 6 Myr ago, orbiting the central SMBH from one or more infalling and colliding clouds (Bonnell & Rice 2008; Hobbs & Nayakshin 2009; Alig et al. 2013). It is also possible that secular perturbation by an intermediate-mass companion (Zheng et al. 2020, 2021; GRAVITY Collaboration et al. 2023) and dynamical relaxation (Rauch & Tremaine 1996; Freitag et al. 2006; Alexander et al. 2007; Kocsis & Tremaine 2011, 2015; Madigan et al. 2011) heated a common-natal disk, generating all the diverse kinematic properties in each subgroup, including the spheroidal SSC. We will discuss this scenario in greater detail in a companion paper (X. Zheng et al. 2024, in preparation).

Another, maybe complementary, scenario for the origin of the S-stars is the Hills mechanism, whereby massive stellar binaries, generated in the larger environment of the SMBH were deflected by some process onto an orbit that brought them so close to the SMBH that they broke up (Hills 1975, 1988; Gould & Quillen 2003; Yu & Tremaine 2003; Ginsburg & Loeb 2006; Amaro-Seoane et al. 2012). Typically the more massive star will lose orbital energy and go into a tight, bound orbit around the SMBH while the other star gains energy and might even be ejected as a high-velocity star or a hyper-velocity star that is unbound to the Milky Way. Stars with such high velocities have been detected in the galactic halo. Observations indicate that these stars have their origin in the galactic center, consistent with this theory. One caveat of the Hills mechanism is, however, the fact that it alone can only explain the current orbits of the most bound S-stars. One would still need secular evolution scatter processes to account for the S-stars at larger distances. In addition, this mechanism does not provide an answer to the question of why almost all S-stars are B-stars, whereas the Hills mechanism should also work equally well for O-stars.

The size and geometry of the SSC provide insight into its origin. More detailed information is, however, available from measurements of the orbital distribution of its stars. By now, we have up to 25 yr of uninterrupted monitoring of stellar orbits, and the number of stars with known orbits (Schödel et al. 2002; Ghez et al. 2003a) has increased to ∼50 (Eisenhauer et al. 2005; Ghez et al. 2005; Gillessen et al. 2009, 2017). One of the most prominent members of the SSC is S2 which passed its first observed pericenter in 2002 and again in 2018, providing the most precise measurement of the black hole mass and allowing for the first time to detect gravitational redshift (GRAVITY Collaboration et al. 2018; Do et al. 2019) and first-order Schwarzschild precession (GRAVITY Collaboration et al. 2020), in excellent agreement with the predictions by general relativity. Unfortunately, as pointed out by Capuzzo-Dolcetta & Sadun-Bordoni (2023), S2 travels too far from the SMBH in order to pick up the contribution of higher-order terms to Schwarzschild precession. S2 and other S-stars with small pericenter distances are key objects in order to test the predictions of general relativity. In order to understand the formation of the SSC, one however needs to investigate the physical property of the whole cluster, which requires a detailed understanding of biasing and selection effects.

In this paper, we therefore focus on the distribution of the orbital parameters of the SSC. In Section 2, we demonstrate that the S-stars follow a tight correlation between their orbital eccentricity e_⋆ and their pericenter distance r_p. Interestingly the G-clouds follow the same distribution, which might indicate a common origin. In Section 3, we show that this correlation differs strongly from a thermal distribution, which is usually adopted to characterize the SSC or other theoretical models that explain the SSC as the result of secular heating of a previous disk of stars. We also demonstrate that the lack of stars with high e_⋆ and large r_p can be understood as a result of selection effects. Stars could indeed exist in this regime. One would, however, need observations on a longer timeline to determine their orbital parameters. The zone of avoidance, on the other hand, is real. It is defined by a sharp edge of r_p versus e_⋆, below which no S-stars exist. A summary follows in Section 4.

The red points with error bars in Figure 1 show the logarithm of r_p versus e_⋆ for the S-stars, compiled by Gillessen et al. (2017). The pericenter increases with increasing (1 − e_⋆). In addition to the stars, cyan triangles mark the currently known G-gas clouds (Gillessen et al. 2012; Eckart et al. 2013; Plewa et al. 2017; Ciurlo et al. 2020; Peißker et al. 2020) which follow the same correlation. This might indicate a common origin of S-stars and G-clouds or a common secular redistribution process that generated the currently observed configuration. One should however note that the number of G-clouds is still small. In addition, scattering processes, if important, should be more disruptive for the G-clouds if they are diffuse, pressure-confined gas clumps, with negligible self-gravity (Murray et al. 1993; Burkert et al. 2012).

Zoom In Zoom Out Reset image size

One might expect a dependence of r_p on e_⋆ if all the stars have a similar semimajor axis a_⋆ and a spread in e_⋆, as for given a_⋆ the pericenter distance r_p = a_⋆ × (1 − e_⋆) decreases with increasing e_⋆. In this case, r_p should however increase linearly with (1 − e_⋆), which is different from the trend, seen in Figure 1. To demonstrate this more clearly, the two dashed black lines in Figure 1 depict the expected correlation for a semimajor axis of a_⋆ = 2 × 10³ AU and a_⋆ = 2 × 10⁴ AU, respectively. One can see that the more circular the orbits, the larger the semimajor axis. This is demonstrated more clearly in Figure 2, which shows the eccentricity dependence of a_⋆ (left panel) and of the apocenter distance r_a = a_⋆ × (1 + e_⋆). Both quantities are not randomly distributed but depend on (1 − e_⋆). This indicates that the dependence of pericenter distance on eccentricity is not just a result of a universal semimajor axis, independent of eccentricity. One should however note that the number of data points is still small and the error bars are large.

Zoom In Zoom Out Reset image size

One possible explanation of the correlations seen in Figures 1 and 2 is selection effects. It might, for example, be harder to determine the complete orbital parameters of stars in the upper left and lower right corners of Figure 1. In order to test this conjecture, we have performed a detailed Monte Carlo study. We start by focusing on a pair of r_p and (1 − e_⋆). Stars are then generated with this set of orbital parameters on randomly oriented orbital planes. Each star is placed somewhere on its orbit with the starting point at time t₀ chosen randomly between zero and the orbital period. The star is then advanced by 20 yr which corresponds to the average observational period and its projected position (x_i , y_i ) at time t_i = t₀ + (i − 1)Δt in the x − y plane is recorded with Δt = 1 yr. Here we assume that the z-axis represents the line of sight to the observer. Next, a Gaussian random error of 10 AU is added to each position which is the expected mean observational error. From this, we determine the time dependence of the distance to the central black hole ${r}_{i}({t}_{i})=\sqrt{{x}_{i}^{2}+{y}_{i}^{2}}$. A second-order regression line r(t) = a + b × t + c × t² is now fitted through the 21 data points r_i (t_i ). Repeating this procedure a large number of times allows us to determine the mean value 〈c〉 and its error δ c. It is this curvature 〈c〉 in the orbit that is required in order to determine the full orbital parameters of a star. We now follow the selection criterion of Gillessen et al. (2017) and assume that the orbital parameters are known if δ c/〈c〉 ≤ 10⁻³. The fraction of cases that fulfill this requirement then defines the probability P to determine the orbital parameters for a given of r_p and (1 − e_⋆). For P = 1, orbital parameters can be determined for all configurations. The smaller P, the less frequently can these parameters be specified. Regions with small values of P should not be populated in Figure 1, even if stars exist on these orbits, because the curvature c and by this, the orbital parameters r_p and e_⋆ cannot be determined reliably for a large fraction of orbital configurations.

We find that, to good approximation, the probability P to measure the orbital parameters of a star with pericenter distance r_p in units of AU and orbital eccentricity e_⋆ is given by the following fit formula

$$\begin{eqnarray}&&\mathrm{log}P=-2+\left(\displaystyle \frac{-{q}_{2}+\sqrt{{q}_{2}^{2}-4{q}_{1}{q}_{3}}}{2{q}_{1}}\right)\end{eqnarray} \tag{ 1 }$$

with

$$\begin{eqnarray}\begin{array}{rcl}{q}_{1} & = & 0.095\times x\\ {q}_{2} & = & 1.266-0.95\times x-0.075\times x\times \mathrm{ln}x\\ {q}_{3} & = & \mathrm{ln}({r}_{p}/\mathrm{AU})-11.4+0.855\times x\\ & & +0.075\times x\times \mathrm{ln}x-\mathrm{ln}x\end{array}\end{eqnarray} \tag{ 2 }$$

and x = (1 − e_⋆).

The dotted lines in Figure 3 show contours of constant P. Above a given contour, the probability is smaller than the corresponding value. Red points show the S-stars with known orbital parameters. One can see that the unpopulated upper left region with high eccentricities and high r_p is characterized by very small probabilities of a few percent or lower. Indeed, almost no stars lie in this region. And those few stars that are found there have large error bars. We, therefore, cannot rule out that this region is populated with S-stars. However, determining their orbital parameters will require a longer epoch of observations and/or higher observational precision. On the other hand, the orbital parameters of S-stars in the empty lower right triangular region below the dashed line should by now have been measured. We call this region the zone of avoidance. It is characterized by

$$\begin{eqnarray}&&\mathrm{log}({r}_{{\rm{p}}}/\mathrm{AU})\leqslant 1.57+2.6(1-{e}_{\star })\end{eqnarray} \tag{ 3 }$$

with the dashed line showing its upper boundary. That this zone of avoidance is not populated by stars therefore indicates that no S-stars exist with these combinations of orbital parameters. The zone of avoidance overlaps with the P = 1 line (solid black line) close to e_⋆ = 0. The orbits of stars on circular orbits above the zone are, therefore, harder to determine as the probability to measure their parameters within on average 20 yr is small. This might account for the fact that not many stars with eccentricities close to unity have been reliably detected yet.

Zoom In Zoom Out Reset image size

To more quantitatively evaluate the effect of the detection probability on the distribution of observed orbital parameters, Figure 4 shows a Monte Carlo test where we randomly draw values of e_⋆ in the regime of 0 and 1 and values of $\mathrm{log}({r}_{{\rm{p}}}/\mathrm{AU})$ within a range of 1 and 5. We then discard this set of parameters if it lies in the zone of avoidance (Equation (3)). Otherwise, we accept the set with a probability as given by Equation (2). Note that this procedure makes the result insensitive to the exact values of the upper and lower boundary of the adopted $\mathrm{log}({r}_{{\rm{p}}})$ as long as the lower boundary lies within the zone of avoidance and the upper boundary lies in the region where the probability of determining the orbital parameters is small. As soon as 41 stars have been selected this way, we determine their cumulative eccentricity distribution. We repeat this procedure 1000 times. The cyan points in Figure 4 show the mean of the cumulative (1 − e_⋆) distribution with the error bars depicting the 1σ mean square variation. The distribution within the error bars is very similar to the observed distribution shown by the red points. This also indicates that a lot of stars could hide in the upper left corner at high pericenter distances and large eccentricities. We have repeated this test, adopting a Rayleigh distribution which is characterized by the probability

$$\begin{eqnarray}&&P({e}_{\star })=\left(\displaystyle \frac{{e}_{\star }}{\sigma }\right)\times \exp \left[0.5\times \left(1,-,{\left(\displaystyle \frac{{e}_{\star }}{\sigma }\right)}^{2}\right)\right]\end{eqnarray} \tag{ 4 }$$

The best fit in this case corresponds to σ = 0.6 (yellow points). Despite this additional free parameter it cannot fit the data as well as the random distribution. Finally, the open points show an isothermal distribution P(e_⋆) ∼ e_⋆, which cannot fit the data at all.

Zoom In Zoom Out Reset image size

We demonstrate that the Galactic center S-star cluster shows a narrow dependence of pericenter distance r_p on (1 − e_⋆), with e_⋆ the orbital eccentricity. Based on their age and luminosity-function similarities, we adopt the assumption that all the young stars, including the S-stars, CWSs, and ODSs, formed in a common disk with nearly circular orbits (see Section 1). In this scenario, their present-day kinematics would imply that the S-stars have lost a substantial fraction of their total initial angular momentum or have markedly increased the magnitude of their total (negative) energy. Although such a large amount of angular-momentum deficit can be acquired through classical relaxation due to uncorrelated two-body distant scattering and close encounters (Binney & Tremaine 1987), it would occur on timescales longer than the age of the S-stars (Freitag et al. 2006).

Mutual gravitational interaction between stars in a nuclear cluster around a SMBH leads to the precession of their pericenter longitudes on timescales much longer than their orbital periods. This precession is correlated with angular-momentum transfer and eccentricity changes. Between successive pericenter passages, the net amount of angular-momentum exchange accumulates, with an insignificant amount of energy modifications. Moreover it modulates on a secular resonant-relaxation timescale over many orbits (Rauch & Tremaine 1996). With their orbital periods in the range from decades to centuries, the S-stars undergo ≳10⁴ orbits around Sgr A* such that resonant relaxation has sufficient time to proceed (Kocsis & Tremaine 2011). However, its effect on stars with nearly circular orbits is negligible since their angular momenta (per energy) are already saturated and cannot be significantly changed. Thus, resonant relaxation alone cannot lead to the excitation of eccentricity among an isolated population of stars with a negligible angular-momentum deficit (i.e., all on nearly circular orbits).

In our previous studies (Zheng et al. 2020, 2021), we showed that some stars may have lost a substantial fraction of their initial angular momentum to Sgr A^⋆'s intermediate-mass companion (IMC), IRS 13E (on an eccentric orbit which is inclined to these stars' natal disk) under the influence of the IMC's von-Zeiple-Lidov-Kozai and sweeping secular resonance during the depletion of their natal disk (Zheng et al. 2020, 2021). In an accompanying paper (X. Zheng et al. 2024, in preparation), we verify that some stars, formed in the proximity of the IMC's orbit, lose a large fraction of their initial angular momentum and encroach into the S-star domain during their pericenter passages. These highly eccentric migrant stars inject a substantial angular-momentum deficit into the inclusive S-star population. Subsequently, resonant relaxation among the intruding and indigenous S-stars randomizes their orbital angular momentum, isotropizes the orientation of their orbit-normal vectors, and leads to their present-day kinematic properties within their current age.

The G-clouds follow a similar correlation, indicating that they might be of similar origin or have experienced similar dynamical evolution. More detections of G-clouds are required in order to test this conjecture. If scattering and secular relaxation mechanisms were important, the fact that the G-clouds were not disrupted would indicate that they might have an invisible, compact central source which could replenish gas that is tidally stripped and which could also provide a deep enough potential well to gravitationally bind the clouds' gas (Scoville & Burkert 2013; Peißker et al. 2021). Whether the G-clouds would be disrupted without a central gravitational anchor is a question that we plan to explore in a subsequent paper.

We then derive the likelihood to determine orbital parameters for S-stars observationally and show that small values of (1 − e_⋆), combined with large pericenter distances r_p require a larger observational time span than what is available at the moment. This explains why no stars have yet been found in this region. It might be populated with yet undetected S-stars. On the other hand, the orbital parameters of stars characterized by $\mathrm{log}({r}_{{\rm{p}}}/\mathrm{AU})\lt 1.57+2.6(1-{e}_{\star })$ should be easy to determine. That no stars exist in this zone of avoidance indicates that whatever dynamical processes shaped the S-star cluster kinematics did not favor this region. Considering the observational biases and the zone of avoidance, we find that the S-stars follow a random sampling of $\mathrm{log}({r}_{{\rm{p}}})$ with (1 − e_⋆). Our results could provide valuable constraints for theoretical models of the origin of the S-star cluster and its kinematics. As discussed above we will show in an accompanying paper that this correlation arises naturally from the perturbation of an intermediate-mass companion (X. Zheng et al. 2024, in preparation) without invoking any tidal-disruption-of-binary hypothesis (Hills 1988; Chu et al. 2023).

This work was supported by the Excellence Cluster ORIGINS which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy—EXC-2094-390783311.