© The Authors 2025

1 Introduction

Galaxies with a stellar mass of ≲10¹⁰M_⊙¹ typically host a nuclear stellar cluster (NSC) in their most central region (Neumayer et al. 2020, and references therein). At the centre of these structures, there is typically a supermassive black hole (SMBH, Neumayer & Walcher 2012; Nguyen et al. 2019) with a mass between 10^6–9 M_⊙ (see e.g. Kormendy & Ho 2013). The Galactic centre (GC) of the Milky Way is no exception. The radio source Sagittarius A* (Sgr A*) is associated with a SMBH of 4 × 10⁶ M_⊙, which in turn is surrounded by a 2.5 × 10⁷ M_⊙ NSC (see e.g. Schödel et al. 2003, 2014; Eisenhauer et al. 2005; Ghez et al. 2005, 2008; Gillessen et al. 2009b, 2017; Alexander & Pfuhl 2014; Boehle et al. 2016; GRAVITY Collaboration 2018).

Very close to Sgr A* there exists a young population of stars that are usually referred to as the S-cluster, and this population indicates that recent star formation has taken place (e.g. Eisenhauer et al. 2005; Genzel et al. 2010; Schödel et al. 2020). So far, finding binaries among these stars has not proven very successful (Chu et al. 2023). Only a small number of main-sequence (MS) binaries has been detected in the GC², and those that have been identified are all massive (Martins et al. 2006; Pfuhl et al. 2014). In contrast, little is known about a less massive and fainter MS population in the NSC. Recently, Peißker et al. (2024) presented the detection of a lower-mass spectro-scopic binary (2.80 ± 0.50 M_⊙ and 0.73 ± 0.14 M_⊙ components) within the S-cluster. One idea, which is further motivated by this detection, is that these lower-mass binaries could be the progenitors of the G-objects, which appear as dust-enshrouded self-gravitating cloud-like stellar sources (Ghez et al. 2005; Gillessen et al. 2012; Ciurlo et al. 2020). All binary systems in the inner parts of the NSC form a hierarchical triple system with the SMBH. As such, the SMBH acts as a distant perturber that induces eccentricity and inclination oscillations in the binaries, the so-called von Zeipel-Lidov-Kozai (ZLK, von Zeipel 1910; Lidov 1962; Kozai 1962) mechanism. These oscillations can be extreme enough to drive the binaries to near parabolic orbits where they may merge because of the finite size of the constituent stars. This process is one of the leading candidate mechanisms for the formation of the G-objects (Prodan et al. 2015; Stephan et al. 2016, 2019; Ciurlo et al. 2020; see also Burkert et al. 2012; Schartmann et al. 2012; Ballone et al. 2013; Shcherbakov 2014; De Colle et al. 2014; Mapelli & Ripamonti 2015; Ballone et al. 2016; Trani et al. 2016; Calderón et al. 2016; Zajaček et al. 2017; Calderón et al. 2018; Ballone et al. 2018; Trani et al. 2019a). Another mechanism responsible for producing binary mergers is frequent encounters between binaries and surrounding stars (see Hamers & Samsing 2019; Samsing et al. 2019; Michaely & Perets 2019; Young & Hamers 2020, and this study).

While Prodan et al. (2015); Stephan et al. (2016, 2019) have all looked at the dynamical evolution of binaries via the ZLK mechanism, their simulations are confined to the inner 0.1 pc of the NSC because further out, the general relativity precession (GRP) timescale becomes short enough that ZLK may be suppressed (Ford et al. 2000; Liu et al. 2015; Naoz et al. 2013; Naoz 2016). Additionally, far out in the NSC, other relevant dynamical timescales become long enough that they are not too different from the Galactic field (see e.g. Rose et al. 2020). This leaves a small region (~0.1–0.3 pc) where the dynamical evolution of a binary is dominated by its interactions with nearby stars and where associated dynamical processes occur on a timescale shorter than or comparable to the age of the NSC (~10 Gyr).

In this paper, we aim to i) highlight the different outcome drivers that determine the fate of binaries and ii) relate different types of binary remnants to the variety of objects that exist in the NSC, with an extra focus on binary mergers. The rest of the paper is structured as follows: In Section 2 we describe our model for the NSC (Sect. 2.1) and the dynamical processes considered for this study (Sects. 2.2–2.4). In Section 3, we describe the binary population we consider (Sect. 3.1) and how we evolve our binaries (Sect. 3.2; see also Appendix A). Section 4 shows the results we obtained, which are divided into mergers (Sect. 4.1), evaporation (Sect. 4.2), and the remaining outcomes (Sect. 4.3).³ This is followed by a general discussion in Section 5, including a more focused discussion on binary mergers (Sect. 5.1) and comparisons with other relevant studies (Sect. 5.2). Finally, our conclusions are presented in Section 6.

2 Dynamics in the nuclear stellar cluster

The orbits of binary systems in the NSC are modified by various processes. Because all binaries form a hierarchical triple system with the SMBH at the centre of the NSC, both general relativity (GR) and ZLK oscillations affect their evolution (Prodan et al. 2015; Naoz 2016; Stephan et al. 2016, 2019). Moreover, as the stellar density in the NSC itself is very high, encounters with surrounding stars cannot be ignored and drive both evaporation (Heggie 1975a; Rose et al. 2020) and mergers (Hamers & Samsing 2019). In the following sections, we outline our approach to dealing with these mechanisms.

2.1 NSC model

The presence of an SMBH has a number of effects on the surrounding stars and the NSC. As the cluster ages and relaxes, it reaches a steady-state characterised by a cuspy power-law stellar density profile (see Eq. (1)) with an exponent α = 1.75, which is typically referred to as a Bahcall-Wolf profile, (Bahcall & Wolf 1976)⁴ $$\rho \left(r_{•}\right)={\rho }_0{\left(\frac{r_{•}}{r_0}\right)}^{-\alpha }.$$(1)

Here, r_• is the distance to the SMBH, and we adopted the scale density ρ₀ = 1.35 × 10⁶ M_⊙ pc⁻³ and r₀ = 0.25pc following Genzel et al. (2010). The gravitational field in the inner parsec of the NSC is dominated by the SMBH. Therefore, stars inside this radius tend to move on Keplerian orbits. Additionally, the resulting one-dimensional velocity dispersion, σ, related to the density distribution via the Jeans equation scales with the mass of the SMBH, M_•, $$\sigma \left(r_{•}\right)=\sqrt{\frac{G{M}_{•}}{r_{•}\left(1+\alpha \right)}}.$$(2)

Here, G is the gravitational constant. We adopted the Bahcall–Wolf profile of α = 1.75 in these simulations, but we note that observations of stars in the NSC suggest slightly shallower profiles with values between α ~ 1.25–1.5 (Gallego-Cano et al. 2018).

2.2 Dynamical three-body encounters

As a binary interacts with a third star (tertiary), there is an exchange of energy and angular momentum that changes the separation and eccentricity of the binary. The separation is inversely proportional to the binary’s internal binding energy, and as such, if the binary loses (gains) energy, its separation grows (shrinks). Typically, one refers to the growing and shrinking as softening and hardening, respectively, which is related to the binary’s hardness ratio, that is, the ratio of its binding energy to the average kinetic energy of the surrounding stars: $$h=\frac{2E_{\text{bin}}}{{\overline{m}}_3{\sigma }^2},$$(3)

where ${\overline{m}}_3$ and σ are the average mass and velocity dispersion of the surrounding stars, respectively. A binary is considered hard if h > 1 and soft if h < 1. Based on the seminal work of Heggie (1975a), it has been shown that there is a tendency for soft binaries to soften and for hard binaries to harden with time; a trend usually referred to as Heggie’s law. A soft binary is therefore expected to soften until it eventually evaporates. In the inner regions of the NSC, the velocity dispersion is on the order of a few hundred km s⁻¹ (see Eq. (2)), which is more than an order of magnitude higher than for a typical globular cluster (see e.g. Ivanova et al. 2005). For a typical population of binaries, such as those we find in the field (Raghavan et al. 2010; Moe & Di Stefano 2017), the majority would be considered soft in the NSC and would therefore be expected to evaporate (see e.g. Prodan et al. 2015; Stephan et al. 2016, 2019; Panamarev et al. 2019).

2.3 von Zeipel–Lidov–Kozai mechanism

As mentioned previously, all binaries in the NSC form hierarchical triple systems with the SMBH. The SMBH exerts small but consistent gravitational perturbations on the binary’s orbit, which can induce periodic oscillations between the binary’s eccentricity and inclination. It is possible to describe the long-term evolution of such a system with Hamiltonian perturbation theory. We base our implementation on the quadrupole approximation and the equations derived in Naoz et al. (2013), and refer the reader to Naoz (2016) for a review. More details on our implementation can also be found in Appendix B.

The ZLK mechanism operates under resonance; any perturbations of this resonance can work to suppress the ZLK oscillations. GRP is one such mechanism that is highly relevant in the NSC. If the timescale for GRP is shorter than the quadrupole timescale, ZLK is suppressed (e.g. Naoz et al. 2013; Stephan et al. 2019; Ford et al. 2000). The quadrupole timescale is given in Eq. (4), where m_bin is the binary mass, M_• is the mass of the SMBH, P_inner (P_•) is the orbital period for the inner (outer) orbit, and e_• is the orbital eccentricity of the outer orbit: $$t_{\text{quad}}\approx \frac{16}{30\pi }\frac{m_{bin}+M_{•}}{M_{•}}\frac{P_{•}^2}{P_{\text{inner}}}{\left(1-e_{•}^2\right)}^{3/2}.$$(4)

The ratio between this timescale and the GRP timescale can be estimated from Naoz (2016, Eq. (60)) and the dimension-less parameter, ϵ_GR, introduced in Liu et al. (2015). These two expressions are respectively given in Eqs. (5) and (6): $$\frac{t_{\text{GRP}}}{t_{\text{quad}}}={ϵ}_{\text{GR}}^{-1}\left(1-e_{\text{bin}}^2\right)$$(5) $${ϵ}_{\text{GR}}=\frac{3G{m}_{\text{bin}}^2{a}_{•}^3{\left(1-e_{•}^2\right)}^{3/2}}{a_{\text{bin}}^4{c}^2{M}_{•}}.$$(6)

The ratio t_GRP/t_quad depends on the separation and eccentricity of the inner orbit a_bin and e_bin, respectively, the binary mass m_bin, the separation and eccentricity of the outer orbit a_• and e_• respectively, and the mass of the SMBH M_•. G is the gravitational constant and c is the speed of light. ZLK oscillations are only possible if this ratio is greater than one, and as long as the binary-SMBH system is hierarchical: a_inner/a_outer < 0.1.

2.4 Tidal dissipation

The statistical tendency for binaries to thermalise as they interact with surrounding stars leads to more highly eccentric orbits (Heggie 1975a,b). Similarly, the ZLK mechanism also provides a pathway for producing highly eccentric orbits. When stars in an eccentric binary are at periapsis, tidal forces coupling to the internal structure of the stars (e.g. convective motions) leads to orbital angular momentum transfer and energy dissipation to the stars. As a consequence, the binary orbit circularises and the binary becomes harder; this may lead to the binaries either merging or remaining intact for a ~ Hubble time, as they become virtually stable against evaporation. We base the following on Mardling & Aarseth (2001) and refer the reader to this paper for more detailed explanations regarding tidal dissipation. Importantly, two distinct dissipation mechanisms are considered, as referred to by the original authors: normal (or equilibrium) and chaotic dissipation.

