The S stars' zone of avoidance in the Galactic center

Introduction

The Galactic center’s proximity allows for the direct measurement of stellar orbits within the central parsec. There are distinct structures within this region: multiple young stellar disks showing O and Wolf–Rayet stars ([80, 101, 17, 135, 133]), as well as the isotropic S-star cluster with many B-and later-type stars (Gillessen et al. 2017).

These stars have been the subject of intense, decades-long studies. They have established the presence of the central supermassive black hole (SMBH) in the Galactic center ([52, 51, 54]), and the closest S stars can be used for stringent tests of general relativity (GRAVITY Collaboration et al. 2018; Do et al. 2019; GRAVITY Collaboration et al. 2020). The S stars can also be used to probe the extended mass-distribution in the Galactic center (GRAVITY Collaboration et al. 2020, 2022).

The origin of the S stars is a long-standing puzzle, considering that strong tides from Sgr A* are expected to suppress star formation in the S-star region ([118, 53]). Instead, the S stars may have migrated inward from larger scales via tidal disruption of binary stars ([71, 57]). Binaries may either come from relatively large scales ($\sim$ tens of parsecs; [107, 67]) or from one of the young disks ([86, 48, 111]). Remnants from binary stars would initially be on highly eccentric ($e \gtrsim0.98$; [71]) orbits. In contrast, the observed S stars have a thermal eccentricity distribution. However, other stars or remnants can perturb the stars’ orbits following binary disruptions. The stars’ energy only evolves due to uncorrelated two-body encounters (nonresonant relaxation). The star’s angular momentum can evolve on much faster timescales due to coherent torques (resonant relaxation; [73, 106, 78, 9] (Rauch & Tremaine 1996)). The S stars can easily isotropize within their lifetime via vector resonant relaxation ([73]), but the eccentricity distribution relies on the slower scalar resonant relaxation time. Many papers have studied whether this process can reproduce the observed eccentricities within the stars’ lifetimes.

