© The Authors 2025

1 Introduction

The Galactic Centre (GC) is unique for its combination of relative proximity and extreme stellar dynamics. It is, however, also a challenging environment to study due to high line-of-sight extinction and source crowding. The GC is host to the super-massive black hole Sagittarius A* (Sgr A^*) with a mass of 4.3 × 10⁶ M_⊙ (Eisenhauer et al. 2005; Ghez et al. 2008; Schödel et al. 2009; Gillessen et al. 2009; Boehle et al. 2016; Gillessen et al. 2017; Do et al. 2019; GRAVITY Collaboration 2019, 2024) and a nuclear star cluster (NSC). In addition, Sgr A* is orbited by a cluster of mainly B-type main-sequence stars out to ∼0.04 pc called the S-star cluster (see review by Genzel et al. 2010). The origin of this cluster remains unknown, since the strong tidal force from Sgr A* inhibits standard star formation from molecular clouds in this region (Genzel et al. 2010). Nonetheless, the S-star cluster has been intensely studied, due to the possibility of tracing individual orbits of stars and thereby constraining the mass and distance to Sgr A* (e.g. Ghez et al. 2008; Gillessen et al. 2009; Do et al. 2013; Gillessen et al. 2017; GRAVITY Collaboration 2018).

Within the highly complex environment of the GC, exotic gravitational interactions can occur. Of particular interest to this work is the Hills mechanism, in which a stellar binary approaches a massive black hole (MBH) to within the tidal radius where the tidal force of the MBH overcomes the self gravity of the binary (Hills 1988; Yu & Tremaine 2003). The result is that one star is captured into a tight elliptical orbit, while the other is ejected as a hypervelocity star (HVS). HVSs can travel at ≳ 10³ km s⁻¹, making them unbound to the Galactic potential (Hills 1988; Bromley et al. 2006; Rossi et al. 2014). The captured stars on the other hand have been suggested as a possible formation mechanism for the S-star cluster, since they should get deposited at radii similar to those of the S-stars (e.g. Gould & Quillen 2003; Ginsburg & Loeb 2006). In addition, recent observations have provided strong indirect evidence that the Hills mechanism operates in the GC (Koposov et al. 2020).

In this study, we examine whether the observed properties of both the S-star cluster and HVSs are consistent with predictions by the Hills mechanism. This will tell us whether the entire S-star cluster might have formed through the Hills mechanism or whether a more complex assembly history is required to explain the observed properties of the S-stars. What makes this study unique compared to earlier studies (e.g. Madigan et al. 2009; Zhang et al. 2013; Generozov & Madigan 2020) is that we simultaneously simulated both the S-star and HVS populations and compared both with state-of-the-art observations.

In Section 2, we summarise our observational knowledge of the S-star cluster. Section 3 describes our methods for forward modelling the population of S-stars and HVSs. In Section 4, we investigate two possible origins for progenitor binaries. The discussion on our results and potential future work are presented in Section 5. Finally, in Section 6 we give our conclusions.

2 The S-star cluster

We start by giving an overview of the observed properties of the S-star cluster in Section 2.1, before discussing the effects of observational selection in Section 2.2.

2.1 Overview

The S-star cluster has isotropically distributed stellar orbits with a radial extent of roughly ∼0.04 pc from Sgr A^*; it is formed by mainly B-type massive main-sequence stars, and it is surrounded by at least one disc of massive young stars containing many O/WR stars out to ∼0.5 pc (Levin & Beloborodov 2003; Paumard et al. 2006; Bartko et al. 2009; von Fellenberg et al. 2022). No O/WR stars are detected in the inner 0.04 pc (Genzel et al. 2010).

The S-star cluster has been observed for nearly three decades (Eckart & Genzel 1996; Ghez et al. 1998). During this time, a range of instruments have been used to study its properties with techniques including adaptive optics-assisted astrometry, integral field spectroscopy, and optical/infrared interferometry. Due to the high line-of-sight extinction (A_V ∼ 50; Fritz et al. 2011), the study of the S-star cluster is done around 2 μ m, where extinction is lower (A_Ks ∼ 2.5; Schödel et al. 2010; Fritz et al. 2011).

The number of parameters needed to fully describe an orbit around the GC is six. Usually, these are sky position, proper motion, radial velocity, and acceleration in the plane of the sky (Gillessen et al. 2017). A catalogue containing orbits for 40 stars is published in Gillessen et al. (2017), most of which are in the central arcsecond. Furthermore, a catalogue of 36 young stars with orbital solutions is given in von Fellenberg et al. (2022). More recently, orbital parameters were determined for 20 stars at 2–7 arcsec from Sgr A* in Young et al. (2023).

Observations indicate that S-stars follow a thermal eccentricity distribution (n(e) ∝ e) and travel on orbits with periods as short as decades, with the shortest observed period being about 13 yr (Gillessen et al. 2017). The most studied S-star is S2: a bright (m_K = 14, without extinction correction), massive star (M_ZAMS ∼ 14 M_⊙) on a 16 year orbit (e.g. Habibi et al. 2017). The metallicity of the S-stars is currently unknown (Habibi et al. 2017).

The young disc surrounding the S-star cluster, often referred to as the clockwise disc (CWD), has a top-heavy initial mass function (IMF) with a slope of around −1.7 (Lu et al. 2013; Gallego-Cano et al. 2024). The CWD is consistent with a single star formation episode about 2–6 Myr ago (Paumard et al. 2006). There is also evidence for additional discs at varying radii (e.g. von Fellenberg et al. 2022).

In Fig. 1, we show orbital solutions in the 2D space of eccentricity against semi-major axis, for stars belonging to either the CWD or the S-star cluster (data from Gillessen et al. 2017). Our figure clearly shows that these two stellar structures occupy different ranges of this parameter space.

2.2 Observational selection function