2.4.1 Equilibrium tidal evolution

Equilibrium tides arise from a quasi-hydrostatic deformation of the star (Zahn 1977). This causes a net torque, in which the orbital angular momentum and orbital energy is dissipated into the stars, leading to more circular orbits with smaller separations (Triaud et al. 2017). In the equilibrium model, the eccentricity changes as $$\begin{array}{l}\frac{de}{dt}=-\frac{2}{4{\pi }^2}\left(\frac{C_1{Q}_1^2}{w_1^4}q\left(1+q\right)+\frac{C_2{Q}_2^2}{w_2^4}\frac{1}{q}\left(1+\frac{1}{q}\right){\left(\frac{r_2}{r_1}\right)}^8\right)\\ \times \frac{f_e\left(e\right)}{{\left(p_0\left(1+e_0\right)\right)}^8},\end{array}$$(7)

where $$f\left(e\right)=\frac{9\pi }{5}e{\left(1-e^2\right)}^{3/2}\left(1+\frac{15}{4}{e}^2+\frac{15}{8}{e}^4+\frac{5}{64}{e}^6\right).$$(8)

We note that 𝒞_i = w_i/2πQ_i, Q = 10⁴ is some nominal value used for the damping timescales defined in Mardling & Aarseth (2001), w_i is the oscillation frequency, and Q_i are the associated overlap integrals (Press & Teukolsky 1977). These are all structure constants dependent on the polytrope of choice, which we take to be n = 3 as a reasonable model for the internal structures of low-mass MS stars. The values are tabulated in Table 1 of Mardling & Aarseth (2001); for the l = 2 f −modes, we have w₁ = w₂ = 8.175 in units of $G{m}_{*}/r_{*}^3$ and Q₁ = Q₂ = 0.2372. q = m₂/m₁ is the mass fraction where m₁ is the dominant star; since all our stars are of the same polytrope, the dominant star that absorbs most of the tidal energy is always the one with the largest radius, r₁, (Mardling & Aarseth 2001). We assume conservation of orbital angular momentum $p=p_0\frac{1+e_0}{1+e}$ to obtain the separation $a=\frac{p{r}_1}{1-e}$.

2.4.2 Chaotic tidal evolution

Low-mass MS stars and giants have convective envelopes. At very close periapsis separations, a coupling between the equilibrium tide and the internal convective motions of the star(s) may result in damping via a turbulent viscosity (Mardling 1995a,b; Mardling & Aarseth 2001). The theory presented in (Mardling 1995a,b) suggests that this tidal evolution becomes chaotic; the exchange in energy between the orbit and the stars is unpredictable and allows the tides to build up to a significant fraction of the binding energy of the star.

Binaries that are very eccentric with a periapsis radius of only a few stellar radii lie in what Mardling & Aarseth (2001) refer to as the chaotic region (see their Section 2.3). As the binary loses orbital energy and circularises, once enough energy has been lost, the binary leaves the chaotic region and proceeds to circularise normally (see Section 2.4.1 and e.g. Figures 3 and 8 in Mardling & Aarseth 2001).

For the purposes of this study, we are mostly interested in how much the binary circularises in the chaotic region. We assumed that the chaotic phase (typically ≳100 yr) is short compared to the time between encounters (on average ≳10⁴ yr for the binaries that undergo chaotic tides in our simulations) and hence that a binary that has entered the chaotic tidal regime can be treated as exiting instantaneously. Following Mardling & Aarseth (2001), we obtained the final eccentricity, e, from⁵ $$\frac{{\left(1+e\right)}^{2.4}}{1-e}{\left(1+e\right)}^5={\left(\frac{p_0\left(1+e_0\right)}{C}\right)}^5{\left(\frac{{\omega }^2}{1+q}\right)}^2,$$(9)

where the initial conditions are the eccentricity e₀ and the periapsis separation scaled by the stellar radius r₁, p₀ = a₀(1 − e₀)/r₁. Here q is the mass ratio of the secondary to the primary, 𝒞 = 1.836, and ω is the oscillation frequency structure constant for a poly trope of n = 3, values can be found in Mardling & Aarseth (2001, Table 1). The solution to Eq. (9) is obtained numerically, and we get the semi-major axis directly from the initial state by assuming constant orbital angular momentum.

3 Method

3.1 Initial stellar and binary population

The stars in the NSC are predominantly old, with roughly 80% being formed more than 5 Gyr ago (Pfuhl et al. 2011; Schödel et al. 2020; Chen et al. 2023). Furthermore, Pfuhl et al. (2011) found that the old population of giant branch stars mainly consists of stars corresponding to zero-age MS masses of 1–1.2 M_⊙, which would put solar mass stars right at the MS turn-off. We are interested primarily in binaries that may be contributing to the old visible population today, which from expectations would be binaries with m₁ ≲ 1 M_⊙ (Pfuhl et al. 2011; Schödel et al. 2020) and located near the hard/soft boundary (H/SB, i.e. h ~ 1, see Eq. (3). See, also e.g. Panamarev et al. 2019; Stephan et al. 2019; Rose et al. 2020). Additionally, the evolution of binaries near the H/SB is the most uncertain (consider, for example, Heggie’s law), which makes the employment of direct N-body methods very useful if not necessary. Consistent with the turn-off mass, we take the age of the NSC to be 10 Gyr and the primary to be between 0.9–1.0 M_⊙.

The mass ratio, period (thereby also separation), and eccentricity for the binaries are all sampled from the distributions in Moe & Di Stefano (2017) using the publicly available package COSMIC (Breivik et al. 2020). We consider secondary masses between 0.1–1.0 M_⊙. We know very little about the binary population in the NSC and especially how a primordial population might have looked like. A binary population in the NSC with properties similar to the binaries found in the solar neighbour-hood (see Raghavan et al. 2010; Moe & Di Stefano 2017) would, for the majority, be considered soft in the NSC due to the high velocity dispersion (see Eq. (3)) and should inevitably evaporate (Prodan et al. 2015; Stephan et al. 2016, 2019; Panamarev et al. 2019). More generally, binary period distributions in dense stellar systems are poorly constrained even for present day distributions (see, e.g. Ji & Bregman 2015; Müller-Horn, Johanna et al. 2025). The present day population in the NSC should have a narrower distribution than in globular clusters because of the higher number density and velocity dispersion, but not much can be said regarding the primordial population. Because of the lack of observational constraints and theoretical guidance, we assume the primordial binary population to resemble the field population of short to relatively long period binaries. We put the upper limit on the binary periods at 1000 days and let the lower limit remain at the default value of 1.4 days (for a 2 M_⊙ binary this roughly corresponds to 2.5 au and 0.03 au, respectively); harder binaries are going to be extremely rare and will likely merge within a relatively short time, and softer binaries should evaporate rather quickly (see e.g. Stephan et al. 2016). Figure 1 shows the initial distribution of hardness ratios at 0.1 and 0.3 pc in black and red solid lines, respectively. The vertical dotted line marks the H/SB.

The stellar and binary properties ensure that we can avoid complex stellar evolution. Instead, we set up a mass-radius interpolation scheme using stellar tracks produced with STARS (Eggleton 1971; Pols et al. 1995) of masses between 0.1–3.0 M_⊙ to also accommodate potential mass gain in collisions. In cases of collisions, we assume mass and momentum conservation and no mixing; the remnant’s MS lifetime is determined by interpolating the helium mass fraction at the time of merger for the binary star to the corresponding age for the remnant. The helium fraction interpolation scheme is also produced using masses and ages from stellar tracks produced with STARS.

Fig. 1 Initial (final) distribution of hardness ratio for binaries originating at 0.1 and 0.3 pc in solid (dashed) black and red lines, respectively. The vertical dotted line marks the H/SB. Each bin has been normalised to the total sum.

3.2 Evolving a binary

3.2.1 Stellar dynamical Interactions

We work under the assumption that the primordial binary fraction is sufficiently small that we can neglect binary-binary and higher-order encounters. In other words, we are only considering three-body encounters, in which the following outcomes are possible: (1) collisions, (2) scatterings, (3) exchanges, and (4) captures, where we form transient triples. It is in turn possible for the binary to be destroyed following one of these encounters, in which it can (5) ionise/evaporate or (6) merge. Additionally, because we only consider the mass of the SMBH and by extension Keplerian orbits, and no stellar evolution beyond the MS, we are forced to terminate the dynamical evolution of the binary if it (7) unbinds from the SMBH (ejected), (8) crosses the SMBH Roche limit, and (9) the primary evolves. We outline the conditions for each of the possible outcomes 1–9 to be fulfilled below. Subscripts 1 and 2 refer to the primary and secondary of the initial binary, respectively, and subscript 3 denotes the tertiary.

1. A collision occurs when r_p/(r_i + r₃) < 1, where r_p is the periapsis and r_i, r₃ are the radii of the two colliding stars (i = 1, 2), respectively. All collisions are assumed to conserve mass and momentum. While destructive collisions do occur in the GC, they are mostly confined to the inner 0.01 pc (Rose et al. 2023).

2. A scatter (also referred to as a flyby in the literature) is defined as resulting in a bound binary (E_1,2 < 0) while the tertiary is both unbound (E_12,3 > 0) and moving away (x₃ ⋅ v₃ > 0) from the binary centre of mass. If these two constraints are fulfilled and the tertiary is further than 10a, a being the instantaneous separation, we consider the scatter to be concluded.

3. An exchange is similar to a scatter but instead of the tertiary being unbound, one of the binary constituents (i) becomes unbound from the common centre of mass of the remaining two stars, E_i,j3 > 0, where i = 1, 2 and j = 3 − i. As before, if x_i ⋅ v_i > 0 and the distance between star i and the common centre of mass of the other two stars is greater than 10a, the exchange is fulfilled.

4. When a tertiary is captured by the binary we form a transient triple, which require all three stars to be bound to their common centre of mass, e.g. E_i,COM < 0 for i = 1, 2, 3. If we produce a transient triple, and the tertiary orbit is longer than 10³ days we assume the encounter to be a scatter; the tertiary will spend most of its time at apoapsis where it very easily should be unbounded. If the orbit is shorter than 10³ days we resolve the encounter by an additional 10 tertiary orbital periods, if at after this point nothing has happened to the binary we also assume a scatter (however, we find no such instances).

5. In an ionisation event all stars must be unbound, e.g. E_i,j > 0 for i = 1, 2 , 3 and j = 1, 2, 3 ≠ i. Furthermore, they also have to be moving away from the tertiary centre of mass x_i ⋅ v_i > 0 for i = 1, 2, 3. In order to reduce computational requirements we also consider binaries that have a > 100 au to inevitably evaporate (ionise) which is therefore also considered as a stopping condition.

6. A merger can be caused by the hardening of an encounter (e.g. scatter or collision), due to ZLK oscillations or due to Roche Lobe overflow. In the two former cases, the same stopping condition as in a collision is used, while for the Roche lobe overflow we instead require R_1,2 > f(q)a(1 − e), where f (q) denotes Eggleton’s formula (Eggleton 1983) for the Roche lobe and is given by⁶ $$f\left(q\right)=\frac{0.49q^{2/3}}{0.6q^{2/3}+\ln 1+q^{1/3}}.$$

7. Ejection implies that the binary itself unbinds from the SMBH and simply requires $E_{COM}^{\text{bin-SMBH}}>0$.