For example, Merritt et al. [91] and Perets et al. [106] used N-body simulations to provide the first detailed theoretical study of the S stars’ eccentricity evolution due to an intermediate-mass black hole (IMBH) and stellar-mass black holes, respectively. Antonini & Merritt [9] were the first to study resonant relaxation in the Galactic center with both relativistic and Newtonian precession. Initial studies used approximate treatments for resonant relaxation. Subsequently, a self-consistent formalism was developed (Bar-Or & Alexander 2014; Sridhar & Touma 2016; Fouvry et al. 2017; Bar-Or & Fouvry

2018. and applied to the S stars by Generozov & Madigan [48] and Tep et al. [129].

Generally, these studies find that a population of stellar-mass black holes is needed to reproduce the S stars’ eccentricity within their lifetimes. The exact number and mass required depend on the modeling assumptions, the observational dataset, and the uncertain ages of the S stars. The other orbital elements provide complementary constraints on the S stars’ origin. For example, Perets & Gualandris [105] found that many theoretical models have difficulty explaining the populations of B stars on small and large scales.

We revisit the problem of the S stars’ orbital evolution, in light of the most recent S-star observations and constraints on the background mass distribution. For example, recently it was shown that there is a dearth of S stars with pericenters of

even though they would be observable ([22]). The S stars’ eccentricity distribution in the latest GRAVITY data is thermal, within observational uncertainties. The zone of avoidance constrains the semimajor axis distribution of the progenitor binaries, which depends on their intrinsic distribution and the rate of loss cone refilling. The eccentricity distribution of the S stars is sensitive to the properties of the background population, including the mass distribution of stellar-mass black holes.

We show that these orbital distributions can be explained via a combination of binary disruptions and relaxation in the Galactic center. If binaries are continuously disrupted in the Galactic center, the observed eccentricity distribution can be reproduced with a population of $\sim10M_{\odot}$ black holes in the S-star cluster. In the case of impulsive injection, as is expected in the young disk scenario, the background black holes would have to be at least $\sim30M_{\odot}$, considering existing constraints on the total enclosed mass in the Galactic center.

Our model accounts for many processes that to our knowledge have never been studied together in the literature. In particular, we have (i) an observationally motivated progenitor binary population, (ii) self-consistent diffusion coefficients for resonant relaxation, (iii) both eccentricity and semimajor axis evolution, (iv) stellar evolution, and (v) a loss cone. We also model the K-band luminosity function (KLF) of the S stars and incorporate this into our observational comparisons.

The remainder of this paper is organized as follows. In Sect. 2, we discuss our methods and the relevant physical processes. In Sect. 3, we show that most S-star observables can be reproduced by continuous disruption of binaries in the Galactic center. In Sect. 4, we discuss additional effects that are not accounted for in our fiducial models. In Sect. 5, we discuss models where binaries are all disrupted at a fixed lookback time. We summarize our main results in Sect. 6.

Relevant physical processes and methods

Binary disruptions

We assume that all the S stars are initially sourced by the disruption of binaries that come near their tidal radius:

where $M_{\mathrm{bh}}$ and $m_{\mathrm{bin}}$ are the mass of the black hole and binary, and $a_{\mathrm{bin}}$ is the binary semimajor axis. During a tidal interaction with an SMBH, a binary can (i) collide, (ii) survive (with some change to the orbital elements), or (iii) disrupt ([72, 56, 8, 20]). In the latter case, the binary is split into two stars. In the limit of a parabolic center-of-mass orbit, one of the stars is always bound to the SMBH, while the other is unbound. The semimajor axis of the remnant star is (e.g., Sari et al. [119])

where $m_{\mathrm{ej}}$ is the mass of the ejected star and $\chi$ is typically a factor of order unity. The pericenter of the remnant star, $r_{\mathrm{p,s}}$, is comparable to that of the progenitor binary. Thus, $r_{\mathrm{p,s}}\lesssim r_{\mathrm{t}}$. The eccentricity of the remnant star, $e_{\mathrm{s}}$, is close to one. More precisely,

Binaries can have a few different sources: the clockwise disk ([86, 48, 111]), scattering from massive perturbers such as large-scale molecular clouds or stellar clusters $\sim5$–$100$ pc away from the Galactic center ([107]), or large-scale gaseous spiral arms at similar scales Hamers & Perets [67]. We focus on the massive perturbers’ origin here.

To model binary disruptions in the Galactic center, we first constructed a model field-like binary population such that

1. The primary mass distribution follows a Kroupa mass function ($dN/dm_{\ast}\propto m_{\ast}^{-2.3}$ for $m_{\ast}>0.5M_{\odot}$) (Kroupa [79]).

2. The binary properties (semimajor axis, eccentricity, and mass ratio) follow the empirical distributions in Moe & Di Stefano [92]. We only considered mass ratios above 0.1, considering the distribution is unconstrained for smaller mass ratios. The higher multiplicity fraction of massive stars is accounted for.

In principle, dynamical processes such as binary evaporation can affect the binary distribution. In practice, this is not relevant for binaries that would populate the S-star region. Specifically, binaries would disrupt prior to evaporating or hardening (see Appendix A).

For each binary, we simulated a close encounter with a pericenter drawn from a uniform distribution (up to three times the maximum tidal radius of the binary population), implicitly assuming binaries are in the full loss cone regime as might be expected for the massive perturbers scenario. Furthermore, the S stars imply high rates of disruptions that are difficult to reproduce in the empty loss-cone regime ([107, 105]). If the pericenter of the encounter is greater than three times the binary tidal radius, the binary is unlikely to be disrupted (see [119] and Figure 5 in Generozov & Madigan [48]), and thus unlikely to deposit stars in the S-star region. Such binaries were excluded from further analysis. Otherwise, we simulated a close encounter between the binary and the central SMBH with the Fewbody code (Fregeau et al. [41]).

For simplicity, we assumed a constant disruption rate. Model binaries were disrupted at a random look-back time between $10^{9}$ years ago and the present day. Furthermore, the binary was disrupted at a random point along the main sequence of the primary star. We discuss the implications of a disruption rate rising toward the present day in Sect. 4, motivated by recent observational results on the star formation history in the Galactic center. Such models may be able to better reproduce the age distribution and the KLF of the S stars. However, more massive background black holes may be required to reproduce the S stars’ eccentricity distribution. Finally, we collected the remnants of all disrupted binaries with semimajor axes $< 0.05$ pc, as S-star progenitors, and modeled their subsequent orbital relaxation.

Relaxation

Binary disruptions deposit stars onto highly eccentric orbits (see equation (4)), and cannot account for the nearly thermal S-star eccentricity distribution Perets et al. [106]. However, the orbits of remnant stars can evolve due to nonresonant and resonant relaxation ([73, 106]). In the former, stars change their energy and angular momenta due to uncorrelated two-body encounters. In the latter, stars change their angular momenta due to coherent torques in potentials with a high degree of symmetry (like the nearly Keplerian potential in the Galactic center).

The nonresonant relaxation time at Galactocentric radius $r$ is

where $P(r)$ is the orbital period, $N_{\ast}(r)$ is the number of stars enclosed within $r$, $Q$ is the ratio of the SMBH mass to a characteristic stellar mass, and $\gamma$ is the power law index of the 3D density profile. It should be noted that the presence of stellar-mass black holes or other massive objects can significantly reduce the relaxation timescale (see review by Alexander [3]).

In the limit in which the precession is dominated by extended stellar mass, the resonant relaxation time is

For highly eccentric orbits, general relativistic precession dominates and resonant relaxation is suppressed (Merritt et al. [90]). These timescales depend on the background density profile of stars and remnants. We modeled the background with two components: (i) low-mass stars and remnants ($\sim1M_{\odot}$), and (ii) stellar-mass black holes ($\sim10M_{\odot}$). This is a reasonable approximation for modeling relaxation in an evolved galactic nucleus, and has been used extensively in the literature (Merritt et al. [90]). In principle, the high-eccentricity S stars should also contribute to the background density. In practice, we assume that the enclosed mass is dominated by the older, relaxed population that dominates the central $10^{\prime\prime}$ (Schödel et al. [122]).

Well inside the sphere of influence, a relaxed stellar profile will be between $r^{-1.5}$ and $r^{-1.75}$, depending on the relative abundance of stellar-mass black holes (Bahcall & Wolf [12]; Alexander & Hopman [4]). The black hole profile will be steeper than $r^{-2}$ on larger scales, and flatten to $r^{-1.75}$ in the innermost parts of the Galactic center ([42, 74, 4, 110, 1, 131, 136]). We used the following power law approximations for the density profiles of stars and black holes:

This is a reasonable approximation for different Galactic center models. (7) is within $\sim10\%$ of the density profiles from the Fokker-Planck models Vasiliev [131] and within $\sim40\%$ of the ‘Fiducial’ model in Generozov et al. [50] between 0.01 and 0.1 pc after 10 Gyr of evolution. The final density profiles of these models are similar, despite differences in the initial conditions. In Vasiliev [131], the stars and black holes all form together in the distant past. In Generozov et al. [50], the black holes are completely sourced by ongoing star formation at the present-day location of the clockwise disk. Recent measurements of the diffuse starlight and resolved stars in the Galactic center indicate a stellar density profile between $r^{-1.1}$ and $r^{-1.4}$ ([120, 46]). There is some evidence for stellar-mass black holes in the central parsec from observations of X-ray transients in this region (Mori et al. [93]). There is also a population of quiescent black hole X-ray binary candidates ([66, 94]), though the identification of these sources as black holes is controversial (Maccarone et al. [85]).

In this model, the total mass enclosed within the apocenter of S2 is $\sim1.5 \times10^{3}M_{\odot}$, comparable to the $1$-$\sigma$ upper limit from the GRAVITY collaboration ($\sim1000M_{\odot}$; Collaboration et al. [28]). We note that, in principle, the present-day upper limit does not rule out a larger number of black holes in the past.

We calculated angular momentum diffusion coefficients for this model, using JuDOKA (Tep et al. [129]). Then we evolved the remnants’ orbital elements forward in time, using the Monte Carlo procedure described in Bar-Or & Alexander [15]. For every timestep, $\Delta t$, the reduced angular momentum, $j^{2}$, and the semimajor axis, $a$, are incremented by

respectively, where $D_{j}$, $D_{jj}$, $D_{a}$, and $D_{aa}$ are the diffusion coefficients, and the subscripts RR and NRR denote resonant and non-resonant relaxation, respectively. The unprimed $\gamma$ coefficients were independently drawn from a Gaussian distribution of unit variance, and

accounting for the covariance between energy and angular momentum diffusion. We used an adaptive timestep for the evolution:

where $\max(D_{jj,\mathrm{RR}}(j))$ is the maximum of $D_{jj,\mathrm{RR}}$ at fixed energy. Gravitational wave emission was included, following Peters [108], but was negligible except for a small minority of stars (e.g., for 99% of stars, the semimajor axis changes by less than 1% with the semimajor axis diffusion artificially turned off). To speed up our calculations, we computed the diffusion coefficients for a predefined grid of semimajor axes and angular momenta.

At each step, we linearly interpolated the resonant diffusion coefficients in $a$ and $j$. The nonresonant diffusion coefficients also depend on the S-star masses. More precisely, they are a linear combination of integrals of the background distribution function, with mass-dependent coefficients. We linearly interpolated each integral term in $a$ and $j$ and combined them. Our interpolation tables extend between $j = 0.001$ and $j = 0.999$ and between $a = 7 \times10^{-4}\ \mathrm{pc}$ and $0.5\ \mathrm{pc}$. These are the boundaries of our integration domain.

The outer boundaries are purely reflective. The inner semimajor axis is purely absorptive.

At the inner $j$ boundary, stars can be disrupted if their pericenter falls below the stellar tidal radius (Rees (1988)),

Here, $m_*$ and $r_*$ are the stellar mass and radius, respectively. These were evolved as is described in Sect. 2.3, so that the inner boundary depends on stellar properties.⁴ To avoid stars overshooting and crossing into negative angular momenta, we implemented a second, reflective inner boundary at $j = 10^{-3}$.

The true tidal radius may differ by a factor of order unity from equation (12), based on the results of hydrodynamic simulations (Ryu et al. [116]). We used equation (2) in any case, as there is currently no fitting formula for the tidal radius of evolved stars that is calibrated to hydrodynamic simulations.

Stellar evolution

As time progresses, the stars evolve. We used both PARSEC isochrones ([21, 26, 25, 128, 88, 100, 99]) and MIST ([35, 27, 102, 103, 104]) stellar evolution tracks⁵ to evolve the stellar masses, radii, and types. We find similar results with both approaches. For conciseness, we only present the results from MIST.

We present results for solar metallicity, though we find similar results for [Fe/H]$=0.5$. There are no existing metallicity constraints for the young stars in the Galactic center.

Summary of model

We now summarize our forward model for the formation of the S-star cluster, starting from binary disruptions. The key steps are (see Figure 1):

1. Generate model binaries following the assumptions in Sect. 2.1. (a) We only consider binaries with semimajor axes up to 15 au, as wider binaries are unlikely to deposit stars into the S-star region (see equation (3)). (b) The primary mass is between 1 and $100M_\odot$. Lower-mass binaries would not produce the presently observed S stars.

2. Generate close encounters with the SMBH following a Monte Carlo procedure: (a) The look-back time for the encounter is sampled uniformly between $10^9$ yr ago and the present day. (b) The pericenter of the encounter is sampled uniformly between 16 gravitational radii (∼ the tidal radius of a single star) and three times the maximum tidal radius of the binary population (equation 2). Thus, we implicitly assume binaries are in the full loss cone regime, such that the pericenter of the center-of-mass orbit is uniformly distributed (see Merritt [89] for a discussion of loss cone physics). (c) If the pericenter of the encounter exceeds three times the binary tidal radius (equation 2), the binary is unlikely to disrupt, and we proceed to the next one ([119, 48]).

3. We simulate tidal encounters between each binary and the SMBH using the FEWBODY code (Fregeau et al. [41]). In these encounters, the binary is initialized at 50 tidal radii on a nearly parabolic orbit with respect to the SMBH. The binary is generally split into two angular momentum bound and unbound stars (see Sect. 2.1).

(a) We discard cases where the binary stars collide. This happens only $\sim1\%$ of the time.

4. We collect all bound remnant stars with semimajor axes $< 0.05\ \mathrm{pc}$, and follow their angular momentum and energy relaxation in a fixed background, as is described in Sect. 2.2. The effect of the remnant stars on the background density profile is neglected [40].

5. Stars are allowed to evolve and die (see Sect. 2.3). Furthermore, stars may be tidally disrupted by the central SMBH (see equation (12)).

6. We collect the properties of all surviving stars with semimajor axes $< 0.05\ \mathrm{pc}$ and $K$-band magnitudes $< 18$, and compare them to observations.

7. For each model star, we compute the probability to measure its orbital elements according to equation 1 in Burkert et al. [22]. To account for observational bias, we exclude stars according to this probability.

Results

The initial conditions of bound stars with semimajor axes of $a_s \le0.05\ \mathrm{pc}$ are shown in the top panel of Figure 2. As was expected, all stars start on highly eccentric orbits. The semimajor axis distribution is close to uniform, due to the approximately log-uniform semimajor axis distribution of the progenitor binaries, and our assumption of a full loss cone (see also Perets & Gualandris [105]). The minimum semimajor axis is $\sim5 \times10^{-4}\ \mathrm{pc}$, corresponding to remnants of near-contact binaries (see equation (3)). After $\sim200\ \mathrm{Myr}$, a quasi-steady state is reached in the S-star region. The final pericenters and eccentricities of surviving stars are shown in the bottom panel of Figure 2. Purple and red stars correspond to early- and late-type stars, respectively.⁶ We also show the observed orbits of early and late-type S stars with blue squares and orange circles, respectively. We group observed stars of unknown type with early stars. If they were late-type stars, they would likely have been identified as such, considering CO absorption lines are prominent, and therefore easily observed with current instruments. Unless otherwise specified, here and in other observational comparisons, we exclude (i) disk S stars⁷, (ii) unbound S stars, and (iii) S stars with semimajor axes greater than $1.25^{\prime\prime}$. The disk and the rest of the S stars differ in their eccentricity distribution, suggesting these are distinct populations that formed separately. The presence of O and WR stars in the disk also supports this. We exclude stars at large semimajor axes, as they will be strongly affected by completeness effects, complicating any comparisons to our model. Furthermore, our assumption of a full loss cone will break down at large radii (see Perets & Gualandris [105]), further complicating the comparison.

Overall, we find a reasonable agreement with the zone of avoidance. Approximately $\sim14\%$ of model stars are in the zone of avoidance, as it is defined in Burkert et al. [22]. Similarly, in the latest orbital data, 3 of the $\sim40$ observed stars are in the zone of avoidance (S14, S18, and S23). If we adjust the zone of avoidance so that all observed stars are outside (dashed gray line in Figure 2), then only $\sim3.9\%$ of model stars are inside the modified zone. For 40 stars, this would correspond to a $\sim20\%$ probability of no stars being in the zone of avoidance.

Figure 3 shows the one-dimensional pericenter and eccentricity of our model stars, compared to the observed S stars. Both distributions are consistent with the observations.

Figure 4 shows the orbital properties of early S stars with $K \le16$ and $K > 16$. This magnitude roughly divides the observed sample in half. The model eccentricity and pericenter distributions of each group are consistent with observations, though the bright stars have a somewhat superthermal eccentricity distribution. The faint stars’ semimajor axis distribution is too skewed toward large semimajor axes, compared to observations. In part, this is due to the assumed model for the progenitor binaries: as the primary mass decreases, the period distribution gradually shifts from approximately log-uniform to log-normal (with a peak near $10^5$ days). We also note that faint stars at large semimajor axes are the most difficult to detect, and this comparison is subject to observational bias.

We now compare the stellar properties of model stars with observations. In the model, $\sim20\%$ of stars are late-type. This is consistent (within Poisson uncertainties) with the observed fraction of $\sim26\%$. However, we note that we have not included disruptions of old ($\gtrsim10^{10}$) stars that make up the bulk of the nuclear stellar disk, and that can increase the population of evolved stars in the S-star region.

We also compare the K-band magnitudes of our model stars to observations (Figure 5). In practice, we used the JWST F210M magnitude from MIST isochrones ([35, 27, 102, 103, 104]) to approximate the K-band magnitude of our model stars. We also assumed an extinction of 2.42 ([44]) and a Galactocentric distance of 8.3 kpc ([61]). We neglected corrections from differential extinction ($\sim0.3$ mag; [121] (Do et al. 2013)).

The top panel of Figure 5 compares the overall KLFs from theory and observations. We rescaled the luminosity functions to match at $K = 16.2$. There are significant discrepancies between the model and observed S-star KLF at both the faint and bright ends. The former can plausibly be explained by completeness effects. The latter is due to $\sim3$ bright giants with $K < 14$ that are not observed, and is of marginal significance. We note that collisions with black holes cannot remove these bright giants (see Sect. 4.2).

The bottom panel of Figure 5 compares the model and observed KLF of early-type stars (blue and orange). Here, there are no large differences at the bright end, while the model has many more stars with $K > 16$. This discrepancy may again be due to completeness effects. Unfortunately, this is difficult to quantify, as the observational selection function is highly nontrivial (the orange KLF only includes stars with measured orbits). The model more closely matches overall KLF from Schödel et al. [122] (light red points; see Appendix B for details), which should be more complete, as even stars without orbits are included.

However, the progenitor binaries are more likely to follow a present-day mass function than an initial mass function (as assumed here). A present-day mass function would steepen the model KLF, and exacerbate the tension with observations (see Sect. 4).

The top panel of Figure 6 shows the K-band magnitude as a function of age for the eight bright ($K \le15.5$) S stars from Habibi et al. [65] and for our model stars, which appear significantly older. To better compare the age distributions, we selected a subset of model stars with comparable K magnitudes. First, we selected all model stars within 0.1 mag for each of the Habibi et al. [65] stars. Then, we sampled 1000 stars (with replacement) from each of these subsets. The gray line in the bottom panel of Figure 6 shows the stacked age distribution of the selected stars. On average, model stars are significantly older, even after accounting for the observational bias toward bright stars. This may suggest that the brightest S stars have a distinct origin or that the star formation rate (and delay time distribution between star formation and disruption) is nonuniform. Finally, there is a degeneracy between the age and metallicity (see Appendix § 9), such that the S stars may in fact have subsolar metallicity and be significantly older. Currently, there are no measurements of metallicity for the S stars, but future observations of the S stars with ERIS and NIRSpec will constrain this. A priori, this explanation is less likely, considering the nuclear stellar disk and nuclear star cluster are mostly metal-rich with small metal-poor populations ([37, 113, 96, 125, 122, 124, 97, 98] (Do et al. 2015)).

Finally, Figure 7 shows a corner plot of $a$, $e$, and $K$. Once again, we see a good agreement between the observed and theoretical eccentricity distributions. At the same time, the model produces too many faint stars, and too many stars at large semimajor axes (though again we caution that faint stars at large semimajor axes are the most difficult to detect, and the latter discrepancy may be due to observational bias). In the observations, many of the lowest-eccentricity stars are in a narrow band of semimajor axes ($3.3 \lesssim\log(a/\mathrm{au}) \lesssim3.55$). This feature is not apparent in the model. However, there are not enough stars to conclude that this is a statistically significant difference. Furthermore, the absence of this feature may simply reflect the differences in the one-dimensional semimajor axis distributions. To test this, we selected the model star with the closest semimajor axis for each observed S star. The structure of these model stars in the $a-e$ space is more similar to the observed one, as is shown in the bottom panel of Figure 7. Interestingly, there is a band of low-eccentricity stars near the Galactocentric radius where the scalar resonant relaxation time⁸ is minimized ($\log(a/\mathrm{au}) \approx3.24$). This is where angular momentum relaxation is most efficient.

Critical discussion of model assumptions

Here, we critically examine a few assumptions about the progenitor binary population, and how they affect observational comparisons. Firstly, the primary mass of each binary is drawn from an initial mass function. However, a steeper present-day mass function is more realistic, considering the shorter lifetimes of massive stars. Secondly, binaries disrupt at a constant rate and are allowed to disrupt immediately after formation. In reality, binaries may need some time to reach disruption, especially if they form on circular orbits. As is discussed below, delays and a present-day mass function exacerbate the tension with the observed KLF.

We also discuss the effect of stellar collisions and repeated, individually nondestructive tidal encounters between the stars and the central SMBH. These effects are not included in our base model, but do not affect the results.

Mass function and disruption times

In our base model, the primary mass of each binary is drawn from a Kroupa initial mass function. However, the mass function is of disrupting binaries will likely be more bottom-heavy, due to the longer lifetimes of low-mass stars.

Here, we account for this effect by re-weighting probability distributions. In particular, each remnant star is weighted by the main-sequence lifetime of the primary star in the progenitor binary. This causes the KLF to steepen, in tension with observations, as is shown in Figure 8.

This result depends on the star formation history in the Galactic center, particularly within the region of the nuclear stellar disk, where disrupting binaries likely originate (see Figure 3 in Perets et al. [107]). Schödel et al. [123] recently showed that the star formation history in the nuclear stellar disk is highly nonuniform: $\gtrsim70\%$ of the stars formed more than 10 Gyr ago, $\sim15\%$ formed $\sim1$ Gyr ago, and up to 10% formed in the last tens to hundreds of millions of years. This suggests a star formation rate rising toward the present day (since a comparable number of stars were formed within the last $\sim10^8$ yr and $\sim10^9$ yr). Furthermore, the star formation rate is likely correlated with the binary disruption rate, since both are correlated with the molecular gas density.

Motivated by these observations, we discarded stars that formed more than $10^8$ yr ago to see the effect of a star formation history rising to the present day. Figure 9 shows the KLF after this cut, with and without re-weighting for stellar lifetimes. In the former case, massive stars receive a weight equal to their main-sequence lifetime divided by $10^8$ yr. Stars with main-sequence lifetimes longer than $10^8$ yr receive a weight of unity. This ameliorates, but does not entirely remove, the tension with observations, as is shown in Figure 9. Furthermore, the model eccentricity distribution is in tension with observations, though this can be resolved if the background black holes are each $20\,M_\odot$ instead of $10\,M_\odot$.

The results also depend on the delay time distribution between star formation and disruptions. In the base model, binaries can be disrupted immediately after formation. More realistically, binaries will take some time to reach the loss cone, especially if they form on circular orbits. Imposing a minimum age at disruption tends to steepen the KLF, as is illustrated in Figure 10. This shows the KLF after (i) discarding stars older than $10^8$ yr and (ii) discarding stars injected into the Galactic center less than $10^7$ yr after formation.

The better agreement between theory and observation in Figure 9 suggests that the early S stars are primarily sourced by binaries formed within the last $\sim100$ Myr. Nonetheless, a slight

tension remains, especially for realistic delay times between star formation and disruption, as in Figure 10. In this case, there is a $\sim2\sigma$ tension between the observed (orange) and theoretical (turquoise) KLF between $K = 14$ and $K = 16$.\ ⁹ Potentially, this tension could be resolved via a bursty star formation rate that is more skewed toward the present day. This may also help with the tension between model and observed ages among the bright S stars (see Figure 6).

Alternatively, the tension with the observed KLF may be resolved with a top-heavy initial mass function. The mass function required itself depends on the binary injection history. In the last case (Figure 10), an initial mass function of $dN/dm \propto m^{-0.2}$ reduces the tension to $1\sigma$.

Collisions

In principle, stars may suffer mass loss or destruction via collisions with other stars or compact objects. Alternatively, star-star collisions can lead to mergers and the buildup of massive stars ([126, 36, 95, 52, 30, 109, 64, 43, 29, 114]).

We estimated the impact of destructive black-hole star collisions. In particular, we computed the orbit-averaged collision rate between the remnant stars and black holes in our Monte Carlo models:

where $m_{\ast}$ and $m_{\mathrm{bh}}$ are the stellar and BH masses, respectively, $r_{\ast}$ is the stellar radius, $v(r) \approx\sqrt{GM_{\mathrm{bh}}/r}$ is the local velocity dispersion, and $n_{\mathrm{bh}}(r)$ is the BH density (see Eq. (7)). The orbit average was calculated from $10^{4}$ points along the orbit, equally spaced in eccentricity anomaly ($EA$). Each point is weighted by $(1-e\cos(EA))$ in the average.

For each Monte Carlo timestep, $\delta t$, we generated a random number between 0 and 1, and recorded a collision if it was less than the expected number of collisions, $\langle\Gamma_{\mathrm{coll}}\rangle\delta t$.

Removing stars that experience collisions has a minor impact on the results. Most notably, the fraction of late-type S stars decreases from $\sim20\%$ to $\sim16\%$. In fact, the impact of collisions may be overestimated, considering that deeply penetrating encounters are necessary to significantly dim giant stars (Dale et al. [29]).

Rose et al. [114] have shown that star-star collisions in the S-star region can lead to the buildup of massive stars (but this would require nondisruptive collisions). We leave this effect to future studies, though it may have a major impact on the comparisons to the observed KLF and ages. Additionally, stellar collisions can lead to transients that look like tidal disruption events or supernovae ([13, 7, 115, 31]).

Repeated tidal encounters with the SMBH

So far, we have only considered the removal of stars due to full tidal disruptions by the central SMBH. However, stars can also be removed due to multiple partial tidal disruptions. Furthermore, stars can become swollen and disrupted due to tidal heating ([5]). To experience such effects, the star needs to pass within a factor of a few of the central SMBH ([63, 87, 117, 47] (Li & Loeb 2013; Lu et al. 2021)). To test how these effects would change the results, we multiplied the tidal radius of each star by two. We found that this has no significant effect on the comparisons previously discussed.

Fixed age model (eccentric disk)

So far, we have considered models in which stellar ages and disruption times can vary freely. We now discuss models in which all stars are injected at similar times in the recent past. Such an impulsive burst of disruptions may be triggered by a secular gravitational instability in the $\sim2.5$–$8$ Myr old ([80, 101, 17, 135] (Lu et al. 2013)) clockwise disk if it starts in a lopsided, eccentric configuration ([86, 48, 111]). This configuration can arise from the tidal disruption of a molecular cloud ([19, 49] (Gualandris et al. 2012)), or via an SMBH merger in the Galactic center ([2]).¹⁰ Interactions between the misaligned young disks can also trigger a burst of binary disruptions (Löckmann et al. [82]).

Previously, Generozov & Madigan [48] and Tep et al. [129] analyzed the relaxation of the S-stars in such scenarios. The former found that $10M_{\odot}$ BHs would be able to reproduce the observed S-star eccentricity distribution, by modeling angular momentum relaxation at a fixed semimajor axis of $0.01$ pc. However, the resonant relaxation time will be significantly longer in the outer part of the S-star cluster. Tep et al. [129] accounted for the S stars’ semimajor axis distribution and found that reproducing the observed S-star eccentricity distribution requires black holes above $100M_{\odot}$ at $2\sigma$ confidence, though this constraint depends on the S-star sample (i.e., this is the constraint including late-type and disk S stars).

We performed a similar analysis, assuming that stars are injected into the Galactic center all at once. First, for consistency with Tep et al. [129], we turned off (i) the loss cone, (ii) stellar evolution, and (iii) semimajor axis diffusion. Effects (i) and (ii) will be negligible for the S stars over million-year timescales. However, we show that the loss cone can affect black hole mass constraints by a factor of $\sim2$.

We also matched the initial conditions in Tep et al. [129], assuming that all stars were injected into the Galactic center 7.1 Myr ago, the mean age of eight bright S stars from Habibi et al. [65]). The initial angular momentum distribution is a Gaussian with a mean of $0.2$ ($e=0.98$) and a standard deviation of $0.02.