The stars for which orbital parameters have been determined cannot be seen as an unbiased dataset for the properties of the S-star cluster: accelerations are more difficult to measure for faint stars with relatively long orbital periods. To understand this bias, we need a prescription of the observational selection function. Using mock observations, Burkert et al. (2024) finds an empirical relation that connects a star’s pericentre distance and eccentricity to the chance that an orbital solution would have been found in Gillessen et al. (2017). However, Burkert et al. (2024) did not consider the apparent magnitude of a star to be part of the selection function. Since they base their selection on our ability to measure the acceleration in the plane of the sky within their mock observations, it is implicitly assumed that a radial velocity measurement was also available. To account for this additional factor, we assumed the orbital solutions to be complete down to Ks = 16, based on the detection limit of SINFONI for a typical run (Gillessen et al. 2009). Since the Gillessen et al. (2017) catalogue, only a single additional orbital solution for a Ks < 16 star within the Burkert et al. (2024) completeness limit has been published (S5; von Fellenberg et al. 2022). This star was reported as being consistent with the CWD, but given its high eccentricity and small semi-major axis (typical for S-stars) we instead treated it as an S-star. For Ks < 16, there are respectively twelve and four orbital solutions for mainsequence and evolved S-stars in our dataset. We did not consider S111, since it has a hyperbolic orbital solution (Gillessen et al. 2017). Throughout this work, our stellar magnitudes are not corrected for extinction.

Fig. 1 Orbital solutions for stars at the GC found in Gillessen et al. (2017), with the exception of S5, which is from von Fellenberg et al. (2022). Stars that are considered part of the young disc are coloured blue, which are S66, S67, S83, S87, S91, S96, S97, and R44. S111 is omitted, because of the hyperbolic solution. The grey points show S-stars fainter than our completeness limit, while black points show those that are brighter.

3 Method

Now that we have discussed the properties and observations of the S-star cluster, we will describe how the Hills mechanism could explain the formation of both the S-star cluster and HVSs. We briefly review the stellar dynamics underlying this model in Section 3.1, and discuss the implementation of this theory in our simulations in Section 3.2.

A summary of our approach is as follows. For a given progenitor binary star formation history, IMF, and Hills mechanism rate, we predicted the number of observed main-sequence and evolved S-stars, the semi-major axes of their orbits around Sgr A^*, and their luminosity function. Additionally, we predicted the number and properties of the related HVSs observable in the Verberne et al. (2024, hereafter V24) survey. Fig. 2 shows a schematic overview of our modelling procedure. Finally, these predictions are simultaneously compared to (i) the observed properties of the S-star cluster, and (ii) the number of observed HVSs in the V24 survey.

3.1 Hills mechanism

In the Hills mechanism (Hills 1988), a stellar binary approaches a massive black hole within the tidal radius, where the tidal force of the black hole separates the binary, ejecting one star as an HVS and capturing the other star in a tight, elliptical orbit (for a review, see Brown 2015). The captured stars have been suggested as the possible origin of the S-stars (e.g. Gould & Quillen 2003; Ginsburg & Loeb 2006).

Progenitor binaries are expected to be on approximately parabolic orbits (Kobayashi et al. 2012), for which the ejection probability is equal for both stars and does not depend on the mass ratio (Sari et al. 2009). In this three-body encounter, there is an energy exchange, where this amount $$\mathrm{\Delta }E={\alpha }^2\frac{m_{\mathrm{e}\mathrm{j}}{m}_{\text{cap}}G}{a_{\mathrm{b}}}{\left(\frac{M_{\mathrm{b}\mathrm{h}}}{m_{\mathrm{b}}}\right)}^{1/3}\text{and}{m}_{\mathrm{b}}=m_{\text{cap}}+m_{\mathrm{e}\mathrm{j}},$$(1)

is gained by the ejected star at the expense of the captured one. In the expression for ΔE, α is a prefactor of order unity that depends on the geometry of the encounter (the inclination of the binary orbital plane, the binary phase at closest approach, the pericentre radius, and the binary eccentricity), m_ej and m_cap are the mass of the ejected and captured star respectively, G the gravitational constant, and a_b the semi-major axis of the progenitor binary. Thus the ejection velocity of the ejected star is $$V_{\mathrm{e}\mathrm{j}}=\sqrt{\frac{2\mathrm{\Delta }E}{m_{\mathrm{e}\mathrm{j}}}}=\alpha \sqrt{\frac{2G{m}_{\text{cap}}}{a_{\mathrm{b}}}}{\left(\frac{M_{\mathrm{b}\mathrm{h}}}{m_{\mathrm{b}}}\right)}^{1/6}$$(2)

(e.g. Rossi et al. 2014). The captured star, on the other hand, will travel on an orbit with a semi-major axis equal to $$a_{\text{cap}}=\frac{G{M}_{\mathrm{b}\mathrm{h}}{m}_{\text{cap}}}{2\mathrm{\Delta }E}=\frac{a_{\mathrm{b}}{M}_{\mathrm{b}\mathrm{h}}}{2{\alpha }^2{m}_{\mathrm{e}\mathrm{j}}}{\left(\frac{m_{\mathrm{b}}}{M_{\mathrm{b}\mathrm{h}}}\right)}^{1/3}$$(3)

after the separation of the binary system. In Fig. 3, we show the semi-major axis of the captured star as a function of the total binary mass and the mass ratio of the captured over the ejected star, where we take α = 1. This figure clearly demonstrates that Hills mechanism disruptions of binaries will result in stars captured into orbits with semi-major axes similar to those of the observed S-stars.

The initial pericentre distance (r_p) of the captured star is equal to the pericentre distance of the binary around the MBH and is set by the penetration factor (β) and tidal radius (r_t) defined as $$\beta =\frac{r_{\mathrm{t}}}{r_{\mathrm{p}}}\text{and}{r}_{\mathrm{t}}=a_{\mathrm{b}}{\left(\frac{M_{\mathrm{b}\mathrm{h}}}{m_{\mathrm{b}}}\right)}^{1/3}.$$(4)

We can thus write the initial eccentricity of the captured star as $$e_{\text{init}}=1-\frac{r_{\mathrm{t}}}{\beta a_{\text{cap}}}=1-\frac{a_{\mathrm{b}}}{\beta a_{\text{cap}}}{\left(\frac{M_{\mathrm{b}\mathrm{h}}}{m_{\mathrm{b}}}\right)}^{1/3}.$$(5)