8. If the binary ventures inside the SMBH Roche limit it evaporates, with a possibility of a tidal disruption event (TDE). We require $\frac{a_2}{a_1}<{\left(\frac{3M_{•}}{M_{1,2}}\right)}^{1/3}\frac{1+e_1}{1-e_2}$.

9. If the binary becomes older than its primary’s expected MS lifetime, we terminate its evolution to avoid computing complex stellar evolution. The star’s MS lifetimes are interpolated from stellar tracks produced with STARS.

All termination constraints naturally work as criteria to run three-body interactions, ZLK, and tidal evolution.

3.2.2 Three-body encounters

The numerical integration of three-body encounters is performed using tsunami (Trani & Spera 2023). In total, we simulate more than 6 × 10⁷ encounters, with roughly 1 in every 150 000 encounters leading to something different than a scatter. The SMBH is not included in the three-body encounters, and as such it is only the interactions between the binary and the tertiary we simulate. Furthermore, we do not follow the ZLK evolution in tsunami, but rather integrate the quadrupole equations of motion described in Section 2.3. Finally, the implementation of the equilibrium tides described in Section 2.4 (see also Mardling & Aarseth 2001), is consistent with that implemented in tsunami. As such, even though the duration of any meaningful (equilibrium) tidal evolution typically far exceeds the time a binary spends in a tsunami encounter, we are not “missing” any contributions to its tidal evolution. The order in which these three mechanisms (three-body encounters, ZLK and tides) are considered is explained in Section 3.2.3.

For simplicity, we neglect the influence of the SMBH’s tidal field on the dynamics of three-body encounters in this work. The SMBH can impact encounter durations, merger rates, and the final orbital properties of surviving binaries, as demonstrated in studies on three-body scatterings of stellar-mass black holes in the presence of a SMBH (Trani et al. 2019b; Trani 2020; Trani et al. 2024). A comprehensive investigation of its effects on stellar encounters is left for future studies. Nevertheless, the binary’s orbit around the SMBH is allowed to change: we updated the outer orbit with the velocity kicks received in any three-body encounter. Hereafter, we refer to this process as migration.

3.2.3 Routine

We introduce a hierarchical setup that is governed by the timescales of the different processes. The ZLK quadrupole timescale t_quad is given by Eq. (4), the instantaneous timescale for encounters t_enc, is given in Appendix A by taking the inverse of Eq. (A.2), and for a tidal dissipation timescale t_circ we refer the reader to Eqs. (20)–(25) in Mardling & Aarseth (2001).

The basis of our routine lies in the three-body encounters; only these are resolved with tsunami, and the secular evolution (ZLK and tides) is computed in-between encounters. The hierarchy of the setup is as follows:

• Chaotic tidal evolution: if a binary lies within the chaotic region, its tidal evolution is almost always fast enough that tidal processes dominate (Mardling & Aarseth 2001). We neglected the time taken for these dynamical tides to damp out and followed the procedure laid out in Section 2.4.2 to obtain the post-chaos parameters.

• Secular evolution: next we compared the timescales for the ZLK mechanism (t_quad; see Section 2.3) and the equilibrium tide (t_circ; see Section 2.4). The secular process with the shortest timescale was selected.

• Close encounters: we obtained the instantaneous encounter timescale t_enc that encompasses all dynamically significant encounters (see Appendix A).

We first sampled the binaries with COSMIC and placed them on an orbit around the SMBH with orbital elements drawn according to the distributions in Table 1. The impact parameter, velocity, and mass for the tertiary were drawn as described in Appendix A.1. We calculated a time step δt = min(t_quad, t_circ, t_enc) and computed the secular process with the shortest timescale for a duration of δt given that⁷

• i)

  In the case of circularisation, t_circ < 1 Gyr and e_inner > 0.002. If t_circ > 1 Gyr and the eccentricity ≤0.002, so the binary is declared circularised (similar to Mardling & Aarseth 2001).

• ii)

  In the case of ZLK, the mutual inclination must be large enough for the torques to arise (≳40 deg), and GRP cannot dominate (e.g. Eq. (5) must be smaller than unity).

At fixed intervals of 0.1 ∆t we sample a three-body encounter with the accept-reject method described in Appendix A: if rejected, the secular evolution continues and we take another step 0.1 δt; if accepted, we resolve the interaction with tsunami and determine the outcome based on the conditions defined in Section 3.2.2. For outcomes 1–4, we update the binary’s inner and outer orbits, and the process is then restarted. For outcomes 5–9, the dynamical evolution of the binary is terminated. A flow chart of the routine is included in Figure 2, which also shows a sketch of the setup and the different outcomes.

Fig. 2 Flow chart of our routine described in Section 3.2.3. The central sketch shows the overall setup of the three-body encounters; a binary (m1, m2) on an orbit at a distance of d• around the SMBH (MSMBH) interacting with a tertiary (m3). To the right of this sketch we show the possible outcomes from a three-body interaction. Group C contains interactions in which the binary remains intact and bound to the SMBH and as a consequence will experience continued evolution. Conversely, the evolution is terminated by outcomes in group T. To the left of the sketch we illustrate how a binary is evolved. A timestep, ∆t, is set by the shortest timescale of ZLK (quadrupole approximation), encounters, and circularisation (tides). The secular process with the shortest timescale min[τquad, τcirc] is then run for ∆t. At every step 0.1 ∆t, a three-body encounter is sampled with the accept-reject method described in Appendix A. If accepted, the three-body encounter is carried out with tsunami. If rejected, we proceed to the next 0.1 ∆t and sample again. After a time, ∆t, has passed or after a three-body encounter, a new ∆t is calculated. In this example, the encounter occurs after ∆t1 + 0.1 ∆t2 + tenc. This process is then continued until an outcome in group T occurs or until the simulation time exceeds 10 Gyr.

4 Results

Figure 3 shows an overview of the outcomes from our simulations, and we base the following discussion on these results. The upper (bottom) row refers to the binaries originally set at 0.1 pc (0.3 pc) of the SMBH. The left-hand side column shows the inner semi-major axis while the right-hand side column shows the outer semi-major axis. The black lines correspond to the initial conditions; for the outer orbits these are vertical lines with a height equivalent to the total number of binaries considered at respective distance (1654 at 0.1 pc and 4082 at 0.3 pc). Also, see Table 2 for the total number number of binaries for each outcome. In the subsections below, we go through Figure 3 line by line, starting with mergers and evaporation (orange and green lines, respectively), followed by the ejected (magenta), intact (red) and evolved (blue) populations.

A merger is the most common outcome for binaries starting with a_• = 0.1 pc (46%) and approximately equally as common as evaporation at a_• = 0.3 pc (29%)⁸. This is perhaps the most clear in Figure 4, which shows the cumulative fraction as a function of time for the various outcomes. The orange boxes in Figure 3 show the distributions of the semi-major axes of merging binaries just before the merger. It can be seen that the merging process is not dominated by hardening of the binaries – the most common semi-major axis for merger is in both cases similar to the initial most common value around 1 au. Although a fraction of mergers occur in very tight binaries in which the semi-major axis a ~ R_1,2, where R_1,2 is the sum of the stars’ radii, the majority of the merging binaries are much wider. More explicitly, it is not a decrease in separation that drives mergers, but rather a rapid increase in the binary’s eccentricity. The cumulative distribution in eccentricity for the different outcomes in Fig. 5 shows this quite clearly. The initial distribution is weighted towards more circular orbits with a median initial eccentricity of around one-third. The corresponding best-fit eccentricity distribution P(e) ∝ e^η has η < 0. However, the mergers are strongly superthermal, with only around a third of merging binaries having e < 0.9 (best fits for the eccentricity distributions of all outcomes are tabulated in Table C.1).

Fig. 3 Left: inner semi-major axis distribution of the initial population (black) and the corresponding final outcomes (colours). Yellow lines refer to binary mergers, green lines to binaries that evaporate, magenta lines to binaries that are ejected (unbound from the SMBH), red lines to intact binaries (i.e. binaries that survive to 10 Gyr), and blue lines to evolved binaries (stellar evolution of the primary occurs before binary destruction and by 10 Gyr). In the cases of evaporation and mergers, the properties shown in the figure refer to the binary just before either destruction mechanism takes place. Ejection corresponds to the region the binaries were ejected from and not where they end up. The top row corresponds to binaries starting with a• = 0.1 pc from the SMBH, while the bottom row corresponds to a• = 0.3 pc. Right: outer semi-major axis distributions. The vertical black lines represent the initial populations, and the boxes of the final outcomes are coloured in the same way as in the left panel.

Fig. 4 Cumulative distribution of outcomes as a function of time. The outcomes are colour-coded the exact same way as in Figure 3. The solid (dashed) lines refer to binaries initially 0.1 pc (0.3 pc) away from the SMBH.

4.1 Mergers

This result is not unexpected, since the binaries suffer ZLK oscillations which drives the eccentricity to values near unity (Prodan et al. 2015; Stephan et al. 2016, 2019). However, we also find that on average, repeated scatters contribute about as much as ZLK does in driving binaries to merge. The left-hand side panel of Figure 6 shows a normalised histogram of the difference between the total change in eccentricity due to ZLK and that due to repeated scatters for merging binaries. That is, for a given binary, we compute the total change in eccentricity solely caused by ZLK, and then subtract the total change solely caused by scatters. The result indicates that ZLK and scatters on average contribute equally in driving up the binary’s eccentricity until it eventually merges. The right-hand side panel instead shows a scatter plot of these two terms; ZLK along the x-axis and repeated scatters along the y-axis. A typical binary will increase its eccentricity through both processes. In some cases, however, the two processes work against each other; one process causes a net decrease (the binary becomes more circular) while the other causes a correspondingly larger net increase.

We also observed some clustering around the zero lines from the right-hand side panel. These are binaries that on average undergo almost no change from one of the processes and as such are dominated by the other. Figure 7 shows the initial (top row) and final (bottom row) separation versus 1 − e for these binaries. For scatter-dominated mergers (|ΣΔe_ZLK ≤ 0.05, shaded area along the vertical zero line), the binaries are initially quite hard and circular. Although many of these binaries become softer and more eccentric, the bulk of them remain hard and relatively circular. They do not undergo ZLK oscillations because GRP suppresses ZLK for hard and tight binaries, and as such, their mergers are driven only by scattering.

Conversely, the ZLK dominated mergers (|ΣΔe_Scatter| ≤ 0.05 – shaded area along the horizontal zero line) are initially quite wide and soft with moderate eccentric orbits and also remain wide and soft, though naturally become more eccentric due to ZLK oscillations. Wide binaries should (and do) suffer scatters at a higher rate than tight binaries, but the timescale at which these scatters induce a merger is typically an order of magnitude longer than ZLK. Figure 8 shows the time distribution of when the mergers shown in Figure 7 occur. On average, a scatter-dominated evolution (left-hand side panel) leads to mergers after 10⁹ years, while for a ZLK dominated evolution, it only takes <10⁸ years.

To demonstrate that there are two distinct groups of mergers, we plot the two-dimensional distributions of the inner semi-major axis α and eccentricity (1 − e) grouped by outcome in Figures 9 and 10.