We present results for both the density profiles in Tep et al. [129]:

and our fiducial density profile (Eq. (7)), which is in better agreement with the latest constraints from GRAVITY (Collaboration et al. [28]) and which has a factor of 3 fewer BHs within 0.05 pc.

Using the same S-star sample and background profile as Tep et al. [129], we find good agreement with their results: black hole masses $<100M_{\odot}$ can be ruled out. However, the full sample includes late-type stars and stars at larger scales that may be unrelated to binary disruptions. For example, eight of the S stars are members of the clockwise disk. For the 25 early, non-disk S stars with semimajor axes within 0.05 pc, only black holes with masses $<17M_{\odot}$ can be ruled out. For the same subsample and our fiducial density profile, masses $<50M_{\odot}$ can be ruled out. Finally, the lower limit is reduced to $\sim24M_{\odot}$ with a loss cone.

Overall, explaining the entire S-star population via binaries from the young clockwise disk requires a background with more massive black holes. There is also tension between the mass functions of the disk and S stars that may require an unusual mass ratio distribution within the disk to reconcile (Generozov [47]). Moreover, it is not clear whether the conditions for the eccentric disk instability could arise with steep background density profiles (e.g., Eq. (7); see the discussions in Madigan et al. [86] and Generozov et al. [49]). On the other hand, efficient angular momentum relaxation of the S stars likely requires a steep black hole cusp. Detailed modeling of disk binary disruptions suggests that they can reproduce the zone of avoidance, but remnant stars will be skewed toward large semimajor axes compared to the observations.¹¹