We refer to this pericentre distance and eccentricity as initial, because the angular momentum of a star deposited in the S-star cluster is expected to change rapidly due to dynamical resonant relaxation processes (e.g. Rauch & Tremaine 1996; Hopman & Alexander 2006; Binney & Tremaine 2008; Antonini & Merritt 2013; Generozov & Madigan 2020). The timescale on which the energy of the orbit changes is much longer and happens on the two-body relaxation timescale. The precise timescales on which the angular momentum and energy change are still an active field of research (see review by Alexander 2017). Additionally, stellar collisions might play an important role near Sgr A* (e.g. Sari & Fragione 2019; Rose et al. 2023; Rose & Mockler 2024; Ashkenazy & Balberg 2025).

Fig. 2 Schematic overview of the steps of our forward model that predicts the numbers of observed HVSs, main-sequence S-stars, and evolved S-stars.

3.2 Simulations

Having discussed the theory behind the Hills mechanism and its potential to explain the S-star cluster, we now describe our simulations.

In our model, progenitor binaries are characterised by mass, age, metallicity, mass ratio, and binary separation. We assumed stars with solar metallicity (Z_⊙ = 0.0142; Asplund et al. 2009) for simplicity and evaluated the impact of this assumption in Section 5.3. Other properties of our progenitor binaries were drawn from probability distributions. In particular, the distribution of the mass ratio is f(q) ∝ q^γ, where 0.1 ≤ q ≤ 1 and we used an orbital period (P) distribution of the form log(P/[1 sec]) ∝(logP/P₀)^π, where P₀ is 1 sec. Our prior ranges are −2 ≤ γ ≤ 2 and 0 ≤ π ≤ 2. The mass and age of the binaries were determined given a star formation history and an initial mass function (IMF). Using the star formation history and the IMF, we can compute the ‘present-day’ mass function at any look-back time. For a Hills mechanism disruption at look-back time t, we calculated the number of stars alive per unit mass m and unit age A from $$\frac{\mathrm{d}N}{\mathrm{d}m\mathrm{d}A}(m,A\mid t)=\mathrm{SFR}(A+t)C(A)P(m\mid A),$$(6)

where SFR(A + t) is the star formation rate at time A + t, C(A) the fraction of stars alive at age A given an IMF, and P(m∣A) the mass function for stars alive of age A, which is the IMF truncated at the maximum mass of a star alive with age A. If a star survives till the present day, it has an age A + t. The fraction of stars alive at age A is calculated by $$C(A)=\int P(m)(m<m_{max}(A))\mathrm{d}m,$$(7)

where P(m) is the normalised IMF, and m_max(A) the maximum mass star alive at age A. We assumed the rate of binary disruptions at a point in time is proportional to the number of stars alive at that time. The Hills mechanism rate we report is always the current rate.

Unless stated otherwise, we assumed that the progenitor binaries are disrupted isotropically, which means that the resulting ejected and captured stars have random orientations. We calculated the energy that is exchanged in the tidal separation from Equation (1). We used the results from the simulations in Sersante et al. (in prep.) and randomly sampled over the prefactor α, using the distribution shown in Fig. 4. We assumed here that the progenitor binaries are on circular orbits in the full lose-cone regime (Merritt 2013) such that β (Equation (4)) is distributed log-uniformly. We also investigated a thermal eccentricity distribution in Section 5.5.

Fig. 3 The semi-major axis of the captured star’s orbit around Sgr A* as a function of the progenitor binary’s total mass and mass ratio, for three different binary separations; the mass ratio is defined as the mass of the captured star over that of ejected star. We only show solutions for which the progenitor binary separation is larger than the sum of the stellar radii at zero-age main-sequence.

Fig. 4 Prefactor α that we sample over to account for random orbital phases and inclinations at the time of disruption (Sersante et al., in prep.). The black and red lines show the distribution for circular and thermally distributed progenitor eccentricities respectively.

3.2.1 S-star predictions

After the binary is disrupted, we calculated if the captured star is within the stellar tidal disruption radius of Sgr A^* (Rees 1988). If so, the star was removed from the simulation, if not we assumed that the eccentricity is rapidly thermalised such that n(e) ∝ e while the orbital energy remains constant. The reason for this is that simulations point out that resonant relaxation can thermalise the S-star eccentricity distribution over ∼10⁷ yr (e.g. Perets et al. 2009; Antonini & Merritt 2013; Generozov & Madigan 2020), which is of the order of the age of the CWD. The extended mass required for resonant relaxation to be fast enough is similar to the 1σ upper limit on the enclosed mass within the S2 star orbit (Generozov & Madigan 2020; GRAVITY Collaboration 2024). After the orbits were thermalised, we again removed any stars whose pericentre lies within the stellar tidal disruption radius of Sgr A^*. On the other hand, we kept the semi-major axes of stars fixed, because we assumed negligible changes in the orbital energy over the timescales of interest (see Section 5.4). We followed the stellar evolution of the captured stars using MIST until the present epoch and determined the Ks magnitude (Dotter 2016; Choi et al. 2016). We assumed a uniform foreground extinction of A_Ks = 2.5 (Schödel et al. 2010; Fritz et al. 2011) towards the GC and corrected the photometry accordingly.

We only investigated stars with Ks < 16, because of observational completeness (see Section 2.2). We assigned to each remaining star a probability of having a known orbital solution following Burkert et al. (2024). We used the equivalent evolutionary point (EEP) number in MIST to separate main-sequence stars from evolved stars at the terminal-age mainsequence (EEP = 454). This leaves us with a number of predicted main-sequence and evolved S-stars that we can compare to observations in Gillessen et al. (2017).

3.2.2 HVS predictions

Simultaneously to the analysis of the captured stars presented above, we tracked the ejected stars using Speedystar¹ (Rossi et al. 2017; Marchetti et al. 2018; Contigiani et al. 2019; Evans et al. 2022a,b). This code simulates the ejection, propagation, and evolution of HVSs and provides synthetic photometry and Gaia observables for simulated HVSs. Our procedure is identical to that presented in V24, Section 4.3, except that we now considered realistic star formation histories as described in Section 4. This was incorporated through the prescription in Equations (6) and (7).

V24 performed a dedicated survey of HVSs. They rely on the near-radial trajectories of HVSs to identify candidates and provide follow-up observations. The survey has a clearly defined selection function, allowing us to assess the detectability of any given simulated star. V24 identify a single previously known HVS, S5-HVS1 (Koposov et al. 2020): a 1700 km s⁻¹ star at about 9 kpc from us that originates in the GC, where it was ejected about 4.8 Myr ago. This is the only star unambiguously linked to a GC ejection and provides strong indirect evidence of the Hills mechanism.