In panel a, the total initial distribution is shown in black, mergers are shown in orange in panel b, evaporated binaries are shown in green in panel c, and the remaining panels d–f in the bottom row, show ejected binaries in magenta, intact binaries in red, and evolved binaries in blue, respectively. Each panel also has three contours (10%, 50%, and 90% of the maximum number of counts) corresponding to the initial distribution for the binaries in that panel. The panel for mergers shows two distinct regions. The first region, with about two-thirds of the mergers, contains the initially wide binaries whose mergers are driven by an increase to their eccentricities, mostly due to ZLK. The second is the merger of initially tight binaries, driven mainly by scattering, which then merge over a wide range of eccentricities.

Fig. 5 Cumulative distribution of the eccentricities for the inner (panel a) and outer (panel b) orbit, colour-coded as in the previous figures (outcome). The dotted grey lines follow specific power-law distribution eη where η = −1 correspond to mostly circular orbits, η = 0 to a uniform distribution and η = 1 to a thermal distribution (e.g. mostly eccentric orbits). As in Figure 3, the coloured distributions consider the final state of a binary where it can still be considered a binary (and not destroyed). Best-fit parameters can be found in Table C.1.

Fig. 6 Panel a: distribution of the differences between changes in eccentricity due to ZLK and due to repeated scatters (only looking at binaries that eventually merge). The black line corresponds to binaries initially 0.1 pc away from the SMBH and the red line to 0.3 pc. Panel b: scatter of these two terms. One point corresponds to one binary. The zero-lines are shown as dotted black lines, and the shaded areas correspond to the ±0.05 region around the zero-lines. The two diagonal lines show x + y = 0 (lower) and x + y = 1 (upper).

Fig. 7 Four-panel plot showing the separation versus 1 − e of the inner orbit. Based on the zero-lines in Figure 6. The left-hand side column corresponds to the vertical zero-line (|ΣΔeZLK| ≤ 0.05), and the right-hand side column corresponds to the horizontal zero-line (|ΣΔeScatterl ≤ 0.05). The top row shows the initial distribution at t = 0, while the bottom row shows the final distribution at t = tf before the binaries merge.

Fig. 8 Time of merger distribution for the binaries in Figure 7. Panel a, as before, refers to the vertical zero-line (effectively no ZLK) in Figure 6, while panel b refers to the horizontal zero-line (effectively no scatters). The black (red) line is for binaries initially at 0.1 pc (0.3 pc) away from the SMBH.

4.2 Evaporation

Evaporation (green lines in Figures 3 to 5) is slightly more common than merger at 0.3 pc (30%) and somewhat less common at 0.1 pc (35%). In particular, evaporation dominates for wider inner orbits ≳1 au⁹; beyond ~100 au the binary unbinds as it exceeds its Hill radius due to the SMBH. For smaller initial separations, it is predominantly scatters that unbind the binary, either via a strong encounter that ionises it or through repeated weaker encounters that gradually inject energy. The green contours in Figures 9 and 10 (panel c) show the distribution of binaries just before they evaporate; the smaller the semi-major axis is the stronger the encounter that causes the binary to unbind needs to be. In other words, if a binary evaporates at 0.1 au, it is ionised by a very strong encounter, and if it evaporates at, for example, 10 au, then it is largely due to a gradual process.

The corresponding final outer orbits seen in Figure 3 are centred around the initial distances (0.1, 0.3 pc) and are distributed rather symmetric in log space. The ratio of evaporated and merged binaries is roughly constant throughout the simulations, as seen in Figure 4. Because both processes are fast, most binaries will therefore not have had time to migrate in either direction before being destroyed, which explains the shape of the distribution in Figure 3¹⁰. For the inner orbits, evaporations follow a very similar distribution to the merged binaries, but the distribution extends to very wide separations. The evaporating binaries have an eccentricity distribution that is close to thermal (η = 1). Eccentric binaries are more vulnerable to being broken up than circular binaries, since they spend much of their time at apoapsis, where a smaller velocity perturbation is sufficient to unbind them.

4.3 Hard binaries

Binaries that have not merged or evaporated fall into three categories: ejected binaries (magenta in Figs. 3 and 5), intact binaries (red) and evolved binaries (blue). All three groups have similarly distributed inner semi-major axes, eccentricities, and are relatively hard; a < 1 au but peaks at ≳0.1 au.

The ejection of a binary occurs when it becomes unbound to the SMBH. We do not attempt to follow the subsequent dynamical evolution of such binaries, although they will likely remain within the NSC on wider galactic orbits. On these orbits the number densities are smaller, so they do not undergo significant further close encounters. Furthermore, GRP should inhibit ZLK oscillations even for the wider binaries (see Eq. (5)). Mass precession, induced by the extended potential of the cluster, may also suppress ZLK; though only at the octupole order¹¹ as the quadrupole order is independent of the outer argument of peri-apsis (see e.g. Naoz et al. 2013; Stephan et al. 2016; Tep et al. 2021). We therefore expect that these binaries will either evolve, if their primary star is more massive than the turn-off mass at the end of the simulation, or remain intact.

Similarly, the intact binaries are all found to have migrated away from the SMBH. This is counter-intuitive at first sight, since binaries, which on average are more massive than single stars, preferentially migrate inward due to dynamical mass segregation. However, binaries that migrate inward are consistently destroyed, whereas those that migrate outward are not; the quadrupole timescale becomes longer and the interaction rate with surrounding stars decreases (see, Eq. (A.3)). Moreover, and partly for the same reason, surviving binaries are preferentially hard and circular (see dashed lines in Figure 1; virtually all binaries have h ≥ 1. See also Figures 3, 5, 9 and 10). In comparison, three-body encounters involving hard binaries in globular clusters, would typically lead to a thermal or super-thermal distribution in eccentricity (Heggie 1975a,b; Stone & Leigh 2019; Ginat & Perets 2023; Rando Forastier et al. 2025). Coupled with the migration discussed above, very hard and circular binaries do not undergo ZLK nor do they interact much with surrounding stars, and the ones that do get eccentric likely merge or circularise as a consequence. Since the hardness of a binary is a relative term (Eq. (3)), a hard binary in the NSC will have a very small separation¹², which further limits the possible eccentricities that do not lead to a merger or the circularisation of the binary. For example, the most common separation for our binaries at 0.1 pc is roughly 0.05 au ≈ 10 R_⊙ (see Figure 3), which means that for an eccentricity >0.8 the binary is likely to merge and otherwise likely to circularise (Mardling & Aarseth 2001). As such, we do not expect a thermal distribution for old primordial binaries.

The evolved population is essentially just a subset of the intact population for a few reasons. In these simulations we follow the evolution of sub-solar mass stars, meaning all of them are on the MS up until the very end of the simulation (≳9 Gyr as seen in Figure 4). We do not compute any complex stellar evolution during the simulation (see details on stellar evolution in Section 3.1), and as such this distinction between evolved and intact binaries is essentially just separating the most massive stars (≳0.97 M_⊙) from the rest. Since the majority of binaries are destroyed within the first few billion years, we would not expect a substantial difference between these two populations. The main thing that separates the evolved population from the intact population is its small fraction of binaries that have evolved prematurely (≲9 Gyr) due to stellar collisions.

Fig. 9 Contours of the inner semi-major axis versus 1 − e for various outcomes at 0.1 pc from the SMBH. Upper row, from the left, initial distribution in black (a), merged binaries in orange (b), and evaporated binaries in green (c). Bottom row, from the left, the intact population (after 10 Gyr) in red (d), the evolved population in blue (e), and ejected binaries in magenta (f). For panels (b), (c), and (d), the distribution refers to the binary’s properties the moment just before the final outcome.

Fig. 10 Same as Figure 9 but at 0.3 pc.

4.4 Collisions

With the high number densities of stars, stellar collisions are relatively common in the NSC. Here we present the sample of collisions from our simulation in more detail.

The upper panel of Fig. 11 shows the impact parameter b versus the velocity at infinity v_∞ for the tertiary in cases where it leads to a collision. The v_∞ velocities range from a few tens to a few hundreds of km s⁻¹. As such, none of these collisions should be able to strip the binary star with which it collides (Rose et al. 2023). Therefore, treating collisions with no mass loss and conservation of momentum should be a reasonable approximation (see Section 3.2.1). The lower panel of the Figure shows the relation between the impact parameter and the binary separation for the collisions; we find a scatter around the 1:1 ratio. This is unsurprising since the impact parameter is defined as the distance from the binary COM, b ~ a roughly corresponds to where one of the binary stars is expected to be.

A consequence of the binaries’ tendency to migrate inwards or outwards from the SMBH, is that the collision rate changes significantly compared to a fixed orbit around the SMBH. The solid lines in Fig. 12 show the average number of collisions per binary as a function of time. The dotted lines show the expected number of collisions for a single 1 M_⊙ star at a fixed galactocen-tric radius – e.g. after 10 Gyr we would expect a 1 M_⊙ star to have suffered at least one collision at 0.3 pc and between 2–3 at 0.1 pc. We are almost always an order of magnitude below this expectation, which is caused by the binaries’ migration (e.g. the dynamical recoil from three-body encounters, see Section 3.2). The dashed lines show the same expectation as the dotted lines, but where the distance to the SMBH changes as the average of the binaries’ outer orbits (normal black/red for semi-major axis, and pale black/red for apoapsis). The way in which the average outer orbit varies as a function of time is shown in Figure 13. After a few Gyr, the yet intact binaries have on average migrated well beyond their starting points, and as a consequence, most collisions with a binary (which there are fewer and fewer of) also occur further out in the NSC. The transparent dashed lines in Figure 12, which follow the average apoapsis of the outer orbits, is the only expectation that translates to the actual rate of collisions per binary; binaries spend most of their (outer) orbits at apoapsis, leading to a much lower collision rate.

Figure 14 shows the time of the final outcome t_life of binaries that have suffered at least one collision versus the time of collision t_coll. The black dotted line marks where t_life = t_coll; for a binary on this line the outcome is a direct consequence of the collision. The percentage for each outcome is also indicated in the legend of the Figure; most commonly, a collision leads to a merger between the remaining stars, and in most cases this also happens as a direct consequence of the collision. We refer to these mergers as three-body pile-ups (3BPUs).

The second most common outcome is the intact + evolved + ejected outcomes, which are covered in Section 4.3. These collisions have had little effect on the subsequent evolution of the binaries and so the final outcome is usually significantly later than the collision. The same goes for the small fraction of binaries that evaporate – while they are not a direct consequence of the collision, they tend to happen very shortly afterward.

Fig. 11 Properties for collisions between a tertiary and one of the binary stars. Panel a shows the relation between the tertiary’s velocity at infinity and its impact parameter. The points are colour-coded based on their initial position from the SMBH. Panel b shows the impact parameter versus the binary’s semi-major axis before the collision. To help guide the eye, the dotted black line corresponds to ainner = b.

Fig. 12 Average number of collisions per binary as a function of time (solid lines). The dotted lines refer to the expected number of collisions a single 1 M⊙ star would have at 0.1 pc (black) and 0.3 pc (red). Because the binaries tend to migrate, the dashed lines show the same expectation as the dotted line, but where the distance is not fixed to the starting positions and instead follows the average outer orbit of the binaries. The pale dashed lines follow the apoapsis, while the non-transparent dashed lines follow the semi-major axis.

Fig. 13 Mean outer semi-major axis for all remaining binaries as a function of time. Black lines show binaries starting at 0.1 pc from the SMBH; red lines show binaries starting at 0.3 pc. Each point shows the time and location of a collision.

Fig. 14 Time of final outcome versus the time of collision for binaries that have suffered at least one collision. The dotted black line shows where tlife = tcoll; points on this line meet their final outcome as a direct consequence of the collision. The points are colour-coded based on the outcome similarly to Fig. 3.