Summary and discussion

We show that many properties of the observed S-star cluster in the Galactic center, including the recently discovered zone of avoidance, can be explained by a combination of isotropic binary disruptions, and relaxation in the Galactic center. Our focus is on explaining the classic, isotropic S-star cluster. In our model, the young disks would have a different origin that we leave to future work.

Our main conclusions are summarized as follows.

1. In the full-loss-cone regime, the remnant stars from massive binary disruptions will have an approximately uniform semimajor axis distribution. This naturally explains the lack of S stars in the zone of avoidance. Empty loss cone models would place too many stars in this region.

2. The overall orbital eccentricity and pericenter distribution of the S-star cluster is well reproduced in our models, via relaxation processes. We have only considered nonresonant and resonant relaxation from a background population of stars and $10M_\odot$ black holes. Potentially, an AGN disk may also perturb the S-star orbits, dramatically reducing the angular momentum relaxation time [24].

3. There is a slight tension with the observed semimajor axis distribution of fainter S stars ($K < 16$), potentially indicating the progenitor binary population is more skewed to small separations than in Moe & Di Stefano (2017).

4. The model KLF is in tension with observations, particularly with a Kroupa present-day mass function and realistic delay times between star formation and disruption. Potentially, this tension may be resolved with a top-heavy initial mass function or with a bursty star formation rate that is skewed toward the present day. However, reproducing the observed thermal eccentricity distribution would be more challenging in such scenarios.