For a given star formation history, we forward modelled the expected number of HVSs in the dedicated survey presented in V24. The forward modelling leaves us with a predicted (simulated) population of HVSs that we can compare to the one HVS found observationally in V24. Therefore, when we mention the predicted number of HVSs, we are referring to the predicted number of HVSs in the V24 survey. The main constraining power of these observations is in providing an upper limit to the Hills mechanism rate, since the detection of only one HVS is very informative.

We consider HVSs to be the best comparison point since, following their definition in V24, they are almost certainly formed through the Hills mechanism. Although different formation mechanisms might contribute to the observed S-star population, we know that Hills mechanism disruptions will put stars into S-star-like orbits (Section 3.1).

4 Progenitor populations

In this section we will investigate two specific progenitor populations for binaries disrupted through the Hills mechanism.

4.1 First scenario: binaries originated from the NSC

We start by considering if binaries originating from parsec scales in the NSC (half light-radius around ∼5–10 pc) can explain the observed properties of the S-star cluster and HVSs.

The origin of progenitor binaries on which the Hills mechanism operates is uncertain. If stellar relaxation dominates the influx of disrupted binaries, the binaries will tend to originate at or beyond the sphere of influence (∼3 pc) of the MBH (Perets et al. 2007). Additionally, Penoyre et al. (in prep.) find that the axi-symmetry of the central potential of the Milky Way produces a collisionless flux of binaries from around the sphere of influence radius that dominates over the collisional rate. Thus, it appears unavoidable that stellar binaries from the NSC are set on orbits that result in their tidal separation, although the rate is still uncertain. We therefore assumed here that all binaries originate from the NSC and have the same star formation history as the NSC.

The NSC stellar population is consistent with having a canonical IMF (Kroupa 2001). On the other hand, its star formation history is still debated, but it is known that it has not been constant (e.g. Pfuhl et al. 2011; Schödel et al. 2020; Chen et al. 2023), with most of the population having formed billions of years ago. The model of the star formation history we used is the super-solar metallicity (2Z_⊙) model excluding the very young stars in the CWD from Schödel et al. (2020, Table 3). In Fig. 5, we plot the fraction of the total stellar mass formed as a function of look-back time. We can see that the NSC is dominated by old stars, with only a few percent having been formed in the last billion years. We assumed that the star formation rate was constant within the intervals reported in Schödel et al. (2020).

We only considered captured stars with M > 0.8 M_⊙, since lower mass stars are not bright enough to be observed near the GC at Ks < 16. Evolved 1 M_⊙ stars, on the other hand, might still be observed at those magnitudes according to MIST. The sample is dominated by relatively low-mass binaries, since the bulk of the population is very old.

In the left panel of Fig. 6, we show our results for binaries originating from the NSC. Since the predicted numbers linearly scale with the Hills mechanism rate, their ratio remains constant: for every observed HVS we expect 0.15 and 3.7 × 10² main-sequence and evolved S-stars respectively within our brightness limit. This is incompatible with observations because

• for every observable main-sequence S-star, around 2.5 × 10³ evolved S-stars should be visible unless they are somehow depleted (see Section 5.2),

• to account for the number of observed main-sequence S-stars, the Hills mechanism would deplete the NSC in ∼100 Myr (M_NSC ∼ 2 × 10⁷ M_⊙; Feldmeier-Krause et al. 2017),

• and there should be ∼10² HVSs in the V24 survey footprint.

Fig. 5 Look-back time against the fraction of stellar mass formed over each star formation epoch, which encapsulates the star formation history of the nuclear star cluster as determined in Schödel et al. (2020) for a 2Z⊙ model. The bars indicate the bin sizes of Schödel et al. (2020) within which we assume the star formation rate was constant.

4.2 Second scenario: binaries originated from the CWD

We showed in the previous section that Hills mechanism disruptions of binaries originating from the overall population of the NSC cannot simultaneously explain the observed properties of the S-star cluster and HVSs. Here we will investigate an alternative progenitor binary population, originating from the CWD (see also Madigan et al. 2009, 2014; Koposov et al. 2020; Generozov & Madigan 2020; Generozov 2021). This alternative scenario is largely inspired by the overlap between the orbital plane of the progenitor binary of S5-HVS1 around Sgr A* and that of the CWD (see Fig. 10 in Koposov et al. 2020), which is suggestive since HVSs are ejected in the plane of the progenitor binary around the MBH. The age of the CWD is uncertain, but likely less than 10 Myr (Lu et al. 2013).

We assumed a single star formation episode 10 Myr ago, and we modeled it with a Gaussian distribution centred at that epoch, with a standard deviation of 1 Myr, and additionally a top-heavy IMF with a slope of −1.7 (Lu et al. 2013; Gallego-Cano et al. 2024). Only stars more massive than about 7.5 M_⊙ are bright enough to pass our brightness limit given the age of the cluster. A clear difference with respect to the NSC scenario is the relative abundance of massive binaries, due to the younger ages of the stars and their top-heavy IMF.

Because HVSs are ejected in the plane of the progenitor binary around the MBH, we included this selection effect in our analysis of the expected number of HVSs in the V24 survey. This was done by only considering HVS ejections in the plane of the CWD, described by an inclination and longitude of the ascending node of (i, Ω) = (130^∘, 96^∘), and a half-width at half-maximum of about 15^∘ (Yelda et al. 2014). For a description of this coordinate system, see Lu et al. (2009).

In the case of the NSC, we assumed that the progenitor binaries were on parabolic orbits around the MBH. Because the CWD is much closer to Sgr A^*, the centre of mass of the progenitor binary might approach SgrA^* on a trajectory with significantly lower total (negative) energy (i.e. lower eccentricity). To account for this, we sampled from the radial surface density of the CWD (∼r⁻²; e.g. von Fellenberg et al. 2022) between 1 and 8 arcseconds from Sgr A*. The sampled semi-major axis sets the initial energy of the binary. We then calculated the semi-major axis of the captured star around Sgr A* taking into account this initial energy. In practice, this has a limited impact, since the initial energy is much smaller than the exchanged energy by a factor of ∼10⁻⁵ for observable S-stars in our simulations.