5 Discussion

In the previous section, we give an overview of the many different outcomes a binary in the NSC may suffer, and the various processes that inevitably steer the binary’s evolution towards one of the outcomes.

The parameter space that we consider for the binaries in these simulations is quite stringent, which makes the varied distribution of outcomes interesting; the path a binary takes seems to be largely stochastic. While processes such as ZLK oscillations are well behaved (in the sense that you can predict the oscillations), encounters with surrounding stars are not. The consequences of this are particularly apparent for binaries at (or near) the H/SB. For example, as seen from Figure 4, at 0.3 pc a binary (from our sample) is essentially as likely to remain intact as it is to merge or evaporate. Figures 9 and 10 also show that the initial distribution in separation and eccentricity is very similar between the outcomes; the contours are virtually the same for binaries that eventually merge and evaporate. The 90% contours, which show two distinct regions in both cases, only differ in extent between the two outcomes. Evaporated binaries are shifted towards greater separations (these are soft and not near the hard/soft limit) while the merged binaries are shifted towards smaller separations and look virtually identical to binaries that remain intact, evolve or get ejected. As discussed in Section 4.3, the only real difference between the intact and evolved binaries is the mass, but in the end, the outcome will still be the same. Although unpredictable, binaries that are moderately soft to soft tend to either merge or evaporate (latter is more likely for softer binaries), binaries that are moderately soft to moderately hard merge, remain intact (evolve), or get ejected. What outcome that occurs seems to be largely due to chance.

5.1 Binary mergers and the S-cluster

Binary mergers are perhaps the most interesting aspect of the outcomes in a binary’s evolution in the NSC. The leading theory for the origin of the G-objects (Ciurlo et al. 2020) is that they form as a consequence of binary mergers (Prodan et al. 2015; Stephan et al. 2016, 2019; Ciurlo et al. 2020; Peißker et al. 2024), and that they may be related to the lower-mass population of the S-cluster (Melamed & Peißker 2024).

Figure 15 illustrates various properties of the post-merger binaries. The solid black lines represent mergers of binaries initially located at 0.1 pc from the SMBH, while the solid red lines correspond to those initially positioned at 0.3 pc. The top panel depicts the distribution of the merger products’ total lifetimes; that is, the time spent in the binary plus the additional time gained from the merger ‘rejuvenation’¹³. The dotted black line marks 10 Gyr – the time considered for our simulations, which roughly corresponds to the age of the NSC. A small fraction (~1%) of merger products are located to the right of this line, suggesting a few mergers originating from the early populations of binaries in the NSC could be observed at present day either as an MS star or as a recently evolved star. Of course, the shape of the distribution(s) is largely due to our choice to focus on binaries where the primary is near the MS turn-off mass of the NSC (Pfuhl et al. 2011; Schödel et al. 2020; Chen et al. 2023); considering more massive stars would push the distribution towards the left¹⁴, and less massive stars would push it to the right. In other words, the same mechanism can make massive blue straggler stars (BSSs) at earlier times and less massive ones today.

The middle panel shows the distribution of the merger products’ orbits around the SMBH, where the black dotted line and red dotted line correspond to the initial positions for the binaries: 0.1 pc and 0.3 pc, respectively. Most of the post-merger orbits remain at these distances, as can be seen from the two distinct peaks. Not surprisingly, the 0.1 pc peak is much narrower; the relevant timescales are shorter, which, apart from causing the binaries to merge more quickly (see e.g. Figure 4), also means that the binaries do not have a lot of time to migrate (in either direction). Still, there is a small fraction (~1%, all of which originate from 0.1 pc) of binaries that have migrated inward towards the SMBH, obtaining an outer semi-major axis smaller than about 1 arcsecond (≈0.04 pc), see also Figure 3.

This is roughly the region that contains the orbits of the S-cluster, including the G-objects and dusty sources (Habibi et al. 2017; Schödel et al. 2020; Peißker et al. 2024; Melamed & Peißker 2024). The results therefore support the idea that the S-cluster and G-objects, may originate further out in the NSC (see also, for example, Verberne et al. 2025). However, in our case, there is a mismatch in timescales and age: the S-cluster and the G-objects appear to be young ≲10–15 Myr (Habibi et al. 2017; Schödel et al. 2020; Peißker et al. 2024; Melamed & Peißker 2024), which is comparable (in the context of the G-objects) to the ZLK merger timescale for binaries on orbits a_• ≲ 0.1 pc (see e.g. Stephan et al. 2016, 2019). Compared to our simulations, the binary mergers end up within the inner arcsec after at least ~300 Myr with the mean being 660 Myr, see Appendix E. Although there is a degeneracy between stellar age and metallicity for S-cluster stars (Generozov et al. 2025), it is not strong enough to account for the differences seen here. Furthermore, with the recent observation of a lower-mass spectroscopic binary system in the NSC (Peißker et al. 2024), estimated to be only a few Myr old and most likely a progenitor of the G-objects, another mechanism that is not in tension with the recent star formation in the GC (Lu et al. 2013; Schödel et al. 2020) is required (see for example Genzel et al. 2010; Trani et al. 2019a; Generozov et al. 2025; Verberne et al. 2025).

However, Figure 15 shows that the mergers of binaries originating between 0.1–0.3 pc from the SMBH can contaminate the region where the S-cluster is. Although most stars in the NSC are old (≳5–8 Gyr), smaller fractions of stars are associated with star formation events at different epochs (Schödel et al. 2020). The second most prominent event in terms of the mass fraction formed per Gyr is believed to have occurred 2–4 Gyr ago (see, for example Schödel et al. 2020). As seen in the top panel of Figure 15; the peak of the total expected lifetime for the mergers is around 2 Gyr for binaries originating both at 0.1 pc (black) and 0.3 pc (red), suggesting that there could be a population of BSSs at or just beyond the MS turn-off scattered in the NSC (see middle panel) associated with this star formation event (SFE).

Assuming the NSC is ~10 Gyr old (dotted black line in the top panel), virtually all mergers from the oldest population of stars in the NSC would have evolved at present day due to single star stellar evolution. Thorsbro et al. (2023) reported metallicities for three ~Gyr old stars in the NSC. One of these, which has an apparent age of 1.5 Gyr and a very low metallicity of [Fe/H] = −1.2, was suggested to be a potential BSS. In their hypothesis (which was not presented as a definite model), this BSS would have formed from an initially hard (a ~ 6 R_⊙) and circular binary with stars of 0.9 and 0.8 M_⊙ – very similar to the binaries considered in our simulations. Following a phase of stable mass transfer at 8.3 Gyr, the binary would eventually merge at an age of 9 Gyr. In their test to this hypothesis, the binary is evolved in isolation, e.g. not in the environment of the NSC. Based on Rose et al. (2020) and their analytical study of binaries in the inner ≲0.5 pc, such a binary could remain intact for roughly that time period outside ≳0.3pc¹⁵.. Figures 9 and 10 show that such hard and circular binaries can survive both from 0.1 pc and 0.3 pc for up to 10 Gyr in the NSC despite its harsh environment, thanks to outward migration. Alternatively, the binary could originally have been both wider and more eccentric, before circularising. As such, the premise that this star may originate from a binary merger associated with an older population of the NSC seems perfectly reasonable. Additionally, it also suggests that more stars similar to it should be present in the NSC; at present day such stars could be located along the asymptotic giant branch on the HR-diagram (where the star in Thorsbro et al. 2023, was observed).

The bottom panel in Figure 15 shows the distribution of masses for the merger-products. A small fraction of these have masses above 2 M_⊙; since the maximum mass for each of our sampled stars are 1 M_⊙, the only possible explanation is that an additional star (tertiary) is included via a collision. We find that when a physical collision occurs between a star in a binary and an incoming single star, in 80% of the cases the binary merges as a consequence. As seen in Fig. 14, most of these mergers are also instantaneous in the sense that they occur as a direct consequence of the collision. These 3BPUs serve as a way to form more massive BSSs.

The top panel of Figure 16 shows the rate $\frac{N_{3\text{BPU}}\left(t\right)}{N_{\text{bin}}\left(t\right)\times t}$ of 3BPUs, normalised to the number of intact binaries, as a function of time. Although the mean 3BPU rates over 10 Gyr are 0.0034 Gyr⁻¹ and 0.0032 Gyr⁻¹ for binaries originating at 0.1 pc and 0.3 pc, respectively, most 3BPUs occur within ≲1−2 Gyr. Given the recent star formation events (Schödel et al. 2020), it would not be implausible to observe a 3BPU merger product in the more central regions of the NSC. Directly observing a 3BPU would be unlikely, though less unlikely than observing a normal merger or collision, which is already extremely rare. Also, when so close to the SMBH, ZLK-induced collisions may deplete the available sample of binaries too quickly for slower scatter-induced 3BPUs to occur.

The middle panel of Figure 16 shows the difference between the actual age t_f (e.g. the time of the 3BPU) and the apparent age τ_age of the 3BPU product, using the helium core mass fraction to infer the age. This assumes there is no mixing, and that the new age is always determined by the most massive star. The difference should be an estimate of how much younger (or older) the 3BPU product would appear if observed. Although we estimate that most of the 3BPUs would appear a few Gyr younger, we also estimate differences almost as high as 10 Gyr. Such a discrepancy could suggest observations of very young and metal-poor stars in the NSC that are similar to, if not even more extreme than, the BSS candidate in Thorsbro et al. (2023).

Finally, the bottom panel of Figure 16 shows the distribution of the final masses. Apart from the more massive stars (>2 M_⊙), the 3BPUs populate the region just below 2 M_⊙ as well, which in part explains the peak seen around this value in Figure 15. The details of the 3BPU merger products from these simulations are rather limited. We assumed that the initial collision conserves mass and angular momentum, which is not a bad approximation since we do not expect to have any collisionally stripped stars (see Rose et al. 2023, and Figure A.5), but means that the mass distribution is an upper limit. Secondly, grazing collisions, where the periapsis approach is similar to the sum of the stellar radii (r_p ~ R_ij), do not necessarily result in the coalescence of the two stars; there would simply be no 3BPU, and as such the fraction of 3BPUs is also an upper estimate. We also made the assumption of no mixing for stellar collisions and mergers. The amount of mixing that occurs as a result of the coalescence of two stars depends on various mechanisms: the relative velocity of the colliding stars, the impact parameter, the stars’ initial masses and chemical compositions (hence also evolutionary stage), the mass loss of the collision etc. (see e.g. Sills et al. 1997; Gibson et al. 2025). Head-on collisions between stars in globular clusters, with masses similar to those considered in our simulations, have been shown not to induce additional convection, thereby limiting the amount of fresh burning material transported into the core region (Sills et al. 1997; Sills et al. 2009). Sills et al. (2001) showed that BSSs that rotate typically remain longer on the MS than their non-rotating counterparts. Off-axis collisions can increase the rotational speed of the collision product, which in turn leads to rotation-induced mixing.

The masses and stellar lifetimes of our collision and merger products are undoubtedly affected by these mechanisms, but should not affect our results in the broader context. The middle and bottom panels of Figure 16 are only rough estimates. The precision of our results could be enhanced by making hydrody-namic simulations of the stellar mergers and incorporating stellar evolution.