5. There is also tension with the measured ages of the bright S stars, with the model predicting older ages. Again, this may point to a bursty star formation rate. Alternatively, the measured ages may be underestimated if the S stars have subsolar metallicities.

Acknowledgements. We thank the anonymous referee for helpful comments. We are very grateful to our funding agencies. AG was supported at the Technion by a Zuckerman Fellowship. RS acknowledges financial support from the Severo Ochoa grant CEX2021-001131-S funded by MCIN/AEI/10.13039/501100011033, grant EUR2022-134031 funded by MCIN/AEI/10.13039/501100011033 and by the European Union NextGenerationEU/PRTR, and grant PID2022-136640NB-C21 funded by MCIN/AEI/10.13039/501100011033 and by the European Union. We thank Pau Amaro Seoane, Andreas Burkert, Sebastiano von Fellenberg, Jean-Baptiste Fouvry, Karamveer Kaur, Fabrice Martins, Nicholas Stone, Odelia Teboul, and Kerwann Tep for helpful comments and discussions. This work made use of the following software packages: Astropy (Astropy Collaboration et al. 2013; The Astropy Collaboration et al. 2018; Astropy Collaboration et al. 2022), Fewbody (Fregeau

Binary evaporation

In principle, the binary population in the Galactic center region may be shaped by dynamical processes. For example, soft binaries can be evaporated by encounters with other stars. Here we quantify this effect and show that in our models binaries are more likely to disrupt rather than evaporate, and thus binary evaporation may be neglected.

The local evaporation timescale is (see [see 6])

where $v_{12}$ is the internal velocity of the binary, $n$, $\sigma$, and $\langle m^{2}\rangle$ are the number density, velocity dispersion of the perturbing population, and $\Lambda_{12}=2\sigma^{2}/v_{12}^{2}$. In our model, binaries come from outside the central 5 pc, where $n\approx5\times10^{3}\,\mathrm{pc}^{-3}$, and $\sigma\approx100\,\mathrm{km\,s^{-1}}$. For binaries with semimajor axes $\leq15$ au in this region, the evaporation timescale is of order 10 Gyr or more.

In comparison, the timescale for binaries to disrupt is

where $j_{\mathrm{lc}}=\sqrt{2GM_{\mathrm{bh}}r_{\mathrm{t}}}$. For an enclosed mass of $10^{7}M_{\odot}$ inside 5 pc,

Binaries on eccentric orbits will cross the central region, where the evaporation timescale is shorter. In particular, the evaporation timescale will be between $r^{-1.3}$–$r^{-1.5}$.{^12} The average evaporation time over the binary population is

where $r_{\mathrm{p}}$, $r_{\mathrm{a}}$, $v_{\mathrm{r}}$, $P$, $e$, and $f(e)$ are the pericenter, apocenter, radial velocity, orbital period, eccentricity, and eccentricity distribution respectively. Note that the inner integral corresponds to a harmonic average.

In steady state $f(e)$ is expected to be thermal. Also, since we are considering binaries that are not too far outside the SMBH sphere of influence, we assume Keplerian scalings with radius for all the dynamical quantities. Overall, we find eccentric orbits give a $\sim10\%$ correction to the evaporation time.

Although we have focused on soft binaries, the hardening rate of hard binaries has a similar dependence on cluster and binary properties [69, 70]. Therefore, we also expect hard binaries to disrupt before they can significantly harden.

K-band luminosity function

Here, we describe how we derive the KLF of early stars inside $0.5^{\prime\prime}$ (the red line in Figure 5) from the data in Schödel et al. [122].

First, we compute the KLF in two radial rings with $r_{\mathrm{proj}}<0.5^{\prime\prime}$ and $4^{\prime\prime}<r_{\mathrm{proj}}<5^{\prime\prime}$. These KLFs include contribution from both early and late stars, with surface density profiles $r^{-0.8\pm0.08}$ and $r^{0\pm0.08}$, respectively (Table 6 in Gallego-Cano et al. [45]). The surface densities are equal at $\sim0.7^{\prime\prime}$, such that late (early) stars dominate on larger (smaller) scales.

The KLF of early stars within $0.5^{\prime\prime}$ is

where $K_{0.5}$ is the total KLF inside $0.5^{\prime\prime}$ and $K_{4}$ is the KLF between $4^{\prime\prime}$ and $5^{\prime\prime}$, which is dominated by late stars. The weight, $W$, is chosen to reproduce the ratio of late to all stars inside $0.5^{\prime\prime}$, $X$. From the surface density profiles in Gallego-Cano et al. [45], $X\approx0.3$ and

where the numerator is the number of stars with projected radius, $r_{\mathrm{proj}}$ less than 0.5″ and $K \le16$. Similarly, the denominator is the number of stars with 4″ < $r_{\mathrm{proj}}$ < 5″ and $K \le16.15$. The cut near $K = 16.15$ is motivated by exclusion of fainter stars from the surface density fits in [45].

From equations (B.1) and (B.2), we obtain the $K_{\mathrm{early},0.5}$, the red line in Figure 5. The uncertainty is calculated by adding the uncertainties from $K_{0.5}$ and $K_4$ in quadrature. The results change by less than the uncertainty, if we use the KLF from 2–3″ or 3–4″ in lieu of $K_4$. In fact, the difference between $K_{0.5}$ and $K_{\mathrm{early},0.5}$ are smaller than the uncertainties for $K \lesssim19$.

S-star ages versus metallicity

Here we show that the S-star ages are degenerate with metallicity. Figure C.1 shows a comparison of MIST isochrones of $\log(g)$ versus $\log(T_{\mathrm{eff}})$ for two different metallicities. Qualitatively, the S stars can either be younger and metal-rich or older and metal-poor. For each star, we estimate the age by identifying the closest main-sequence isochrone point.¹⁴ We summarize these estimates in Table C.1, along with the fits from Habibi et al. [65]. We emphasize that our ages are merely estimates to illustrate the trend with metallicity, rather than detailed fits.

Star	Habibi et al. 2017 [Myr]	MIST, [Fe/H]=0.0 [Myr]	MIST, [Fe/H]=-0.5 [Myr]	MIST, [Fe/H]=-1.0 [Myr]
S1	$4.3^{+4.2}_{-4.2}$	5.0	10.0	14.1
S2	$6.6^{+3.4}_{-4.7}$	4.5	7.9	11.2
S4	$5.9^{+4.1}_{-5.5}$	5.0	10.0	14.1
S6	$11.5^{+7.8}_{-9.3}$	8.9	17.8	25.1
S7	$0.3^{+14.3}_{-0.3}$	4.0	15.8	25.1
S8	$3.1^{+4.4}_{-3.0}$	3.2	7.9	11.2
S9	$3.5^{+13.2}_{-3.5}$	5.0	15.8	25.1
S12	$14.7^{+10.6}_{-13.9}$	14.1	28.2	35.5