In the right panel of Fig. 6, we show the Hills mechanism rate against the predicted number of observed HVSs, main-sequence S-stars, and evolved S-stars. Given this star formation history, we expect, within our brightness limit, 7.2 and 0.05 main-sequence and evolved S-stars, respectively for every HVS in the V24 survey. This means that the number of observed HVSs in V24 is consistent with the twelve observed main-sequence S-stars within our brightness limit. If S5-HVS1 was ejected from a binary originating in the CWD, it is thus likely that most if not all main-sequence S-stars in the S-star cluster were deposited through the Hills mechanism. On the other hand, this model has trouble explaining the evolved S-stars, since for every evolved S-star it predicts around 140 main-sequence S-stars. This further demonstrates that if we assume the S-star cluster is exclusively formed through the Hills mechanism, the ratio of evolved to main-sequence S-stars is sensitive to the star formation history of the progenitor binary population. However, there is no obvious single progenitor population of binaries that can explain the number of HVSs, main-sequence S-stars, and evolved S-stars simultaneously.

Fig. 6 Hills mechanism rate against the predicted number of observed HVSs, main-sequence S-stars, and evolved S-stars. The dashed horizontal lines correspond to the observed number in each of these categories. The black circles highlight where the model and data agree. Left: the model predictions when progenitor binaries for the Hills mechanism originate from the NSC. Right: same as the left panel, but for binaries originating from the young CWD.

4.3 Combined CWD and NSC progenitors

We argued in Section 4.1, that progenitor binaries exclusively originating from the NSC cannot explain the observed properties of HVSs and main-sequence/evolved S-stars. However, we would like to stress that theoretically, one would expect there to be a finite disruption rate >0 yr⁻¹ of binaries originating from the NSC. Furthermore, a young progenitor population is needed to produce sufficient main-sequence S-stars. A natural interpretation of the observed properties of the S-star cluster compared to our simulations is that the main-sequence S-stars were recently deposited on their current orbits by disruptions from a young population of stars, such as the CWD, while the evolved S-stars were deposited through disruptions of mainly old, low-mass binaries originating in, for example, the NSC. Here we investigate if this scenario is able to explain all the observables simultaneously.

If progenitor binaries are disrupted from both the NSC and CWD, we expect there to be two rates that make up the total Hills mechanism rate: a constant rate caused by the inflow of binaries from the NSC (Γ_NSC) and a transient rate from the CWD (Γ_CWD). Our model is a linear combination between these two rates. The three observables we fitted to are the number of HVSs, main-sequence S-stars, and evolved S-stars that are observed. We show the resulting log-likelihood of Γ_CWD and Γ_NSC in Fig. 7 in combination with the confidence regions of the 68th and 95th percentile. We find that the rate of the Hills mechanism for binaries originating from the NSC is ${1.4}_{-0.6}^{+0.9}\times {10}^{-5}{\mathrm{y}\mathrm{r}}^{-1}$ and that from the CWD is ${1.0}_{-0.2}^{+0.3}\times {10}^{-4}{\mathrm{y}\mathrm{r}}^{-1}$. The current rate of Hills mechanism disruptions from the CWD would thus have to be about an order of magnitude higher than the constant background rate of binary disruptions from the NSC in order to fit the observables. Over the past ∼10 Myr, we therefore predict that ∼90% of ejected stars are young (≲10 Myr). Our derived rates are comparable to old full and empty loss-cone rate estimates (Hills 1988; Yu & Tremaine 2003). In this scenario, S5-HVS1 is about 150 times more likely to have been ejected from the young progenitor population; this, combined with the age estimate for S5-HVS1 of less than 10⁸ yr (Koposov et al. 2020) and the alignment of the CWD plane with that of its progenitor binary, leads us to argue that S5-HVS1 came from the CWD and shares the same age.

Now that we have a combined CWD and NSC model, we evaluate if this model can reproduce all the observables of both the S-star cluster and HVSs in Fig. 8. Firstly, we recount that neither progenitor population, when considered individually, was able to simultaneously predict the observed numbers of HVSs, main-sequence S-stars, and evolved S-stars. However, as shown in the top panel of Fig. 8, the combined model successfully accounts for all three of these quantities at once. Moreover, both the K s-band and semi-major axis distributions are consistent between our model and the observations, with KS-test p-values of 0.88 and 0.55 respectively.

Fig. 7 Log-likelihood on the Hills mechanism rate for our combined CWD and NSC model. In orange we show the 68th and 95th percentile confidence regions of our fit.

Stellar tidal disruption events

We find that for both our progenitor populations about 20% of Hills mechanism disruptions result in a stellar tidal disruption event (TDE). The resulting TDE rate of the order 10⁻⁵ to 10⁻⁶ yr⁻¹ is consistent with extragalactic TDE rates (Stone & Metzger 2016). However, this estimate is strongly influenced by our assumption of a full loss-cone. In the empty loss-cone regime, the Hills mechanism contribution to the TDE rate would be about two orders of magnitude lower. Our estimated Hills mechanism contribution to the TDE rate of 10⁻⁵ yr⁻¹ should thus be seen as an upper limit.

5 Discussion

In this work we have investigated whether the Hills mechanism can explain both the observed properties of the S-star cluster and HVSs, and we found a plausible formation scenario that is consistent with data (see Section 4.3 and Fig. 8). In this Section we critically review our model assumptions in the first five subsections. We investigate further the comparison between the semi-major axis distribution of our simulated population and the observed one in Section 5.6. We discuss future prospects in Section 5.7.

5.1 Star formation history

In this work, we have investigated two progenitor binary populations. The most important parameter that sets the ratio of main-sequence to observed S-stars is the star formation history. To demonstrate this, we calculated the amount of time a star in the GC can be above our brightness limit as a function of its initial mass, and we show this in Fig. 9. First of all, we notice that only main-sequence stars more massive than about 4 M_⊙ can be observed out to the GC, which have a main-sequence lifetime of ≲ 170 Myr. Secondly, we notice that evolved, low-mass stars can be brighter than our threshold for about an order of magnitude longer than massive main-sequence stars. For an old population of stars, such as the NSC, the result will thus be that most of the stars above the brightness limit will be evolved, low-mass stars. For a young population, such as the CWD, most of the stars above the brightness limit will instead be massive main-sequence stars. The ratio of evolved to main-sequence stars therefore traces the star formation history.