Fig. 15 Panel a: normalised histograms for the total expected lifetime of the mergers. The dotted black line marks the age of the NSC at 10 Gyr. Panel b: post-merger semi-major axis around the SMBH. The dotted vertical lines mark the initial separations. The solid vertical line marks 1 arcsecond. Panel c: total mass after the merger. For all panels, binaries originating at 0.1 pc (0.3 pc) are shown in black (red).

Fig. 16 Panel a: rate of 3BPUs as a function of time. The dotted lines correspond to the mean rates over a 10 Gyr period. Panel b: normalised distribution of the estimated apparent age for the 3BPU merger products. The apparent age is determined by mapping the Helium-core mass fraction of the primary at the time of the merger to the corresponding age a normal star with the mass of the merger product would be with an identical Helium-core mass fraction. Panel c: normalised distributions of the total masses of the 3BPU merger products. For all panels, binaries originating at 0.1 pc are shown in black, and those from 0.3 pc are in red.

5.2 Comparisons to other studies

We have previously mentioned a few studies that also looked at binaries and their evolution in the NSC. Rose et al. (2020) conducted an analytical where the shortest timescale determines the fate of a binary, and Panamarev et al. (2019) carried out the largest N-body simulation of the NSC to date. Below, we limit the discussion to these two studies, as they are the most relevant to ours. We highlight their most relevant results and how they compare with ours.

Our findings are overall in good agreement with Rose et al. (2020); initially softer binaries evaporate, and those that are not, typically remain stable against evaporation. We find that the rate of evaporation is somewhere in-between their lower and upper limits (see e.g. the evaporation and max-evaporation curves in their Figure 1), with a typical evaporation timescale of a few ~100 Myr. Binaries that do not evaporate, as in those that initially are marginally soft or hard, will most likely merge (or alternatively migrate inward where they can evaporate) on a timescale set by the ZLK quadrupole timescale or the segregation timescale, which agrees with our results, see Figure 8. Rose et al. (2020) suggest that migration plays a secondary role to evaporation in shaping the binary distribution of the NSC. Although we also find that evaporation dominates (along with mergers), we would like to put extra emphasis on the migration’s role in shaping the binary distribution in the NSC. For example, the widest possible separation a 2 M_⊙ binary could have after 1 Gyr at 0.3 pc is estimated to be ~0.1 au (Rose et al. 2020, Fig. 7) and is set by the evaporation timescale. We, however, find a handful of such binaries at the same distance even after 10 Gyr, where some originate from 0.1 pc.

Furthermore, binaries may oscillate stochastically in their distance to the SMBH, which gives them a chance to momentarily harden (soften) further out (in) in the NSC (consider again Heggie’s law, Heggie 1975a). When it comes to the limits, e.g. the maximum separation we can expect at a given distance to the SMBH (which essentially is the purpose of Rose et al. 2020), migration, or these stochastic oscillations, become very important because it significantly changes the correlation between inner and outer orbits as a function of time. Therefore, finding a wider binary than expected does not imply a recent dynamical formation scenario for the binary, but nor does it exclude it.

Panamarev et al. (2019) ran an N-body simulation of the NSC for ~5 Gyr with 10⁶ particles and an initial binary fraction of 5%. Their initial binaries differ from ours in a few ways: They consider stellar masses of 0.08–100 M_⊙ and binary separations between 0.005 and 50 au, and more generally they cover the inner 10 pc of the NSC. Qualitative trends are the same: Wide binaries become less common with time, as they are either destroyed or become harder; the binary fraction has a radial dependence (e.g. more binaries remain intact further out); and binaries ≳1 au become very rare. In particular, they found that initially wide (soft) binaries, apart from evaporating, also become harder (see their Fig. 10). What is surprising, in comparison to our results, is the large number of hard binaries they obtained. In absolute terms, there is roughly a factor of three more hard binaries than in their initial population, with a maximum around a_bin ~ 2 R_⊙ (dominated by low-mass MS binaries). Overall, they find that about half of the initial binary population is destroyed after ~5 Gyr (30–40% are therefore hardened). Although this is dependent on their location in the NSC, this seems to be the case in the inner 1−0.1 parsec (see their Figs. 9 and 12). In contrast, we find a significantly stronger depletion of binaries along with a more pronounced dependence on distance from the SMBH over time, with roughly 80% (60%) originating from 0.1 pc (0.3 pc) being destroyed after ~5 Gyr (see e.g. Fig. 4). There are, however, a few possible explanations for these discrepancies: (i) The binary evolution code BSE (Hurley et al. 2002) used in their code has very efficient tidal circularisation for fully convective stars. (ii) Their paper lacks sufficient detail to determine precisely how binaries are treated, but it appears that the overall encounter rate is likely underestimated due to each N-body particle representing 65 stars. This would lead to an underestimation of binary destruction (e.g. through evaporation) and therefore alter the properties of the surviving binary population. (iii) Other relaxation processes not considered in our simulations could change the statistics of our outcomes (see a discussion in Section 5.3).

5.3 Resonant relaxation

Certain dynamical processes that were not considered in our simulations could have an affect on our results. Non-spherical fluctuations in the overall stellar distribution lead to torques acting on the stellar orbits; these torques drive diffusion in both the direction and the magnitude of the orbital angular momentum of the stellar orbits around the SMBH (Kocsis & Tremaine 2015; Hamers et al. 2018; Tep et al. 2021; Panamarev & Kocsis 2022). The random re-orientation of the direction of the angular momentum is known as vector resonant relaxation (VRR), while changes to the magnitude occur on longer timescales and is known as scalar resonant relaxation (SRR). The main effect of SRR is that the overall eccentricity distribution of the binaries’ outer orbits is driven to be thermal; as we already assume a thermal distribution, this would simply add an ever randomising effect on their outer eccentricities. This is most likely not going to affect the statistics of our results much. ZLK is operating on shorter timescales for virtually all distances considered (including migrated binaries) and there is only a weak dependence on outer orbit eccentricities when it comes to evaporation (Rose et al. 2020).

The effects of VRR are likely more relevant; a constant reorientation of the outer inclination (e.g. direction of angular momentum) has been shown to increase the number of mergers via the ZLK mechanism (the mutual inclination of a binary can sporadically reach ~90 degrees where the ZLK oscillations are the strongest Hamers et al. 2018). At the same time, Dodici et al. (2025) show that the combination of VRR, ZLK, flybys via the impulse approximation and strong tides can be very effective in circularising binaries, making them virtually stable against any continued ZLK as well as evaporation. It is important to note that this primarily affects our initially softer binaries since the harder binaries do not undergo ZLK to begin with. In other words, the binaries initially near the H/SB are unlikely to be affected significantly. How important these processes are therefore depends on the initial distribution one considers. Overall, our conclusion is that VRR likely increases the mergers and the amount of binaries that remain intact, while the number of evaporating binaries decreases. A comprehensive study of the complex interplay is left for a later investigation.

6 Conclusions

In this work, we have followed the dynamical evolution of 5736 lower-mass (≤2 M_⊙) binaries initially at 0.1 and 0.3 pc from the SMBH in the NSC. The dynamical evolution entails close encounters with field stars in the NSC, ZLK oscillations induced by the SMBH, tidal dissipation, and stellar evolution by interpolating pre-computed tracks. We find that binary migration caused by three-body encounters plays a vital role in determining the outcome of binaries. Inward migration leads to a certain death in the form of mergers and evaporation, while outward migration serves as a way to retain binaries within the NSC. We paid extra attention to the stellar collisions and binary mergers in these simulations because of their likely connection to the G-objects in the inner parts of the NSC. Although we find that binaries, as well as their merger products, can end up on orbits confined to this region – which suggests that the S-cluster and the G-objects can originate further away from the SMBH – the time it takes for the binaries to end up there is too long (0.3–1.2 Gyr) compared to the apparent ages of these objects (~10 Myr). Instead, we suggest that the S-cluster may be contaminated by low-mass BSSs that are related to earlier SFEs in the NSC (Schödel et al. 2020).

We also find that in a majority (≲80%) of cases where a surrounding star collides with one of the binary stars, it leads to a subsequent merger with the remaining binary star, resulting in a 3BPU. These are relatively common within the first 1–2 Gyr but stagnate later on. Based on a simple interpolation of the core-helium mass fraction of the 3BPU products, these objects could appear up to 10 Gyr younger than their true age. Finally, with a SFE possibly having occurred ~2 Gyr ago (Lu et al. 2013; Schödel et al. 2020), we would expect to see both recent mergers (e.g. in the form of G-like objects) and more evolved merger products (e.g. in the form of Thorsbro et al. 2023, BSS-like candidates) associated with this SFE scattered throughout the NSC as a consequence of binary migration.

Data availability

The data from this article will be shared upon reasonable request.

Acknowledgements