Our results show that a young progenitor population is needed to produce sufficient main-sequence S-stars to match observations. Our conclusion here is in line with Generozov et al. (2025), who, through a different argument, recently found the need for a young population of stars to explain the main-sequence S-stars. However, we argue that binary disruptions exclusively from a young population such as that of the CWD cannot explain the presence of the evolved S-stars.

Our conclusions most importantly relate to the star formation history of the progenitor populations. Alternative populations of possible old progenitor binaries, such as the nuclear stellar disc (Schödel et al. 2023), result in similar conclusions. This is particularly relevant, since massive perturbers might increase the disruption rate of binaries originating at these radii by several orders of magnitude (Perets et al. 2007). In the case of the young progenitor binary population, we consider the CWD to be the most likely origin, given its close proximity to Sgr A^* and its young age.

As a more agnostic star formation history, we also investigated a constant star formation rate over the past 13 Gyr and a canonical IMF slope of −2.3 for our progenitor population. Based on these simulations, we expect there to be 0.11 and 1.5 × 10² main-sequence and evolved S-stars respectively within the brightness limit for every HVS in the V24 survey. This is qualitatively similar to the results for the NSC progenitor origin and also incompatible with observations.

In the context of star formation history, we would also like to point out that the commonly used single power-law slope for the (present-day) mass function implicitly assumes a constant star formation rate with an IMF that is not explicitly defined. Since the stellar lifetime roughly scales with M^−2.5 (e.g. Ryan & Norton 2010), the IMF can be approximated by subtracting about −2.5 from the power-law slope of the current mass function if star formation was constant over the last several billion years. For a population of stars with a current mass function slope of −2.3, the IMF would thus be about $\frac{\mathrm{d}N}{\mathrm{d}M}\propto M^{0.2}$, unless star formation was not constant. Such a highly top-heavy IMF is not realistic even for the top-heavy CWD (e.g. Bartko et al. 2010; Lu et al. 2013; Gallego-Cano et al. 2024).

Fig. 8 Comparison of our best fit model, including Hills mechanism disruptions from both the NSC and CWD, with observations. Top: number of stars from our simulation compared to the observations for HVSs, evolved S-stars, and main-sequence S-stars. The error bars show the Poisson error on the predicted numbers from the model. Bottom left: cumulative distribution function of the visual Ks magnitude of S-stars for our simulations compared to the observations. Bottom right: cumulative distribution function of the semi-major axis of S-stars for our simulations compared to the observations.

Fig. 9 Initial stellar mass against the time spent with a magnitude of less than 16 in the Ks band, when the star is located at the GC. We utilise the stellar evolution code MIST.

5.2 Stellar collisions

So far, we have ignored the effects of stellar collisions. In general, in a dense cluster of stars and black holes with a high velocity dispersion, like that occupied by the S-stars, stellar collisions might become important. There is currently no consensus as to the importance of collisions for shaping the S-star cluster. Collisions may shape the density profiles of NSCs on a scale similar to that of the S-star cluster, depleting stars at small radii and potentially even leading to mergers (Rose et al. 2023; Rose & Mockler 2024; Ashkenazy & Balberg 2025). Generozov et al. (2025), however, find that collisions only have a minor impact on their results.

It should also be noted that there is no consensus as to the net effect of collisions. Depending on the mass ratio, evolutionary stage of each component involved, impact parameter, and impact velocity the result of a collision will vary drastically (e.g. Benz & Hills 1987; Trac et al. 2007; Mastrobuono-Battisti et al. 2021; Rose et al. 2023; Rose & Mockler 2024).

In Fig. 10, we show the collisional timescale as a function of eccentricity and semi-major axis for a Bahcall & Wolf (1976) cusp following Sari & Fragione (2019). We can see that stellar collisions might be important for a significant part of the parameter space. Particularly for low-mass stars that could otherwise reside in the S-star cluster for several billion years during their main-sequence, collisions might play an important role. The massive main-sequence S-stars will not be significantly affected, since their main-sequence lifetimes tend to be shorter than the collisional timescale shown here.

Collisions are thus relevant for our NSC progenitor population model. To ensure our conclusions are robust, we also calculated our NSC progenitor results for ejections only over the past 100 Myr. The effect of only considering recent ejections is that the predicted number of evolved S-stars decreases. However, even if we only consider these recent binary disruptions, we would still expect about 24 evolved S-stars for every main-sequence S-star. This means that even when only considering recent disruptions, an old progenitor binary population is inconsistent with observations. Because possible depletion will remove a fraction of the simulated S-stars, we expect that the Hills mechanism rate constraint for the NSC in our combined model is a lower limit. The true rate would need to be higher to account for the depletion rate over time.

Fig. 10 Collisional timescale of a 10 M⊙ star with a 5.6 R⊙ radius as a function of stellar eccentricity and semi-major axis. The stars with orbital solutions (see Fig. 1) are overplotted.

5.3 Metallicity

The metallicity of the S-stars is, as of yet, unknown (Habibi et al. 2017). For this reason, it is unclear what metallicity to assume for the progenitor binary population.

In the case of the NSC, most studies find a super-solar mean metallicity of [Fe/H] ∼ 0.2 to 0.3 (e.g. Schultheis et al. 2019; Schödel et al. 2020) with a minor population of a few percent being metal poor (e.g. Do et al. 2015; Feldmeier-Krause et al. 2017; Gallego-Cano et al. 2024). The metallicity of the CWD, like that of the S-stars, is unknown, though in situ star formation in the GC is generally expected to be (super) solar (e.g. Do et al. 2015; Feldmeier-Krause et al. 2017; Schödel et al. 2020).

Fortunately, the effect of the assumed metallicity in this study is minimal, since we only consider if stars are on the mainsequence or evolved, in combination with a brightness limit. If we assume a metallicity of [Fe/H] = 0.3, we find that the number of main-sequence S-stars to evolved stars changes from 7.0 to 7.5 and the number of evolved S-stars from 0.05 to 0.04 for the CWD binary origin. Also for the HVS observations, metallicity is unimportant, because of the large coverage of colour-magnitude space (V24).