We would like to thank the anonymous referee for their comments and suggestions that greatly improved this paper. We give special thanks to Abbas Askar and Alessandra Mastrobuono-Battisti for their comments and valuable suggestions, which also helped improve this paper. The computations and data handling were enabled by resources provided by LUNARC, The Centre for Scientific and Technical Computing at Lund University. AAT acknowledges support from the Horizon Europe research and innovation programs under the Marie Skłodowska-Curie grant agreement no. 101103134. This work made use of Astropy: a community-developed core Python package and an ecosystem of tools and resources for astronomy (Astropy Collaboration 2013, 2018, 2022, http://www.astropy.org); Scipy: (Virtanen et al. 2020, https://scipy.org), Numpy: (Harris et al. 2020, https://numpy.org/); Matplotlib: (Hunter 2007, https://matplotlib.org) and Pandas: (Wes McKinney 2010, https://pandas.pydata.org).

References

Appendix A Accept-reject method

Mastrobuono-Battisti et al. (2021) constructed a method to sample encounters that lead to collisions between two single stars. We modify this method to handle interactions between a binary system and a single star. Binaries predominantly evolve through scatters, so the cross-sections considered for these three-body encounters are not confined to the physical sizes of the stars themselves. We first define the effective two-body periapsis distance r_p,2 as the closest approach between the tertiary and the centre of mass of the binary, assuming that the latter is a point mass. We then consider encounters out to a maximum distance of $r_{p,2}^{\max }=𝒞a\left(1+e\right)$, where a is the binary separation, e is the eccentricity, and 𝒞 is a scaling parameter. Generally, the farther away the periapsis approach of the tertiary is from the binary, the less the binary will be affected. The scaling parameter helps reduce encounters that have no effect on the binary. We define 𝒞 empirically as a function of a binary’s hardness; out to how many apoapsis do we expect encounters to induce changes >1% to the binary’s separation. Therefore, it may be thought of as the parameter that maximises the effective two-body periapsis distance (r_p,2) between the tertiary and the binary while keeping the number of interactions at a minimum. The empirical values and details on how they are obtained can be found in A.2.

For an encounter with total binary mass M_1,2, tertiary mass M₃, relative velocity at infinity v_∞, and impact parameter b, one can show that $$r_{p,2}=-\frac{G\left(M_{1,2}+M_3\right)}{2v_{\infty }^2}\left[1-\sqrt{1+{\left(\frac{2v_{\infty }^2}{G\left(M_{1,2}+M_3\right)}\right)}^2{b}^2}\right],$$(A.1)

where G is the gravitational constant.

With empirically determined values for 𝒞, we construct an accept-reject method for binary-tertiary encounters based on a modification of the approach in Mastrobuono-Battisti et al. (2021). First we sample a “velocity at infinity” from a Maxwellian distribution with a velocity dispersion according to eqn. (2), which at 0.1 pc gives σ ~ 250 km s⁻¹. We then define the encounter rate $R_{ij}^v$ (see eqn. A.2 and Mastrobuono-Battisti et al. 2021, eqn. 13) where n_j is the number density, M_i is the binary mass, M_j is the tertiary mass and R_ij = Ca(1 + e) is the maximised effective two-body periapsis (e.g. $r_{p,2}^{\max }$). We may choose a time step ∆t such that there is an interaction when $X<R_{ij}^v\Delta t$, where X is a randomly chosen number between 0 and 1. The associated ‘average’ timescale for an encounter is given by eqn. (A.3); where we have replaced v∞ with the velocity dispersion σ(r_•): $$R_{ij}^v\left(v_{\infty }\right)=n_i\pi R_{ij}^2{v}_{\infty }\left[1+\frac{2G\left(M_i+M_j\right)}{v_{\infty }^2{R}_{ij}}\right].$$(A.2)

Once we accepted an interaction, we could calculate the interaction radius (see eqn. (A.4)), where $R_{\text{rand}}^2=X/\left(\pi n{}_j{v}_{\infty }\Delta t\right)$ and $R_{\text{esc}}=2G\left(M_i+M_j\right)/v_{\infty }^2$. The interaction radius corresponds to the collision radius in Mastrobuono-Battisti et al. (eqn. 15 2021): $${\overline{t}}_{\text{enc}}=\frac{1}{2\pi }{\left(4\sqrt{\pi \sigma n{}_j{R}_{ij}^2\left(1+\frac{M_{ij}}{2{\sigma }^2{R}_{ij}}\right)}\right)}^{-1}.$$(A.3)

Using R_col and conservation of energy and angular momentum, we can infer the initial impact parameter of the tertiary using eqn. (A.5) $$R_{\text{col}}=\frac{R_{\text{esc}}}{2}\left[\sqrt{1+4{\left(\frac{R_{\text{rand}}}{R_{\text{esc}}}\right)}^2-1}\right]$$(A.4) $$b^2=R_{\text{col}}^2\left[1+\frac{2G\left(M_i+M_j\right)}{v_{\infty }^2{R}_{\text{col}}}\right].$$(A.5)

A.1 Setting up a tertiary

After having sampled a velocity at infinity v_∞ and an impact parameter b (see Appendix A) for the tertiary, we can determine the position and velocity of the tertiary. Although v∞ and the impact parameter are sampled at infinity, these quantities need to be extrapolated to a finite distance from the binary; the binary should also be in a random position in its orbit at the time of closest approach. As such, we placed the tertiary at a distance D where, at a constant velocity v_∞, it takes 10 binary orbital periods to reach the binary. The extrapolation is achieved using conservation of energy and angular momentum in which we treat the binary centre of mass as a point mass. Denoting the extrapolated velocity and impact parameter by v_D and b_D respectively, we obtained $$v_{\text{D}}=\sqrt{v_{\infty }^2+\frac{2G{M}_{ij}}{\sqrt{D^2+b_{\text{D}}^2}}};b_{\text{D}}=\frac{b{v}_{\infty }}{v_{\text{D}}}.$$(A.6)

With v_D and b_D, we constructed the tertiary’s initial 6D phase-space coordinates (position and velocity vectors). We begin by drawing a value for an angle ϕ ∈ [0, 2π] and z ∈ [−1, 1], both from uniform distributions. We then define $R=\sqrt{1--z^2}$ from which it follows that x = R sin ϕ and y = R cos ϕ. The initial position r without any offset is simply given by the three coordinates x, y, z multiplied by the distance D. To set up the offset vector b we first make a coordinate transformation and define $\widehat{\mathbf{\text{i}}}\text{=}\frac{\widehat{\mathbf{\text{x}}}\times \widehat{\mathbf{\text{r}}}}{\left|\widehat{\mathbf{\text{x}}}\times \widehat{\mathbf{\text{r}}}\right|}$¹⁶. We then obtained the impact parameter vector $$\mathbf{\text{b}}=\left(\widehat{\mathbf{\text{i}}}\text{\hspace{0.17em}sin\hspace{0.17em}}{\varphi }_2+\widehat{\mathbf{\text{j}}}\text{\hspace{0.17em}cos\hspace{0.17em}}{\varphi }_2\right)b_{\text{D}},$$(A.7)

where $\widehat{\mathbf{\text{j}}}=\widehat{\mathbf{\text{i}}}\times \widehat{\mathbf{\text{r}}}$ and ϕ₂ is drawn from a uniform distribution ϕ₂ ∈ [0, 2π). The phase-space coordinates with the offset are given by $$\{\begin{array}{cc}\mathbf{\text{x}}& D\mathbf{\text{r\hspace{0.17em}+\hspace{0.17em}b}}\\ \text{v}& -v_{\text{D}}\widehat{\mathbf{\text{r}}}\end{array}.$$(A.8)

Lastly, the tertiary’s mass is drawn from a two-part power law initial mass function (IMF) with α = −1.3 for masses below 0.5 M_⊙ and α = −2.3 for masses above (Kroupa 2001). We restrict the tertiary mass to 0.1 < m₃/M_⊙ < 1.0, where the upper constraint is imposed from the turn-off mass (~1 M_⊙) and all stars considered are on the MS throughout the simulation. For simplicity we adopt a fixed upper limit for m₃ which causes us to neglect a smaller number of encounters in the early life of the binary in which the tertiary could be more massive than the primary.

Fig. A.1 Left: Scaling parameter as a function of hardness. The edges are open to illustrate that any binary with, for example h > 10 has 𝒞 = 6. Right: Scatter plots for the fractional change in separation of binaries (top, panel b) and eccentricity (bottom, panel c), against the distance between the tertiary and binary centre-of-mass in units of the binary apoapsis. Only encounters with with a hardness of 0.01 <h≤ 0.10 are considered in this Figure. Blue (orange) dots correspond to encounters that change the binary separation with ≥ 1% (< 1%), where the black dotted lines show the 1% boundary and the black dash-dotted line the chosen value for the scaling parameter 6.

A.2 The scaling parameter C

In order to reduce the computational cost for these simulations, we have introduced a scaling parameter 𝒞 that maximises the effective two-body periapsis distance (r_p,2) between the tertiary and the binary while keeping the number of interactions at a minimum. By matching the inferred pericentric passage in eqn. (A.1) to the induced changes in a binary’s semi-major axis and eccentricity we can empirically determine values of 𝒞 for different hardness ratios. In total, we ran 10⁶ individual interactions for this empirical test; each encounter samples a tertiary and a binary in the same way as for the real simulations, see Section 3.2.1.

The left-hand side panel in Figure A.1 shows our chosen values as a function of the logarithm of the hardness ratio of the binary. The values of 𝒞 for the smallest and biggest hardness ratios are not closed (e.g. h = 10² gives the same scaling parameter as h = 10⁴). The number of encounters that fall in these brackets is virtually zero and is expected to have no effects on the main simulations. One of the brackets (0.01 <h≤ 0.10) is also highlighted; the two panels on the right show the out-comes of each encounter for these hardness ratios. In the top (bottom) panel, we plot |Δa/a| (|∆e/e|) against r_p,2/a(1 + e) ≡ 𝒞. Distant encounters are more likely to cause significant changes in a binary’s eccentricity than in its separation. This is mainly due to a larger spread, making it more difficult to justify specific values for 𝒞. For example, at 𝒞 = 8 (chosen value for the bracket 0.01 <h≤ 0.10), there are no encounters that induce changes 0.01 < |Δa|/a, but there are many encounters that induce changes 0.01 < |Δe|/e. We make the assumption that it is enough to consider the changes in separation for the scaling parameter. This means that we are neglecting multiple strong encounters with respect to changes in eccentricity. However, as we have roughly equally many encounters that decrease the eccentricity as there are encounters that increase it, we argue it is a reasonable assumption to make. We show this in Figure A.2; in the top panel we have the distributions of encounters (> 𝒞 = 8 for 0.01 <h≤ 0.10) that increase the eccentricity (black curve) and encounters that decrease it (red curve), normalised to the total number of encounters. In the bottom panel, we show the difference between the two distributions. On average, there is an increase of 0.0005 associated with each encounter. We are therefore underestimating the rate at which the eccentricity increases, but only by a very small amount.

Fig. A.2 Top: Distribution of the fractional effect on a binary’s eccentricity for encounters that increase (black line) and decrease (red line) it. Bottom: Residuals for the top panel. The dotted black line shows the average (0.0005) from the residuals.

Appendix B Secular processes

The Hamiltonian of a hierarchical triple can be decomposed into the inner orbit (e.g. the binary), the outer orbit (e.g. the binary centre of mass around the SMBH) and a coupling term describing the weak interaction between the orbits. The corresponding Hamiltonian is given by eqn. (1) in Naoz et al. (2013), and is defined to be the negative of the total energy, following the convention of Harrington (1968). Naoz et al. (2013) adopt the canonical variables known as Delaunay’s elements; these coordinates are chosen to be the mean anomalies l₁ and l₂, the longitudes of ascending nodes h₁ and h₂ and the argument of periapsis g₁ and g₂. Subscripts 1 and 2 denote the inner and outer orbits, respectively. For the full set of differential equations, we refer the reader to Naoz et al. (2013, see their Section 3 and Appendix A). We choose to integrate the differential equations of the Keplerian orbital elements (Naoz et al. 2013, see Eqns. A29–A35). We explicitly state the expression for the evolution of the semi-major axis in eqn. (B.1) which is not given in the original paper. This expression can be obtained by combining the conjugate momenta of the mean anomalies (L₁) and the arguments of periapsis (G₁), see eqns. (3) and (4) in Naoz et al. (2013), respectively. $${\dot{a}}_1=\frac{\sqrt{\left(m_1+m_2\right)a_1}\left({\dot{G}}_1+L_1\frac{e_1}{\sqrt{1-e_1^2}}{\dot{e}}_1\right)}{m_1{m}_2\sqrt{1-e_1^2}}.$$(B.1)

Here m₁, m₂ are the stellar masses of the binary stars, and a₁, e₁ are the separation and eccentricity of the inner orbit, respectively.

Appendix C Power-law fits

Appendix D Tidal disruption events

When binaries evaporate, the now single stars typically find themselves on eccentric bound orbits around the SMBH. There is then the possibility that they venture close enough to the SMBH to be tidally disrupted. A tidal disruption event (TDE) requires the periapsis to be smaller than the star’s tidal radius, and its angular momentum to be smaller than the critical angular momentum associated with the loss cone. Figure D.1 shows the ratio between the periapsis and the tidal radius along the x-axis and the ratio between the star’s angular momentum and the critical angular momentum along the y-axis. We find a total of five possible TDEs from these simulations.

Fig. D.1 Ratio between the periapsis and tidal radius versus the ratio between the specific angular momentum and the critical specific angular momentum for the loss-cone. Blue points correspond to the primary of the binary while orange points correspond to the secondary. The dotted black lines shows where respective ratio is unity – values within the box are inside the loss-cone.

Figure D.2 shows the stellar orbits around the SMBH (origin) in the left-hand side panel. The right-hand side panel shows the penetration parameter of a TDE $\beta \equiv \frac{r_t}{r_p}$ where β − 1 ≲ 1 corresponds to non-violent events that may lead to partial disruption of the stars (e.g. Phinney 1989; Guillochon & Ramirez-Ruiz 2013). While the number of TDEs is far too limited to infer any statistical properties of these events, we can say that the disruption of binaries will lead to TDEs and that they will most likely be non-violent.

Fig. D.2 Left: Orbits for the stars found inside the loss cone in Figure D.1 with the SMBH at the origin. Right: Scatter between the stars’ mass and the penetration parameter $\beta \equiv \frac{r_t}{r_p}$. The orbits and points share the same colour-coding as Figure D.1.

Appendix E Inner arcsec orbits

Figure E.1 shows a 2D projection of all merger products’ (see Figure 15) orbits with a semi-major axis smaller than 0.04 pc which corresponds to the inner arcsec of the GC (colour shows the time at which they merge); this is also more or less the region where the S-stars reside. We find that nine binaries merge with orbits in this region and migrate there within a few hundred megayears to a bit more than 1 Gyr (mean time is 660 Myr). All of the mergers within the inner arcsec originate from 0.1 pc which roughly corresponds to 1% of the mergers from this distance. In total there are 14 and 3 binaries originating at 0.1 pc and 0.3 pc, respectively, that at some point in their lives have an orbit with a_outer < 0.04 pc; those not included in the 9 mergers above, merge at a slightly greater distance. None of the additional binaries migrate into the region quicker than 300 Myr.

The bottom panel of the figure shows the cumulative distribution of eccentricities. The best power-law fit (solid black line) is η ≈ 0.08, which is surprisingly uniform in contrast to the eccentric orbits of the S-stars (Gillessen et al. 2009a,b; Sabha, N. et al. 2012) and the G-objects (Ciurlo et al. 2020). In the same panel, the dashed black line shows a super-thermal distribution for all binary merger orbits with a periapsis within the inner arcsec, which roughly corresponds to 60% of the mergers originating at 0.1 pc. The uniform distribution may therefore, or at least partially, be caused by a low number of orbits.

Fig. E.1 Top: Projection in 2D of the orbits of the nine merger products. The orbits are colour-coded as a function of the time of the merger. The black cross in the origin corresponds to the SMBH, and the dotted line marks the inner arcsec. Bottom: Cumulative distribution of the merger products’ eccentricities around the SMBH (similar to Figure 5). The solid black line corresponds to the merger orbits seen in the top panel, while the dashed black line corresponds to all mergers with a periapsis a(1 − e) ≤ 0.04 pc.

All Tables

All Figures

	Fig. 1 Initial (final) distribution of hardness ratio for binaries originating at 0.1 and 0.3 pc in solid (dashed) black and red lines, respectively. The vertical dotted line marks the H/SB. Each bin has been normalised to the total sum.
In the text

Fig. 2 Flow chart of our routine described in Section 3.2.3. The central sketch shows the overall setup of the three-body encounters; a binary (m1, m2) on an orbit at a distance of d• around the SMBH (MSMBH) interacting with a tertiary (m3). To the right of this sketch we show the possible outcomes from a three-body interaction. Group C contains interactions in which the binary remains intact and bound to the SMBH and as a consequence will experience continued evolution. Conversely, the evolution is terminated by outcomes in group T. To the left of the sketch we illustrate how a binary is evolved. A timestep, ∆t, is set by the shortest timescale of ZLK (quadrupole approximation), encounters, and circularisation (tides). The secular process with the shortest timescale min[τquad, τcirc] is then run for ∆t. At every step 0.1 ∆t, a three-body encounter is sampled with the accept-reject method described in Appendix A. If accepted, the three-body encounter is carried out with tsunami. If rejected, we proceed to the next 0.1 ∆t and sample again. After a time, ∆t, has passed or after a three-body encounter, a new ∆t is calculated. In this example, the encounter occurs after ∆t1 + 0.1 ∆t2 + tenc. This process is then continued until an outcome in group T occurs or until the simulation time exceeds 10 Gyr.

In the text

Fig. 3 Left: inner semi-major axis distribution of the initial population (black) and the corresponding final outcomes (colours). Yellow lines refer to binary mergers, green lines to binaries that evaporate, magenta lines to binaries that are ejected (unbound from the SMBH), red lines to intact binaries (i.e. binaries that survive to 10 Gyr), and blue lines to evolved binaries (stellar evolution of the primary occurs before binary destruction and by 10 Gyr). In the cases of evaporation and mergers, the properties shown in the figure refer to the binary just before either destruction mechanism takes place. Ejection corresponds to the region the binaries were ejected from and not where they end up. The top row corresponds to binaries starting with a• = 0.1 pc from the SMBH, while the bottom row corresponds to a• = 0.3 pc. Right: outer semi-major axis distributions. The vertical black lines represent the initial populations, and the boxes of the final outcomes are coloured in the same way as in the left panel.

In the text

	Fig. 4 Cumulative distribution of outcomes as a function of time. The outcomes are colour-coded the exact same way as in Figure 3. The solid (dashed) lines refer to binaries initially 0.1 pc (0.3 pc) away from the SMBH.
In the text

Fig. 5 Cumulative distribution of the eccentricities for the inner (panel a) and outer (panel b) orbit, colour-coded as in the previous figures (outcome). The dotted grey lines follow specific power-law distribution eη where η = −1 correspond to mostly circular orbits, η = 0 to a uniform distribution and η = 1 to a thermal distribution (e.g. mostly eccentric orbits). As in Figure 3, the coloured distributions consider the final state of a binary where it can still be considered a binary (and not destroyed). Best-fit parameters can be found in Table C.1.

In the text

Fig. 6 Panel a: distribution of the differences between changes in eccentricity due to ZLK and due to repeated scatters (only looking at binaries that eventually merge). The black line corresponds to binaries initially 0.1 pc away from the SMBH and the red line to 0.3 pc. Panel b: scatter of these two terms. One point corresponds to one binary. The zero-lines are shown as dotted black lines, and the shaded areas correspond to the ±0.05 region around the zero-lines. The two diagonal lines show x + y = 0 (lower) and x + y = 1 (upper).

In the text

Fig. 7 Four-panel plot showing the separation versus 1 − e of the inner orbit. Based on the zero-lines in Figure 6. The left-hand side column corresponds to the vertical zero-line (|ΣΔeZLK| ≤ 0.05), and the right-hand side column corresponds to the horizontal zero-line (|ΣΔeScatterl ≤ 0.05). The top row shows the initial distribution at t = 0, while the bottom row shows the final distribution at t = tf before the binaries merge.

In the text

	Fig. 8 Time of merger distribution for the binaries in Figure 7. Panel a, as before, refers to the vertical zero-line (effectively no ZLK) in Figure 6, while panel b refers to the horizontal zero-line (effectively no scatters). The black (red) line is for binaries initially at 0.1 pc (0.3 pc) away from the SMBH.
In the text

Fig. 9 Contours of the inner semi-major axis versus 1 − e for various outcomes at 0.1 pc from the SMBH. Upper row, from the left, initial distribution in black (a), merged binaries in orange (b), and evaporated binaries in green (c). Bottom row, from the left, the intact population (after 10 Gyr) in red (d), the evolved population in blue (e), and ejected binaries in magenta (f). For panels (b), (c), and (d), the distribution refers to the binary’s properties the moment just before the final outcome.

In the text

	Fig. 10 Same as Figure 9 but at 0.3 pc.
In the text

Fig. 11 Properties for collisions between a tertiary and one of the binary stars. Panel a shows the relation between the tertiary’s velocity at infinity and its impact parameter. The points are colour-coded based on their initial position from the SMBH. Panel b shows the impact parameter versus the binary’s semi-major axis before the collision. To help guide the eye, the dotted black line corresponds to ainner = b.

In the text

Fig. 12 Average number of collisions per binary as a function of time (solid lines). The dotted lines refer to the expected number of collisions a single 1 M⊙ star would have at 0.1 pc (black) and 0.3 pc (red). Because the binaries tend to migrate, the dashed lines show the same expectation as the dotted line, but where the distance is not fixed to the starting positions and instead follows the average outer orbit of the binaries. The pale dashed lines follow the apoapsis, while the non-transparent dashed lines follow the semi-major axis.

In the text

	Fig. 13 Mean outer semi-major axis for all remaining binaries as a function of time. Black lines show binaries starting at 0.1 pc from the SMBH; red lines show binaries starting at 0.3 pc. Each point shows the time and location of a collision.
In the text

	Fig. 14 Time of final outcome versus the time of collision for binaries that have suffered at least one collision. The dotted black line shows where tlife = tcoll; points on this line meet their final outcome as a direct consequence of the collision. The points are colour-coded based on the outcome similarly to Fig. 3.
In the text

Fig. 15 Panel a: normalised histograms for the total expected lifetime of the mergers. The dotted black line marks the age of the NSC at 10 Gyr. Panel b: post-merger semi-major axis around the SMBH. The dotted vertical lines mark the initial separations. The solid vertical line marks 1 arcsecond. Panel c: total mass after the merger. For all panels, binaries originating at 0.1 pc (0.3 pc) are shown in black (red).

In the text

Fig. 16 Panel a: rate of 3BPUs as a function of time. The dotted lines correspond to the mean rates over a 10 Gyr period. Panel b: normalised distribution of the estimated apparent age for the 3BPU merger products. The apparent age is determined by mapping the Helium-core mass fraction of the primary at the time of the merger to the corresponding age a normal star with the mass of the merger product would be with an identical Helium-core mass fraction. Panel c: normalised distributions of the total masses of the 3BPU merger products. For all panels, binaries originating at 0.1 pc are shown in black, and those from 0.3 pc are in red.

In the text

Fig. A.1 Left: Scaling parameter as a function of hardness. The edges are open to illustrate that any binary with, for example h > 10 has 𝒞 = 6. Right: Scatter plots for the fractional change in separation of binaries (top, panel b) and eccentricity (bottom, panel c), against the distance between the tertiary and binary centre-of-mass in units of the binary apoapsis. Only encounters with with a hardness of 0.01 <h≤ 0.10 are considered in this Figure. Blue (orange) dots correspond to encounters that change the binary separation with ≥ 1% (< 1%), where the black dotted lines show the 1% boundary and the black dash-dotted line the chosen value for the scaling parameter 6.

In the text

	Fig. A.2 Top: Distribution of the fractional effect on a binary’s eccentricity for encounters that increase (black line) and decrease (red line) it. Bottom: Residuals for the top panel. The dotted black line shows the average (0.0005) from the residuals.
In the text

	Fig. D.1 Ratio between the periapsis and tidal radius versus the ratio between the specific angular momentum and the critical specific angular momentum for the loss-cone. Blue points correspond to the primary of the binary while orange points correspond to the secondary. The dotted black lines shows where respective ratio is unity – values within the box are inside the loss-cone.
In the text

	Fig. D.2 Left: Orbits for the stars found inside the loss cone in Figure D.1 with the SMBH at the origin. Right: Scatter between the stars’ mass and the penetration parameter $\beta \equiv \frac{r_t}{r_p}$. The orbits and points share the same colour-coding as Figure D.1.
In the text

Fig. E.1 Top: Projection in 2D of the orbits of the nine merger products. The orbits are colour-coded as a function of the time of the merger. The black cross in the origin corresponds to the SMBH, and the dotted line marks the inner arcsec. Bottom: Cumulative distribution of the merger products’ eccentricities around the SMBH (similar to Figure 5). The solid black line corresponds to the merger orbits seen in the top panel, while the dashed black line corresponds to all mergers with a periapsis a(1 − e) ≤ 0.04 pc.

In the text