5.4 Energy diffusion

In this paper, we assumed that the energy of any captured stars does not evolve after the tidal interaction. In reality, the energy diffuses on the two-body relaxation timescale (e.g. Bar-Or & Alexander 2016; Sari & Fragione 2019). If we use the formalism from Sari & Fragione (2019), while assuming a Bahcall & Wolf (1976) cusp of stars, we get the two-body relaxation timescale in energy shown in Fig. 11. We can see that only S2 has a shorter energy relaxation time than 1 Gyr, which is still much longer than its main-sequence lifetime (Habibi et al. 2017). We noted in Section 5.1, that only main-sequence stars with masses above ∼4 M_⊙ are expected to pass our brightness limit, for which the main-sequence lifetime is about 170 Myr. This is much shorter than the typical two-body relaxation timescale for the energy shown in Fig. 11. We conclude that for the assumed density profile, our assumption that the captured stars’ energy is constant is safe for main-sequence stars. If, instead, a star was deposited >1 Gyr ago, the energy might have changed significantly. This is the case for very old disruptions from the NSC progenitor population discussed in Section 4.1.

Fig. 11 Two-body energy relaxation timescale as a function of eccentricity and semi-major axis assuming a Bahcall & Wolf (1976) cusp. The stars with orbital solutions are overplotted (see Fig. 1). The two solid lines mark the 1 and 10 Gyr relaxation times.

5.5 Progenitor binary eccentricity

The internal eccentricity of progenitor binaries influences the energy that is exchanged in the Hills mechanism through the prefactor that we plot in Fig. 4. If we assume a thermal eccentricity distribution for our progenitor binaries originating in the CWD, the number of main-sequence and evolved S-stars per observed HVS change from 7.2 and 0.05 to 7.1 and 0.06 respectively. Moreover, if we compare the cumulative distribution functions of the simulated semi-major axes of the S-stars, we find a maximum difference of about 0.02, which is well below statistical uncertainties given the number of observed S-stars. Assuming a thermal eccentricity distribution or a circular one makes thus little difference for our results. This is explained by considering that disruptions that result in observable S-stars tend to be of high energy, which means high values for the prefactor. For high prefactors, we can see from Fig. 4 that the thermal eccentricity distribution is very similar to the circular one.

5.6 The semi-major axis distribution

So far, we have not discussed the semi-major axis distribution of observed and modelled S-stars in-depth. The reason is that it is the most difficult observable to interpret. One needs to model the impact of collisions in combination with energy and angular momentum diffusion and observational biases, before attempting to compare a simulated semi-major axis distribution with an observed one.

We showed in Fig. 8 how our model compares to the observed semi-major axis distribution. We found the two to be consistent, but the relatively limited number of S-stars for which we can assume the observations to be complete limits our sensitivity. The lack of predicted S-stars at semi-major axes greater than about 2 × 10⁴ AU is not surprising, since this is caused by the observational selection function (Burkert et al. 2024). However, our model is able to naturally explain the lack of S-stars at semi-major axes below about 10³ AU, without the need for stellar collisions. The limiting factors for producing stars at such radii are the energy that can be exchanged in the Hills mechanism and the high initial eccentricity that can cause TDEs to occur for small semi-major axes.

In the literature, a ‘zone of avoidance’ was identified in pericentre distance against eccentricity where no S-stars have been found, even though observations should have identified them if they exist. Burkert et al. (2024) defines the region by log(r_p) < 1.57 + 2.6 ×(1−e), where r_p and e are the pericentre distance and eccentricity respectively. In Fig. 12, we show the zone of avoidance in semi-major axis and eccentricity space (which are independent quantities), in combination with the predicted S-star cluster distribution of our CWD origin simulations. The figure demonstrates that our simulations can naturally explain the lack of observed stars in the zone of avoidance. The fundamental reasons being the limit on the energy that can be exchanged in the Hills mechanism (Equation (1)), captured stars ending up as TDEs at small semi-major axis (since their initial eccentricity is close to 1), and subsequent relaxation of the eccentricity distribution. Our results therefore agree with the conclusion of Generozov et al. (2025), that the Hills mechanism can naturally explain the lack of orbital solutions in the zone of avoidance.

Fig. 12 Eccentricity against semi-major axis for our simulation of the CWD origin binaries. The colour bar shows the number of S-stars predicted for an ejection rate of 10−2 yr−1. In white we show the boundary of the zone of avoidance identified in Burkert et al. (2024). The drop in the number of simulated S-stars at large semi-major axes is due to the selection function from Burkert et al. (2024).

5.7 Prospects

Definitively identifying the formation origin of the S-star cluster is very challenging. Key observables of the S-stars that could assist in this are in particular the individual ages, metallicities, and masses of the S-stars. These quantities could help distinguish if the S-stars were formed from multiple populations of stars, or a single star formation episode, which is especially useful for evaluating scenarios. Of particular interest for the Hills mechanism origin hypothesis is the age and metallicity distribution of the evolved S-stars. If the evolved S-stars are old, low-mass stars it would suggest that the main-sequence stars have a different origin. Either the Hills mechanism efficiently disrupts binaries from young structures such as the CWD or alternative scenarios are needed to explain their presence near Sgr A^*, such as those explored by Levin (2007), Rantala & Naab (2024), and Akiba et al. (2024).

A deeper completeness limit of S-stars observed near Sgr A* would also help eliminate formation scenarios, since it would allow for more detailed comparisons using the Ks-band luminosity function and semi-major axis distributions in particular. This deeper completeness limit will be allowed by new instruments such as GRAVITY(+) (Eisenhauer et al. 2011; GRAVITY+ Collaboration 2022), and an increased temporal baseline for the observations.

Specifically for relating the Hills mechanism to the S-star cluster, a larger sample of confirmed HVSs originating in the GC is vital. This would allow for much lower uncertainties in the Hills mechanism rate, even for a single new HVS discovery. Additionally, measurements of the metallicities, ages, and masses of HVSs would provide strong constraints on their progenitor population. Large spectroscopic surveys such as Gaia DR4, DESI (Cooper et al. 2023), WEAVE (Dalton et al. 2014), and 4MOST (de Jong et al. 2019) will provide unprecedented numbers and depth of observations and are expected to discover new HVSs (V24).

6 Conclusions

In this work, we investigated if binary disruptions through the Hills mechanism can explain both the observed population of S-stars and state-of-the-art HVS observations. We forward modelled the S-stars and HVSs and compared both populations to the latest observational results, taking observational selection effects into account. We find that, in particular, the ratio of main-sequence to evolved S-stars is highly informative on the star formation history of the progenitor binary population and investigated different progenitor binary populations.

We show that an old progenitor binary population, where the majority of the stellar content formed billions of years ago such as the NSC, cannot explain the ratio of main-sequence to evolved S-stars, nor the relative number of HVSs. Furthermore, our results show that a young progenitor binary population, such as the CWD, can simultaneously explain the observed population of main-sequence S-stars and HVS observations, but cannot explain the evolved S-stars that are observed.

We argue that the main-sequence and evolved S-star populations, if formed through the Hills mechanism, have different progenitor populations. Only the model in which progenitor binaries come from both the CWD and NSC can successfully the numbers of HVSs, main-sequence S-stars, and evolved Sstars. Furthermore, we show that the Ks-band and semi-major axis distributions predicted by this model are consistent with observations. The current Hills mechanism rate from the young binary population would need to be about an order of magnitude higher than that of the old population at roughly 10⁻⁴ and 10⁻⁵ yr⁻¹ respectively. About 90% of ejected stars over the past ≲ 10 Myr should thus originate from the CWD. In general, we expect most young HVSs to be ejected in bursts coinciding with star formation near Sgr A*.

We additionally argue that, provided S5-HVS1 was formed through the Hills mechanism, most (if not all) S-stars were formed through Hills mechanism disruptions of binaries.

Upcoming large spectroscopic surveys are expected to improve our understanding of the rate of the Hills mechanism, and thereby provide strong constraints on the Hills mechanism’s contribution to the formation of the S-star cluster. In addition, continued observations of the GC will allow more precise measurements of the properties of the S-stars, facilitating increasingly constraining tests of different formation scenarios.

Acknowledgements

The authors thank Stephan Gillessen, Antonia Drescher, Diogo Ribeiro, Matteo Bordoni, Bianca Sersante, and Yuri Levin for helpful discussions. E.M.R. acknowledges support from European Research Council (ERC) grant number: 101002511/project acronym: VEGA_P. S.K. acknowledges support from the Science & Technology Facilities Council (STFC) grant ST/Y001001/1. For the purpose of open access, the author has applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submission. Software: NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), Matplotlib (Hunter 2007), Astropy (Astropy Collaboration 2013, 2018, 2022), isochrones (Morton 2015), Speedystar (Contigiani et al. 2019; Evans et al. 2022a).

References

All Figures

Fig. 1 Orbital solutions for stars at the GC found in Gillessen et al. (2017), with the exception of S5, which is from von Fellenberg et al. (2022). Stars that are considered part of the young disc are coloured blue, which are S66, S67, S83, S87, S91, S96, S97, and R44. S111 is omitted, because of the hyperbolic solution. The grey points show S-stars fainter than our completeness limit, while black points show those that are brighter.

In the text

	Fig. 2 Schematic overview of the steps of our forward model that predicts the numbers of observed HVSs, main-sequence S-stars, and evolved S-stars.
In the text

Fig. 3 The semi-major axis of the captured star’s orbit around Sgr A* as a function of the progenitor binary’s total mass and mass ratio, for three different binary separations; the mass ratio is defined as the mass of the captured star over that of ejected star. We only show solutions for which the progenitor binary separation is larger than the sum of the stellar radii at zero-age main-sequence.

In the text

	Fig. 4 Prefactor α that we sample over to account for random orbital phases and inclinations at the time of disruption (Sersante et al., in prep.). The black and red lines show the distribution for circular and thermally distributed progenitor eccentricities respectively.
In the text

	Fig. 5 Look-back time against the fraction of stellar mass formed over each star formation epoch, which encapsulates the star formation history of the nuclear star cluster as determined in Schödel et al. (2020) for a 2Z⊙ model. The bars indicate the bin sizes of Schödel et al. (2020) within which we assume the star formation rate was constant.
In the text

Fig. 6 Hills mechanism rate against the predicted number of observed HVSs, main-sequence S-stars, and evolved S-stars. The dashed horizontal lines correspond to the observed number in each of these categories. The black circles highlight where the model and data agree. Left: the model predictions when progenitor binaries for the Hills mechanism originate from the NSC. Right: same as the left panel, but for binaries originating from the young CWD.

In the text

	Fig. 7 Log-likelihood on the Hills mechanism rate for our combined CWD and NSC model. In orange we show the 68th and 95th percentile confidence regions of our fit.
In the text

Fig. 8 Comparison of our best fit model, including Hills mechanism disruptions from both the NSC and CWD, with observations. Top: number of stars from our simulation compared to the observations for HVSs, evolved S-stars, and main-sequence S-stars. The error bars show the Poisson error on the predicted numbers from the model. Bottom left: cumulative distribution function of the visual Ks magnitude of S-stars for our simulations compared to the observations. Bottom right: cumulative distribution function of the semi-major axis of S-stars for our simulations compared to the observations.

In the text

	Fig. 9 Initial stellar mass against the time spent with a magnitude of less than 16 in the Ks band, when the star is located at the GC. We utilise the stellar evolution code MIST.
In the text

	Fig. 10 Collisional timescale of a 10 M⊙ star with a 5.6 R⊙ radius as a function of stellar eccentricity and semi-major axis. The stars with orbital solutions (see Fig. 1) are overplotted.
In the text

	Fig. 11 Two-body energy relaxation timescale as a function of eccentricity and semi-major axis assuming a Bahcall & Wolf (1976) cusp. The stars with orbital solutions are overplotted (see Fig. 1). The two solid lines mark the 1 and 10 Gyr relaxation times.
In the text

Fig. 12 Eccentricity against semi-major axis for our simulation of the CWD origin binaries. The colour bar shows the number of S-stars predicted for an ejection rate of 10−2 yr−1. In white we show the boundary of the zone of avoidance identified in Burkert et al. (2024). The drop in the number of simulated S-stars at large semi-major axes is due to the selection function from Burkert et al. (2024).

In the text