EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 703 (November 2025) A&A, 703 (2025) A157 Full HTML
Open Access
Issue
A&A
Volume 703, November 2025
Article Number A157
Number of page(s) 12
Section Galactic structure, stellar clusters and populations
DOI https://doi.org/10.1051/0004-6361/202556753
Published online 13 November 2025

© The Authors 2025

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1 Introduction

Star clusters are often used as laboratories for studying star formation, stellar evolution, and the dynamics of galaxies. When focusing on the old (>2 Gyr) and relatively massive (>104 M⊙) globular clusters (GCs), these gravitationally bound systems have long been known to host multiple populations of stars with distinct chemical compositions and evolutionary histories (see, e.g. Carretta et al. 2010; Gratton et al. 2012; Bellini et al. 2015; Bastian & Lardo 2018; Martocchia et al. 2018; Libralato et al. 2019; Tiongco et al. 2019; Martocchia 2020; Sollima 2021; Vesperini et al. 2021, and references therein). Hereafter, stars from the primordial (or first) population are referred to as ‘P1’, and the enriched second population as ‘P2’. For instance, variations in the abundances of light elements point to complex formation scenarios of P1 and P2, including the accretion of enriched material onto existing stars or a direct formation of P2 stars out of this enriched material (e.g. Carretta et al. 2010; Bekki 2011). In some clusters, P1 and P2 are also characterised by different dynamical evolution, which leads to measurable differences in present-day kinematic properties, such as rotation or anisotropy in velocity dispersion (Cordoni et al. 2020b,a; Dalessandro et al. 2024; Leitinger et al. 2025; Cadelano et al. 2024). Since velocity anisotropy affects the relaxation time scales, mass segregation, and the evolution towards energy equipartition (see, e.g. Tiongco et al. 2016; Bianchini et al. 2017; Pavlík & Vesperini 2021, 2022a,b; Pavlík et al. 2024), it is expected that each stellar population might evolve differently (as shown, e.g. by Livernois et al. 2024).

Most models of GC formation predict that P2 stars are being preferentially born in their cores (e.g. D’Ercole et al. 2008; Lacchin et al. 2022; Yaghoobi et al. 2022a,b). Although the spatial distribution of these populations evolves due to relaxation processes or gas expulsion (e.g. Vesperini et al. 2013; Decressin et al. 2010), dynamically younger GCs are expected to retain their initial structural differences (e.g. Dalessandro et al. 2019, 2024, and references therein). However, a recent observational study by Leitinger et al. (2023) showed that the spatial distribution of P1 and P2 varies among the dynamically young Galactic GCs. Notably, while P2 stars in some GCs (e.g. NGC 5024 or NGC 6809) appear to be more centrally concentrated than P1 (as expected for GCs less than a few relaxation times old), in other clusters which have similar dynamical ages as these ones (e.g. NGC 3201 or NGC 6101), P2 stars primarily occupy the outer regions of the clusters and P1 is preferentially in the central regions. There are also several dynamically young GCs which are fully mixed, i.e. they show no preferred concentration of either P1 or P2 (e.g. NGC 288, NGC 4590, NGC 5053, NGC 5904, NGC 6205, NGC 7078, and NGC 7089). At first glance, these findings contradict most models of GC formation that rely on P2 stars being born in a centrally concentrated manner. A possible explanation for this variety of morphologies may, however, also be that the internal dynamical processes (which might play out differently in each GC) are responsible for an early population mixing or inversion of an initially P2-centrally-concentrated GC, to instead show P2 stars in the outer regions.

We investigate whether the interactions between binary and single stars may be the dynamical origin of population mixing and inversion in GCs. Our preliminary study (Pavlík et al. 2025) with several small and compact clusters supported this idea. However, it only focused on the idealised case with full tangential stellar velocity anisotropy, where any systematic radial motion early in the cluster evolution would be produced by binary–single star scattering in the cluster core. Here, we investigate a larger suite of more realistic initial conditions and also elucidate the effects theoretically. A similar dynamical scenario, where stars are heated by a black hole subsystem, has also been independently proposed in a recent study by Giersz et al. (2025).

Generally, the binary–single star scattering is a prevalent cause of stellar ejections from the GC core (see, e.g. Aarseth 1972; O’Leary et al. 2014; Heggie & Hut 2003) and is usually mostly visible when the GC goes through core collapse later in its dynamical evolution (see, e.g. Tanikawa et al. 2012, 2013; Fujii & Portegies Zwart 2014; Pavlík & Šubr 2018). However, suppose a GC contains massive stars initially (O or B-type of several tens of M⊙). In that case, they will most likely be paired in binaries (Kroupa 1995; Goodwin & Kroupa 2005) and positioned in the core very early on – either the GC is already formed mass segregated, as has been shown for young star clusters (Plunkett et al. 2018; Pavlík et al. 2019; Pavlík 2020) but also for older GCs (Zonoozi et al. 2024), or they quickly sink into the core due to mass segregation (e.g. Binney & Tremaine 2008). Some of those very massive stars would evolve into black holes (BHs) on a time scale of a few Myr (e.g. Hurley et al. 2000, 2002; Belczynski et al. 2008, 2010; Ziosi et al. 2014). If their supernova natal kick velocity is then adequately small (see, e.g. Lyne & Lorimer 1994; Podsiadlowski et al. 2004), they may be retained in the GC (Peuten et al. 2016; Banerjee 2017; Baumgardt & Sollima 2017; Pavlík et al. 2018). Hence, BH–BH binaries might be produced promptly in GCs, subsequently ejecting P2 stars from the core to the outer regions through binary–single star scattering.

Our paper is organised as follows: in Section 2, we present theoretical estimates for the binary–single star scattering, in Section 3, we show numerical simulations designed to test our theory, in Section 4, we focus on the evolution of our models, in Section 5, we compare our key findings with real observed GCs, in Section 6, we discuss our results and possible alternative scenarios, and we summarise our conclusions in Section 7.

2 Theoretical estimates

We are interested in understanding the degree to which the binary–single scattering can eject an initially centrally concentrated P2 stellar population to the outer regions of the GC. In this section, we will estimate the rate of encounters and the typical properties of such binary and single stars. There are three important time scales to consider: (1) the time scale for a given binary to encounter a random single star, likely leading to the ejection of the single star from the core; (2) the time scale on which a significant fraction of the stars within the core encounter a binary and are, therefore, removed from the core; and (3) the time scale on which stars from outside the core will sink into the core. As the core originally contains only P2 stars with the P1 stars outside of the core, an inversion of the stellar population (i.e. P2 ending up largely placed further out than P1) will occur if the second time scale is shorter than the third but also if (for some reason) the binary encounters cease after this time. If, alternatively, the third time scale is shorter than the second, and/or the binary-single scattering process continues, then instead the P1 and P2 populations will be mixed by binary–single scattering.

In a first approximation, we consider equal-mass single stars in a GC (with masses m3) and a binary star with components m1 and m2 (typically two BHs), such that (m1 + m2) ≫ m3. Then the gravitationally focused close encounter timescale is (see, e.g. Eq. (7.196) in Binney & Tremaine 2008) τenc≈32.5 Myr(105/pc3n)(M⊙M)(10 auRmin)(v∞10 km s−1),$\[\tau_{\mathrm{enc}} \approx 32.5 ~\mathrm{Myr}\left(\frac{10^5 / \mathrm{pc}^3}{n}\right)\left(\frac{M_{\odot}}{M}\right)\left(\frac{10 ~\mathrm{au}}{R_{\min }}\right)\left(\frac{v_{\infty}}{10 \mathrm{~km} \mathrm{~s}^{-1}}\right),\]$(1)

where n is the stellar density (104 to 105 pc−3 for the central regions of a typical GC), M is the total mass of the interacting bodies (i.e. M = m1 + m2 + m3), Rmin is the distance of the closest approach (where strong scattering typically occurs when Rmin ≈ abin, the binary semi-major axis), and v∞ is the stellar velocity at infinity (comparable to the velocity dispersion of the m3 stars). Thus, the question we aim to answer is, what are the possible parameters for the binary to have a sufficiently short timescale τenc to produce enough scatterings while the cluster is dynamically young?

The energy transfer scales with the binary properties as follows. The binary binding energy is Ebin=Gm1m22abin,$\[E_{\mathrm{bin}}=\frac{G m_1 m_2}{2 a_{\mathrm{bin}}},\]$(2)

where G is the Newtonian gravitational constant. The energy of a hard binary is Ebin >12m3v∞2$\[E_{\text {bin }}>\frac{1}{2} m_{3} v_{\infty}^{2}\]$ (Heggie 1975; Hut 1983; Heggie & Hut 2003) for the m3 perturber mass. Or in terms of the semi-major axis, the hard/soft limit is abin≈8.87 au(m1 m2m3 M⊙)(10 km s−1v∞)2,$\[a_{\mathrm{bin}} \approx 8.87 ~\mathrm{au}\left(\frac{m_1 ~m_2}{m_3 ~M_{\odot}}\right)\left(\frac{10 \mathrm{~km} \mathrm{~s}^{-1}}{v_{\infty}}\right)^2,\]$(3)

where we used G ≈ 887 au (km s−1)2 M⊙−1$\[M_{\odot}^{-1}\]$. The change in binding energy after an encounter with an incoming single star is ΔEbin=12(m1+m2)vbin2+12m3v32,$\[\Delta E_{\mathrm{bin}}=\frac{1}{2}\left(m_1+m_2\right) v_{\mathrm{bin}}^2+\frac{1}{2} m_3 v_3^2,\]$(4)

where vbin and v3 are the typical velocities per encounter between the binary’s centre of mass and a single star in the GC, respectively. Using the conservation of momentum in the binary’s centre-of-mass frame, i.e. (m1 + m2)vbin = m3v3, Eq. (4) becomes ΔEbin=12m3(m1+m2+m3)m1+m2v32.$\[\Delta E_{\mathrm{bin}}=\frac{1}{2} \frac{m_3\left(m_1+m_2+m_3\right)}{m_1+m_2} v_3^2.\]$(5)

The average change of binding energy during an encounter can be rewritten as a fraction of the binary binding energy, i.e. ΔEbin = αEbin, where α is a distribution function describing the strength of the encounters. The mean stellar speed per encounter then becomes v3=α Gm1m2(m1+m2)m3(m1+m2+m3)1abin.$\[v_3=\sqrt{\alpha ~G \frac{m_1 m_2\left(m_1+m_2\right)}{m_3\left(m_1+m_2+m_3\right)} \frac{1}{a_{\mathrm{bin}}}}.\]$(6)

Assuming, for instance, that the cluster follows a Plummer model (Plummer 1911), the stellar escape velocity from radius r is vesc(r)≈2 GMcl[r2+(rh1.3)2]−1/4,$\[v_{\mathrm{esc}}(r) \approx \sqrt{2 ~G M_{\mathrm{cl}}}\left[r^2+\left(\frac{r_{\mathrm{h}}}{1.3}\right)^2\right]^{-1 / 4},\]$(7)

where rh is the half-mass radius and Mcl is the cluster mass, i.e. Mcl = Nm3 + Nbin(m1 + m2) ≈ Nm3 (if the mass of the binary stars is negligible compared to the cluster mass). Thus, to eject single stars from the core region, i.e. v3 ≳ vesc(0), Eqs. (6) and (7) give the required binary star semi-major axis abin≲0.38 α rhm1m2(m1+m2)m3(m1+m2+m3)1Mcl.$\[a_{\mathrm{bin}} \lesssim 0.38 ~\alpha ~r_{\mathrm{h}} \frac{m_1 m_2\left(m_1+m_2\right)}{m_3\left(m_1+m_2+m_3\right)} \frac{1}{M_{\mathrm{cl}}}.\]$(8)

Using Eqs. (1) and (8) and taking Rmin = abin, the time scale for this binary to encounter a random single star (and get it ejected) is τenc≳1.26⋅103 Myrα(rhpc)2(m32 M⊙m1m2(m1+m2))(v∞10 km s−1).$\[\tau_{\mathrm{enc}} \gtrsim \frac{1.26 \cdot 10^3 ~\mathrm{Myr}}{\alpha}\left(\frac{r_{\mathrm{h}}}{\mathrm{pc}}\right)^2\left(\frac{m_3^2 ~M_{\odot}}{m_1 m_2\left(m_1+m_2\right)}\right)\left(\frac{v_{\infty}}{10 \mathrm{~km} \mathrm{~s}^{-1}}\right).\]$(9)

For an example of a BH–BH binary (m1 = m2 = 30 M⊙) and a single star (m3 = 1 M⊙) in a cluster with rh = 0.5 pc and v∞ = 10 km s−1, Eq. (9) yields τenc ≳ α−1 5.83 · 10−3 Myr. We also plot other estimates derived from Eqs. (8) and (9) in Fig. 1.

To achieve mixing (or radial inversion) of P1 and P2, we need to remove a significant fraction of P2 stars (i.e. NencP2 ~ NP2) from the core where they were born. Therefore, the time scale for a binary star to encounter and eject most/all of them is τej≡τencNencP2.$\[\tau_{\mathrm{ej}} \equiv \tau_{\mathrm{enc}} N_{\mathrm{encP} 2}.\]$(10)

However, a single binary can only interact with a fraction of nearby stars whose trajectories are within its cross-section. Thus, to allow most stars in the core to have the opportunity to interact with the binary over τej, their orbits must be sufficiently perturbed by relaxation. This implies that the time scales are τrc ≲ τej, where τrc is the core relaxation time (see, e.g. Heggie & Hut 2003), defined as τrc(t)=0.34 σc3(t)⟨mc(t)⟩ ρc(t) ln⁡ [0.4N(t)],$\[\tau_{\mathrm{rc}}(t)=\frac{0.34 ~\sigma_{\mathrm{c}}^3(t)}{\left\langle m_{\mathrm{c}}(t)\right\rangle ~\rho_{\mathrm{c}}(t) ~\ln~ [0.4 N(t)]},\]$(11)

where σc is the core velocity dispersion, ⟨mc⟩ is the mean stellar mass in the core, ρc is the core mass-density, and N(t) is the number of stars in the cluster at time t1.

To observe a mixed/inverted dynamically young GC, most of the binary–P2-star encounters need to occur early in the cluster evolution, so τej per binary must be much lower than the half-mass relaxation time of the cluster (e.g. Spitzer & Hart 1971) τrh(t)≈0.138ln⁡ [0.4N(t)]N(t) rh3(t)G m3.$\[\tau_{\mathrm{rh}}(t) \approx \frac{0.138}{\ln~ [0.4 N(t)]} \sqrt{\frac{N(t) ~r_{\mathrm{h}}^3(t)}{G ~m_3}}.\]$(12)

This is because even if we start from completely P2-centrally-concentrated initial conditions, the P1 stars outside of the core will fall inwards as the cluster relaxes. Thus, we should expect the number of binary–P2-star scatterings to reduce over time and the number of binary–P1-star scatterings to increase. Once the probability of scattering P1 and P2 becomes comparable, the binary would only be able to achieve mixing of the two populations because P1 will get scattered just as much as P2. Consequently, to achieve inversion, the majority of P2 ejections via binary–single scattering need to happen early on and then stop.

We must explore how to set these three timescales above to kick most of the P2 outwards very early. We note that in this section, we have assumed that the single m3-stars must escape from the GC completely. However, Leitinger et al. (2023) have shown that the radial distribution of P2 in some GCs (e.g. NGC 6101 and NGC 3201) reaches its maximum at the projected half-light radius2. Thus, mixing and partial radial inversion of P1 and P2 may also be driven by slightly less energetic encounters, and, for instance, the values of the semi-major axis in Fig. 1 and Eq. (8) are rather the lower estimates.

thumbnail Fig. 1

Estimated semi-major axis (top) and the encounter time scale (bottom) of a binary star that can eject single stars from the core in a cluster of a given mass and half-mass radius. In this plot, we used Eqs. (8) and (9), for clusters containing stars of equal masses (m3 = 1 M⊙), one binary star with m1 = m2 = 30 M⊙ at its centre, and v∞ = 10 km s−1. Time in the bottom panel is normalised by the half-mass relaxation time, Eq. (12).
thumbnail Fig. 2

Evolution of the half-mass relaxation time, τrh(t), calculated from Eq. (12) in terms of the initial half-mass relaxation time, τrh(0). The left-hand set of panels is for the 10k-* clusters, the right-hand set is for 50k-*. The models are separated by the primordial binary star numbers and masses (columns), the binary star initial semi-major axes (colours), and the initial cluster half-mass radius (line styles). Each line is averaged over all realisations of the corresponding model. (We note that the models 10k-5M* and 50k-30M* do not show any differences from the models without binaries, which we also integrated for comparison but do not display them here.

3 Numerical models

3.1 Initial conditions of the star clusters

We performed a grid of simulations to test the theoretical estimates above. We used clusters with N = 104 stars (to probe the parameter space with computational simplicity) and N = 5·104 stars (for a more realistically sized system). Some of these stars were paired into Nbin = 1, 3, 10 or 20 massive primordial binaries, and the remaining ones were single equal-mass stars (with m3 = 1 M⊙). The clusters were isolated, had a Plummer (1911) profile with the half-mass radius rh = 0.5, 1.0, 2.0 or 4.0 pc, isotropic velocity distribution, and were initialised in virial equilibrium. The initial conditions were generated with MCLUSTER (Küpper et al. 2011) and integrated with the N-body tree code PETAR (Wang et al. 2020a) in several random realisations each3. The simulation time corresponded to several τrh(0) – see also Fig. 2 with the relaxation times of the models for comparison. We note that in Section 5, we compare these models with significantly older GCs in terms of physical units of time (age of several Gyr, whereas our models range between 100 Myr and 2.5 Gyr), however, in terms of their dynamical age (i.e. the number of τrh), the observed systems are similar to our models.

In each model, we split the single m3-stars into two populations based on their total initial energy, where for the i-th star Etot,i=12mivi2−GMcl[ri2+(rh1.3)2]−1/2.$\[E_{\mathrm{tot}, i}=\frac{1}{2} m_i v_i^2-G M_{\mathrm{cl}}\left[r_i^2+\left(\frac{r_{\mathrm{h}}}{1.3}\right)^2\right]^{-1 / 2}.\]$(13)

We labelled as ‘P2’ the stars with total energy in the lowest quartile, and the rest were labelled ‘P1’4. This separation was based on the assumption that P2 stars are born more concentrated and would have their orbits initially more confined in the central regions. It is usually assumed that P2 would be formed 30–100 Myr after P1 if the cause of enrichment is winds from the asymptotic-giant-branch stars or from interacting binaries (e.g. Bastian & Lardo 2018). However, in other formation scenarios, where the enrichment is caused by very-/super-massive stars (e.g. Gieles et al. 2018) or stellar-wind mass-loss from stellar merger events (e.g. Wang et al. 2020b), P2 may form within a few Myr after P1. Simulating two populations formed at two different times is beyond the scope of this paper; nevertheless, even our simple distinction between P1 and P2 helps us reveal several effects of the binary–single star interactions and spatial redistribution of the stars.

3.2 Primordial binary stars

The binary components had m1 = m2 = 5 or 30 M⊙ depending on the model5. They were generated on circular orbits with semi-major axis abin = 5, 10, 20, or 50 au – this made them very hard6. Even the choice of the widest 50 au binary was substantially small to theoretically drive the radial redistribution of P1 and P2 stars (see Fig. 1).

The binaries were initially randomly placed inside the central region, but they quickly segregated to the core, from where they started ejecting stars. The mass segregation time scale may be estimated from (Spitzer & Hart 1971; Spitzer 1987; Binney & Tremaine 2008) τseg≈m3m1τrh,$\[\tau_{\mathrm{seg}} \approx \frac{m_3}{m_1} \tau_{\mathrm{rh}},\]$(14)

where we approximate the average stellar mass by m3. Plugging in the numbers from each model and using the relaxation time from Fig. 2, we may verify that τenc ≫ τseg, hence the binary does not have much time to interact with the m3-stars while sinking towards the GC core.

The names of the models in this paper follow the convention: [N] k−[Nbin x( if >1)][m1/M⊙] M r[rh/pc] a [abin /au]$\[[N] ~\mathrm{k}-\left[N_{\text {bin}} ~\mathrm{x}(\text { if }>1)\right]\left[m_1 / M_{\odot}\right] ~\mathrm{M~r}\left[r_{\mathrm{h}} / \mathrm{pc}\right] \text { a }\left[a_{\text {bin }} / \mathrm{au}\right]\]$

For instance a model with 104 stars, rh = 1 pc, one binary of (5+5) M⊙ and abin = 20 au is called 10k-5Mr1.0a20. We also use a wildcard ‘*’ to refer to a group of models – e.g. all models with 104 stars, three (30+30) M⊙ binaries and rh = 0.5 pc, regardless of the binary’s semi-major axis, are collectively called 10k-3x30Mr0.5*. We also note that we did not follow all possible combinations of the initial conditions due to the cost of the simulations and their predictive nature.

thumbnail Fig. 3

Evolution of the core relaxation time, τrc(t) in the units of τrh(t), calculated from Eq. (11). The horizontal axis shows time in terms of the initial half-mass relaxation time, τrh(0). The left-hand set of panels is for the 10k-* clusters, the right-hand set is for 50k-*. The models are separated by the primordial binary star numbers and masses (columns), the binary star initial semi-major axes (colours), and the initial cluster half-mass radius (line styles). Each line is averaged over all realisations of the corresponding model.

4 Evolution of star clusters and binary stars

First, we explore the GCs’ evolution, focusing on their 3D structure. The time evolution of the relaxation time is plotted in Fig. 2, calculated from Eq. (12). We note that although our models are isolated, we take into account that some stars might be escaping when calculating τrh. Thus, we iteratively determined the number of stars N (and the enclosed cluster mass, Mcl) within the cluster’s Jacobi radius rJ=xGal(Mcl3MGal)1/3,$\[r_{\mathrm{J}}=x_{\mathrm{Gal}}\left(\frac{M_{\mathrm{cl}}}{3 M_{\mathrm{Gal}}}\right)^{1 / 3},\]$(15)

with the assumed orbit of the cluster at xGal = 5 kpc, corresponding to the Galactic mass of MGal ≈ 5 · 1010 M⊙ (Faber & Gallagher 1979; Bland-Hawthorn & Gerhard 2016). Fig. 2 shows that the models with one (5+5) M⊙ or one (30+30) M⊙ binary (i.e. 10k-5M* and *-30M*) have near-constant τrh(t) and therefore near-constant rh(t). They evolve analogously to single-component models without binaries. In other words, the heating power of a single stellar-mass binary BH is either not high enough to cause a notable expansion of the cluster overall, or the m3-stars do not go through as many scattering encounters on the same time scale as in the models with higher binary counts.

Similarly, we can infer from Fig. 3 that the cores of the clusters with a single low-mass binary (i.e. 10k-5M* and 50k-30M* models) gradually proceed towards core collapse. On the other hand, the cores of the models with massive or multiple binary stars contract more rapidly initially. In the case of the models 50k-10x30M*, 50k-20x30M*, and to an extent also 10k-10x30M, Fig. 3 also shows an increase of τrc after the initial drop, corresponding to the expansion of the core due to binary heating.

Fig. 2 further shows that τrh(t) rises over time in the models with more massive binaries. This expansion is greater in the following scenarios:

  1. In the most compact clusters, because the relaxation processes are more rapid and the binary–single star interactions can occur more frequently in denser environments (compare the different line styles in each panel of Fig. 2).

  2. In the clusters with wider binaries, because the scattering cross-section of each binary is larger (e.g. compare the line bifurcation in the models 10k-10x30Mr0.5* in the fourth panel from the left of the figure).

  3. In the models with more binaries, as their scattering power increases (compare the same line styles across the panels of the figure).

Consequently, the present-day τrh of the GCs with multiple binaries does not correspond to the relaxation time scale at birth. This means that the dynamical age of such GCs might be under-estimated if only inferred from their present-day structure. For comparison, we present some of our results (e.g. in Figs. 4 and 5) in terms of the initial and present-day (i.e. end-of-simulation) τrh simultaneously. And later, when we compare our models with observations, we do so only in terms of the present-day τrh.

We must also emphasise that while the dynamical heating increases with the number of primordial binaries in the GC, it is not exactly proportional to it because these binaries also interact with each other and destroy themselves. We show this in Figs. 4 and 5, where we plot the number of massive binaries in our models over time. For instance, see the panels for 10k-30M* and 10k-3x30M* in Fig. 4 – the latter had three binaries initially, however, most of them are quickly destroyed and, on average, only one or two remain after 1 τrh, while the former started and finished the simulation with one binary. This binary–binary destruction is more prevalent when the binaries have larger semi-major axes – see that in most cases, the lighter-coloured lines (tighter binaries) are above the darker-coloured lines (wider binaries) in these figures. Temporary oscillations of Nbin in Figs. 4 and 5 are mainly due to the binary capturing a single star and entering a triple system (occasionally also a quadruple), and we do not count multiples in these figures.

Comparing the plots in Fig. 2 with the highest number of primordial binaries (i.e. 10k-10x30M* and 50k-20x30M*), we see that even with twice as many binaries, τrh(t) does not grow as much in the latter model. This is due to the smaller mass ratio Nbin(m1+m2)/Mcl in the 50k-20x30M* models than in the 10k-10x30M* models. Moreover, this ratio is further lowered as the number of binary BHs in the 50k-20x30M* models drops to a number between four and ten (see Fig. 5), comparable to the final binary count in 10k-10x30M* (see Fig. 4). A much more numerous BH subsystem or more massive BHs might be necessary to cause an expansion of 50k-20x30M* similar to the 10k-10x30M* model.

Moreover, Fig. 3 further shows that the ratio τrc/τrh in the 50k-20x30M* models is higher than in the 10k-10x30M* models. This signals (as discussed in Sect. 2) that the binaries in the 50k-20x30M* models cannot interact with as many core P2 stars as in the 10k-10x30M* models, and the population mixing in the former models would be less effective.

In Figs. 6 and 7, we show the distribution of α, introduced in Sect. 2 as the change of the binary binding energy per interaction with a single star. We also plot the theoretical lower limit for the interactions that should lead to the escape of the m3-stars. Comparing the models of different initial sizes (panels on the same lines in the figures), we see that as rh(0) of the system increases, the probability of strong interactions leading to escape drops. This is most visible in the models with a single massive binary star (i.e. 10k-30M* in Fig. 6 or 50k-30M* in Fig. 7). While even the wider models with (30+30) M⊙ binaries still produce a population of escaping stars, we may notice that the 10k-5Mr2.0* models with the least massive primordial binary are essentially unable to eject stars (see the top-right panel of Fig. 6), thus any population inversion or mixing in these models should happen at the lowest rate.

While the models with multiple binaries also produce a population of escapers, their distributions of α are more shifted below the limiting value for escape (compare the models in the same column in Figs. 6 or 7, e.g. 50k-*r0.5* in the left-hand panels of Fig. 7). This suggests that as the binary star number in the system increases, a higher number of the single m3-stars that pass through the core region are subject to low-energy binary–single interactions rather than a single encounter that would kick them out of the cluster. Consequently, we may expect that the binary–single interactions leading to the mixing (or inversion) of the radial distributions of P1 and P2 would predominantly happen in the more compact clusters with the more massive binary stars. Or, in other words, that the more compact clusters would mix more rapidly than the extended clusters of the same mass and binaries.

thumbnail Fig. 4

Number of massive binaries in each 10k-* model in time (averaged over each model’s realisations). The models are separated by their initial rh (columns), primordial binary systems (rows), and their semi-major axes (colours). In each plot, the top horizontal axis shows the evolutionary time in the initial half-mass relaxation times, and the lower horizontal axis is in multiples of τrh at the end of the simulation (to mimic the ‘present-day’ value). While we permit the dynamical formation of the low-mass (m3+m3) or mixed (m1,2+m3) binaries in the simulations, the numbers here only show the massive (m1+m2) binaries, including component switching. However, if the primordial binary goes through a state of an unstable triple, it is not counted towards Nbin (which is why we sometimes see the number of binaries decreasing and then increasing again).
thumbnail Fig. 5

Same as Fig. 4 but for the 50k-* models.
thumbnail Fig. 6

Probability density function (PDF) of the binding energy change, α, per interaction, in the 10k-* models. Binary–single and binary–binary interactions are counted together. For reference, we also plot with dashed vertical lines the theoretical values of α, for which a single m3-star should theoretically escape, as derived using Eq. (8). The figure is separated by each model’s initial half-mass radius (columns), number and masses of the primordial binaries (rows), and the binary semi-major axes (colours).
thumbnail Fig. 7

Same as Fig. 6 but for the 50k-* models.
thumbnail Fig. 8

Selected 10k-* models with the highest A+ at 4 τrh for each primordial binary set-up (as labelled in the top left corners). Top panels: normalised projected cumulative radial distributions of stars in two populations (P1 and P2 – distinguished by colours) in the clusters, calculated from Eq. (16). All realisations of the corresponding model are plotted. Columns correspond to different values of Rlim. Bottom panels: the ratio of P2 stars in several radial bins (similarly to Leitinger et al. 2023). Each model realisation is plotted with a black line. The parameter ARlim+$\[A_{R\text {lim}}^{+}\]$ (its maximum range and median with standard deviation) is also displayed for each region.

5 Comparison with observations

To evaluate the effect of population redistribution due to binary heating, we calculate the area parameter (introduced by Alessandrini et al. 2016) which shows the relative abundance of P1 and P2 from the projected normalised cumulative radial distributions in several radial regions: ARlim+=1Rh∫0Rlim[CP1(R)CP1(Rlim)−CP2(R)CP2(Rlim)]dR,$\[A_{R \text{lim} }^{+}=\frac{1}{R_{\mathrm{h}}} \int_0^{R_{\lim }}\left[\frac{C_{\mathrm{P} 1}(R)}{C_{\mathrm{P} 1}\left(R_{\lim }\right)}-\frac{C_{\mathrm{P} 2}(R)}{C_{\mathrm{P} 2}(R_{\lim })}\right] \mathrm{d} R,\]$(16)

where 𝒞(R) is the cumulative radial distribution of stars at a projected radius R for the specified population, and 𝒞(Rlim) ensures its normalisation to one at the limiting radius. Hence, the methodology is similar to Leitinger et al. (2023, their Figure 14 and Appendix A). For our analysis, we assume Rlim = 1 Rh (central region), 4 Rh (to compare with Leitinger et al. 2023, 2025), and also 10 Rh (for the whole cluster), where Rh is the projected half-mass radius.

In Figs. 8 and 9, we plot the radial distributions of P1 and P2 stars in selected models at 4 τrh (each with a different number of binaries) where we found the highest values of A+. All these simulations correspond to rh(0) = 0.5 pc, which is expected because the densest clusters are also the most dynamically active, and the binary–single star scattering there is more effective. In all models, P1 and P2 mix in the central regions up to the half-mass radius (i.e. A1+$\[A_{1}^{+}\]$ ≈ 0, albeit the values are slightly more negative in the models with fewer or less massive binaries). The cumulative distribution functions of the model 10k-10x30Mr0.5a50 in Fig. 8 show a higher level of mixing even up to 4 or 10 Rh, proving that a higher number of binary stars can mix the central and outer populations more than fewer binaries on the same dynamical time scale. This is further supported by the bottom panels of each set of plots in the same figure, which show that the depletion of P2 stars from the core and their migration to the outer regions is dominant in the models with multiple binaries.

The more populous 50k-* models plotted in Fig. 9 do not exhibit as much variation between realisations as the less populous models in Fig. 8 because the distributions are averaged over a larger number of stars in the clusters. On the other hand, the distributions of P2 stars in the bottom panels of Fig. 9 show that the mixing power of even 20 primordial binaries with 30 M⊙ components is not enough to cause significant migration of the P2 stars outwards.

Furthermore, in Figs. 10 and 11, we plot the time evolution of the parameter A+ in the 10k-* and 50k-* clusters, respectively. All models have highly negative initial conditions (P2 is centrally concentrated by definition) and then gradually evolve towards a mixed state (A+≈0). We see that the more concentrated clusters and those with more binaries mix slightly more rapidly; however, all clusters with a given number of primordial binaries reach similar A1+$\[A_{1}^{+}\]$ values at 1–2 τrh. On the other hand, the more compact clusters show higher A4+$\[A_{4}^{+}\]$ at any given time (compare the lines of the same colour across different panels of the figure). Similarly, the models with a lower binary count require significantly longer dynamical evolution to reach the same A4+$\[A_{4}^{+}\]$ values as the more compact ones (compare the differently coloured lines in the same panel)7. We note that our models do not showcase a clear population inversion (A1+$\[A_{1}^{+}\]$ ≫ 0) anywhere up to 4.5 τrh.

In the bottom panels of both figures, we plot several dynamically young GCs from the sample of Leitinger et al. (2023). Focusing only on the 10k-* models in Fig. 10, we see that those with ten primordial binaries and rh(0) = 0.5 pc fit several of the very mixed GCs within their 1-sigma uncertainties. The less-mixed GCs (e.g. NGC 5024, 5272, and 6809) are, in turn, well represented by the models with fewer primordial binaries or the more expanded clusters.

Comparing Figs. 10 and 11, the 50k-* models show lower spread of A+ than the less populous 10k-* models. These fluctuations and inconsistencies between realisations of the same model are caused by the specific dynamical evolution and binary–binary and binary–single star interactions in each cluster, which either lead to their escape or ionisation, temporarily reducing the binary heating engine in the core (as we discussed also in terms of the radial distribution of P2 stars). We also see that the mixing of P1 and P2 is not as effective in the more populous clusters, even in the most compact ones with 20 primordial binaries. On the other hand, the trend we described above (‘more binaries equal more mixing’) is still visible. This suggests that in a more populous cluster, either a more numerous stellar-mass BH–BH subsystem is necessary, or we need to increase the BH masses – to a degree, this is also suggested by the theoretical derivations we did in Sect. 2, specifically in Fig. 1.

While our results do not provide conclusive evidence for the cause of P1/P2 mixing in young GCs, they are highly suggestive that the dynamical heating by binaries plays a role, especially when it is set appropriately to the size of the system. It is, therefore, worth investigating this in subsequent studies.

thumbnail Fig. 9

Same as Fig. 8 but for the 50k-* models.

6 Discussion

We acknowledge that we have only explored numerically the domain of clusters with ten to fifty thousand stars, and we need to map those results to ‘real’ GCs containing up to a few million particles. The distributions of α (presented in Figs. 6 and 7) show that the results on binary – single star scatterings leading to stellar escapes in our models are comparable between the models, although one has five times more stars (e.g. compare the curves for 10k-30M* and 50k-30M*). This also refers back to Fig. 1, where we see that the relationship between the binary semi-major axis and the star cluster size and mass is self-similar. Therefore, this knowledge about the distribution of α from our less populous models can be applied to understanding what fraction of encounters would lead to the ejection from the core, the ejection from the cluster entirely, or placing the core stars in the halo in real GCs.

Despite our designing of star clusters where massive binary stars could efficiently eject the central population of P2 stars, our models do not show the radial inversion of P1 and P2 but only mixing. Nonetheless, we must note that out of the 13 dynamically young GCs in the sample of Leitinger et al. (2023, 2025), four have A4+$\[A_{4}^{+}\]$ < 0, two have A4+$\[A_{4}^{+}\]$ > 0, and seven have A4+$\[A_{4}^{+}\]$ ≈ 0. This indicates – although the sample is small – that high positive A+ values may be rare to achieve. We also note that the two clusters with A4+$\[A_{4}^{+}\]$ > 0 have the following radial distributions. NGC 3201 primarily contains P2 stars within the range up to 0.5 Rh and beyond 1.5 Rh, and P1 stars in the range from 0.5 Rh to 1.5 Rh. On the other hand, NGC 6101 shows a consistent central concentration of P1 stars throughout the GC and is, therefore, technically the only GC with a clear P1 central concentration. Meanwhile, NGC 3201 may be indicative of a GC that began with a P2 central concentration and underwent a semi-successful inversion of its populations. Furthermore, both of these GCs with high positive A4+$\[A_{4}^{+}\]$ are reportedly accreted from the Sequoia dwarf galaxy (Massari et al. 2019), i.e. the formation conditions might have determined their structural properties more than in the other GCs (see also Figure 15 in Leitinger et al. 2025, which connects the progenitor system of the GCs).

More to the point, while the binary heating engine in our models caused the outward migration of the central P2 stars, it also pushed all the sinking P1 stars to higher radii as the cluster relaxed. We introduced four evolutionary time scales in Sect. 2 – those are (1) the time scale for individual encounters between a given binary and a random single star, (2) the time scale to eject a significant amount of P2 through binary–single star encounters from the core, (3) the time scale to sufficiently mix up the stars in the core, such that all core stars have the opportunity to be ejected, which is related to the core relaxation time and (4) the time scale for the infall of P1, which is related to the relaxation time scale of the entire system. To achieve the radial inversion of P1 and P2, the ejection of P2 from the core must be more rapid than the infall of P1, and once most of the P2 stars are gone from the core, it has to stop before P1 builds up in the core. If, however, the binary engine continues to scatter, or the infall of P1 happens on a shorter time scale than the ejection of the core stars, the only possible outcome is P1/P2 mixing. Although we saw in our models that any high initial binary count has been reduced over time, it was never fully stopped. The kicking of stellar-mass BHs out is rather stochastic and cannot be timed or predicted from the initial conditions. Consequently, all our models mixed P1 and P2 in the long term.

An alternative mechanism to quench the binary–single star scatterings is to have a single binary system composed of two intermediate-mass black holes (IMBHs). Let us assume one of our modelled clusters with N = 50 000 stars (where we take 0.1N as the core P2 stars), rh,0 = 0.5 pc, m3 = 1 M⊙, and one BH–BH binary of m1 = m2 = 30 M⊙. Then the ejection time scale, according to Eq. (10) and the typical value of α ≈ 10−1.8 (see Fig. 7), is τej ≈ 1.8 · 103 Myr. This is very long compared to τrh ≈ 16.5 Myr and τrc ≈ 8 Myr (see the 50k-30M* panels in Figs. 2 and 3). However, when we replace the stellar-mass BH with an IMBH of m1 = m2 = 300 M⊙, the encounter time scale from Eq. (1) becomes a thousand times smaller and τej ≈ 1.8 Myr. Comparing this with the 50k-10x30M* panels in Figs. 2 and 3, where the mass contained in the stellar mass BHs is equivalent to this IMBH–IMBH binary, we get τrc ≈ τej ≪ τrh. Therefore, such IMBHs could eject the P2 stars as effectively as ten stellar-mass BH binaries. Unlike stellar-mass BHs, the IMBHs could also rapidly merge due to gravitational wave emission and leave the GC core due to the recoil kick (Campanelli et al. 2007; González et al. 2007) that can have up to thousands of km/s (e.g. Mahapatra et al. 2024). The scattering of P2 stars would also make the IMBH binary more eccentric, helping with merging and further amplifying the kick (Sperhake et al. 2020). For context, the general relativistic merger time of an eccentric BH–BH binary is (see, e.g. Davies et al. 2011) τgr≈1.825⋅1012 Myr(abinau)4M⊙3m1m2(m1+m2)(1−e2)7/2.$\[\tau_{\mathrm{gr}} \approx 1.825 \cdot 10^{12} ~\mathrm{Myr}\left(\frac{a_{\mathrm{bin}}}{\mathrm{au}}\right)^4 \frac{M_{\odot}^3}{m_1 m_2\left(m_1+m_2\right)}\left(1-e^2\right)^{7 / 2}.\]$(17)

We also refer to Giersz et al. (2025), who argue that the change in the external tidal field and tidal stripping may be one of the main reasons for the P1/P2 radial inversion, albeit for all of their simulated clusters, such an inversion was only marginal and transient. Another mechanism that would lower the central potential well and push the central P2 population outward is gas expulsion (e.g. Decressin et al. 2010). However, gas loss would affect all stars in the GC. Thus, it is to be seen whether it could cause clear radial inversion of P1 and P2, since the scattering of P2 needs to be more localised in the centre. Lastly, we note that while it is conventionally assumed that P2 is born in the centres of GCs, there is still a possibility for another formation scenario of P2. All these additional aspects that could lead to faster mixing or the radial inversion of P1 and P2, which we were not able to explore in their entirety in our simple dynamical scenario, warrant future work.

thumbnail Fig. 10

Time evolution of the parameter ARlim+$\[A_{R\text {lim}}^{+}\]$ in the 10k-* clusters – the top row is for Rlim = 1 rh, the bottom one is for Rlim = 4 rh, see Eq. (16) for the definition. Only the models with the initial binary semi-major axis of 5 au (dashed lines) and 50 au (solid lines) are shown, averaged over all realisations of the models. The shaded area corresponds to their combined 1-sigma uncertainties. The A4+$\[A_{4}^{+}\]$ data points of NGC clusters are taken from Leitinger et al. (2025) for comparison.
thumbnail Fig. 11

Same as Fig. 10 but for the 50k-* models.

7 Conclusions

We studied the dynamics of stars inside star clusters to understand why certain observed, dynamically young GCs have their P1 more centrally concentrated than P2, some vice versa, and some GCs appear fully mixed. Using theoretical arguments and their testing in specifically designed numerical simulations of various clusters, we found that massive binary stars may be responsible for mixing of P1 and P2 early in the GC evolution. They heat the cluster core through binary–single star interactions and facilitate population redistribution on a shorter time scale than relaxation processes.

First, we note that while those high-mass binaries can push the centrally born P2 stars outwards, we do not observe a clear radial inversion of the P1 and P2 stars in any of our numerical setups. Consequently, we still cannot explain the nature of systems such as NGC 3201 and NGC 6101, where P1 is more centrally concentrated than P2. Provided, however, that only these two examples of inverted GCs are known to date and that they originated in a galaxy accreted by the Milky Way (Massari et al. 2019; Leitinger et al. 2025), this might suggest some specific evolution, such as tidal stripping (see, e.g. Giersz et al. 2025), which is currently beyond the scope of this paper.

More importantly, we conclude with our numerical analysis that high-mass primordial binaries can mix the GC members up to the half-mass radius and, in more compact systems, even to higher radii. While one massive binary is already able to affect the radial profile of the centrally-born P2 stars, full mixing of the populations is more effective if the models include several (e.g. ten or twenty) primordial binaries. It is also more rapid in the most compact models with the initial rh(0) = 0.5 pc rather than those with rh(0) = 1, 2 or 4 pc. Hence, the radial density profile of P1 and P2 is affected the most 1) in denser GCs, 2) by more massive binaries, and 3) by a higher binary count. The mechanism proposed in this paper can help us explain the observations of some dynamically young GCs, such as NGC 4590 or NGC 5904.

Acknowledgements

This research used computational resources from e-INFRA CZ (project ID:90254), supported by the Ministry of Education, Youth and Sports of the Czech Republic; and the computational cluster Virgo at the Astronomical Institute of the Czech Academy of Sciences. V.P. is funded by the European Union’s Horizon Europe and the Central Bohemian Region under the Marie Skłodowska-Curie Actions – COFUND, Grant agreement ID:101081195 (“MERIT”). V.P. also acknowledges support from the project RVO:67985815 at the Czech Academy of Sciences. E.I.L. acknowledges support from the ERC Consolidator Grant funding scheme (project ASTEROCHRONOMETRY, https://www.asterochronometry.eu, G.A. n. 772293). A.B. acknowledges support from the Australian Research Council (ARC) Centre of Excellence for Gravitational Wave Discovery (OzGrav), through project number CE230100016. M.H. acknowledges financial support from the Excellence Cluster ORIGINS which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC 2094 – 390783311. A.J.W. has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 101104656 and the Royal Society’s University Research Fellowship, reference URF/R1/241791. Last but not least, we thank the referee for valuable suggestions that improved this paper, in particular for pointing out the importance of the relative values of the core and half-mass relaxation time scales together with the time scale to eject stars from the core.

References

  1. Aarseth, S. J. 1972, Binary Evolution in Stellar Systems, ed. M. Lecar (Springer: The Netherlands), 88 [Google Scholar]
  2. Alessandrini, E., Lanzoni, B., Ferraro, F. R., Miocchi, P., & Vesperini, E. 2016, ApJ, 833, 252 [NASA ADS] [CrossRef] [Google Scholar]
  3. Banerjee, S. 2017, MNRAS, 467, 524 [NASA ADS] [Google Scholar]
  4. Bastian, N., & Lardo, C. 2018, ARA&A, 56, 83 [Google Scholar]
  5. Baumgardt, H., & Sollima, S. 2017, MNRAS, 472, 744 [Google Scholar]
  6. Bekki, K. 2011, MNRAS, 412, 2241 [NASA ADS] [CrossRef] [Google Scholar]
  7. Belczynski, K., Kalogera, V., Rasio, F. A., et al. 2008, ApJS, 174, 223 [Google Scholar]
  8. Belczynski, K., Bulik, T., Fryer, C. L., et al. 2010, ApJ, 714, 1217 [NASA ADS] [CrossRef] [Google Scholar]
  9. Bellini, A., Vesperini, E., Piotto, G., et al. 2015, ApJ, 810, L13 [NASA ADS] [CrossRef] [Google Scholar]
  10. Bianchini, P., Sills, A., & Miholics, M. 2017, MNRAS, 471, 1181 [Google Scholar]
  11. Binney, J., & Tremaine, S. 2008, Galactic Dynamics: 2nd edn. (Princeton: Princeton University Press) [Google Scholar]
  12. Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529 [Google Scholar]
  13. Cadelano, M., Dalessandro, E., & Vesperini, E. 2024, A&A, 685, A158 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Campanelli, M., Lousto, C. O., Zlochower, Y., & Merritt, D. 2007, Phys. Rev. Lett., 98, 231102 [Google Scholar]
  15. Carretta, E., Bragaglia, A., Gratton, R. G., et al. 2010, A&A, 520, A95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Cordoni, G., Milone, A. P., Marino, A. F., et al. 2020a, ApJ, 898, 147 [Google Scholar]
  17. Cordoni, G., Milone, A. P., Mastrobuono-Battisti, A., et al. 2020b, ApJ, 889, 18 [Google Scholar]
  18. Dalessandro, E., Cadelano, M., Vesperini, E., et al. 2019, ApJ, 884, L24 [Google Scholar]
  19. Dalessandro, E., Cadelano, M., Della Croce, A., et al. 2024, A&A, 691, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. Davies, M. B., Miller, M. C., & Bellovary, J. M. 2011, ApJ, 740, L42 [NASA ADS] [CrossRef] [Google Scholar]
  21. Decressin, T., Baumgardt, H., Charbonnel, C., & Kroupa, P. 2010, A&A, 516, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. D’Ercole, A., Vesperini, E., D’Antona, F., McMillan, S. L. W., & Recchi, S. 2008, MNRAS, 391, 825 [CrossRef] [Google Scholar]
  23. Faber, S. M., & Gallagher, J. S. 1979, ARA&A, 17, 135 [NASA ADS] [CrossRef] [Google Scholar]
  24. Fujii, M. S., & Portegies Zwart, S. 2014, MNRAS, 439, 1003 [NASA ADS] [CrossRef] [Google Scholar]
  25. Gieles, M., Charbonnel, C., Krause, M. G. H., et al. 2018, MNRAS, 478, 2461 [Google Scholar]
  26. Giersz, M., Askar, A., Hypki, A., et al. 2025, A&A, 698, L11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. González, J. A., Sperhake, U., Brügmann, B., Hannam, M., & Husa, S. 2007, Phys. Rev. Lett., 98, 091101 [NASA ADS] [CrossRef] [Google Scholar]
  28. Goodwin, S. P., & Kroupa, P. 2005, A&A, 439, 565 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Gratton, R. G., Carretta, E., & Bragaglia, A. 2012, A&A Rev., 20, 50 [NASA ADS] [CrossRef] [Google Scholar]
  30. Heggie, D. C. 1975, MNRAS, 173, 729 [NASA ADS] [CrossRef] [Google Scholar]
  31. Heggie, D. C., & Hut, P. 2003, The Gravitational Million-Body Problem: A Multidisciplinary Approach to Star Cluster Dynamics (Cambridge, UK: Cambridge University Press) [Google Scholar]
  32. Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543 [Google Scholar]
  33. Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897 [Google Scholar]
  34. Hut, P. 1983, ApJ, 272, L29 [NASA ADS] [CrossRef] [Google Scholar]
  35. Kroupa, P. 1995, MNRAS, 277, 1491 [Google Scholar]
  36. Kroupa, P. 2001, MNRAS, 322, 231 [NASA ADS] [CrossRef] [Google Scholar]
  37. Küpper, A. H. W., Maschberger, T., Kroupa, P., & Baumgardt, H. 2011, MNRAS, 417, 2300 [Google Scholar]
  38. Lacchin, E., Calura, F., Vesperini, E., & Mastrobuono-Battisti, A. 2022, MNRAS, 517, 1171 [NASA ADS] [CrossRef] [Google Scholar]
  39. Leitinger, E., Baumgardt, H., Cabrera-Ziri, I., Hilker, M., & Pancino, E. 2023, MNRAS, 520, 1456 [NASA ADS] [CrossRef] [Google Scholar]
  40. Leitinger, E. I., Baumgardt, H., Cabrera-Ziri, I., et al. 2025, A&A, 694, A184 [Google Scholar]
  41. Libralato, M., Bellini, A., Piotto, G., et al. 2019, ApJ, 873, 109 [NASA ADS] [CrossRef] [Google Scholar]
  42. Livernois, A. R., Aros, F. I., Vesperini, E., et al. 2024, MNRAS, 534, 2397 [Google Scholar]
  43. Lyne, A. G., & Lorimer, D. R. 1994, Nature, 369, 127 [NASA ADS] [CrossRef] [Google Scholar]
  44. Mahapatra, P., Favata, M., & Arun, K. G. 2024, Phys. Rev. D, 110, 084041 [Google Scholar]
  45. Martocchia, S. 2020, PhD thesis, Liverpool John Moores University, UK [Google Scholar]
  46. Martocchia, S., Cabrera-Ziri, I., Lardo, C., et al. 2018, MNRAS, 473, 2688 [NASA ADS] [CrossRef] [Google Scholar]
  47. Massari, D., Koppelman, H. H., & Helmi, A. 2019, A&A, 630, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  48. Merritt, D. 2013, Class. Quantum Grav., 30, 244005 [Google Scholar]
  49. O’Leary, R. M., Stahler, S. W., & Ma, C.-P. 2014, MNRAS, 444, 80 [CrossRef] [Google Scholar]
  50. Pavlík, V. 2020, A&A, 638, A155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Pavlík, V., & Subr, L. 2018, A&A, 620, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Pavlík, V., & Vesperini, E. 2021, MNRAS, 504, L12 [CrossRef] [Google Scholar]
  53. Pavlík, V., & Vesperini, E. 2022a, MNRAS, 509, 3815 [Google Scholar]
  54. Pavlík, V., & Vesperini, E. 2022b, MNRAS, 515, 1830 [Google Scholar]
  55. Pavlík, V., Jeřábková, T., Kroupa, P., & Baumgardt, H. 2018, A&A, 617, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Pavlík, V., Kroupa, P., & Šubr, L. 2019, A&A, 626, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Pavlík, V., Heggie, D. C., Varri, A. L., & Vesperini, E. 2024, A&A, 689, A313 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  58. Pavlík, V., Davies, M. B., Leitinger, E. I., et al. 2025, in Proceedings of MODEST-25 / IAU Symposium 398 [arXiv:2508.19747] [Google Scholar]
  59. Peuten, M., Zocchi, A., Gieles, M., Gualandris, A., & Hénault-Brunet, V. 2016, MNRAS, 462, 2333 [CrossRef] [Google Scholar]
  60. Plummer, H. C. 1911, MNRAS, 71, 460 [Google Scholar]
  61. Plunkett, A. L., Fernández-López, M., Arce, H. G., et al. 2018, A&A, 615, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  62. Podsiadlowski, P., Langer, N., Poelarends, A. J. T., et al. 2004, ApJ, 612, 1044 [NASA ADS] [CrossRef] [Google Scholar]
  63. Sollima, A. 2021, MNRAS, 502, 1974 [NASA ADS] [CrossRef] [Google Scholar]
  64. Sperhake, U., Rosca-Mead, R., Gerosa, D., & Berti, E. 2020, Phys. Rev. D, 101, 024044 [Google Scholar]
  65. Spitzer, Jr., L. 1987, Dynamical Evolution of Globular Clusters (Princeton, USA: Princeton University Press) [Google Scholar]
  66. Spitzer, Jr., L., & Hart, M. H. 1971, ApJ, 164, 399 [NASA ADS] [CrossRef] [Google Scholar]
  67. Tanikawa, A., Hut, P., & Makino, J. 2012, New A, 17, 272 [NASA ADS] [CrossRef] [Google Scholar]
  68. Tanikawa, A., Heggie, D. C., Hut, P., & Makino, J. 2013, Astron. Comput., 3, 35 [Google Scholar]
  69. Tiongco, M. A., Vesperini, E., & Varri, A. L. 2016, MNRAS, 455, 3693 [Google Scholar]
  70. Tiongco, M. A., Vesperini, E., & Varri, A. L. 2019, MNRAS, 487, 5535 [NASA ADS] [CrossRef] [Google Scholar]
  71. Vesperini, E., McMillan, S. L. W., D’Antona, F., & D’Ercole, A. 2013, MNRAS, 429, 1913 [Google Scholar]
  72. Vesperini, E., Hong, J., Giersz, M., & Hypki, A. 2021, MNRAS, 502, 4290 [NASA ADS] [CrossRef] [Google Scholar]
  73. Wang, L., Iwasawa, M., Nitadori, K., & Makino, J. 2020a, MNRAS, 497, 536 [Google Scholar]
  74. Wang, L., Kroupa, P., Takahashi, K., & Jerabkova, T. 2020b, MNRAS, 491, 440 [Google Scholar]
  75. Yaghoobi, A., Calura, F., Rosdahl, J., & Haghi, H. 2022a, MNRAS, 510, 4330 [Google Scholar]
  76. Yaghoobi, A., Rosdahl, J., Calura, F., Khalaj, P., & Haghi, H. 2022b, MNRAS, 517, 4175 [Google Scholar]
  77. Ziosi, B. M., Mapelli, M., Branchesi, M., & Tormen, G. 2014, MNRAS, 441, 3703 [Google Scholar]
  78. Zonoozi, A. H., Rabiee, M., Haghi, H., & Kroupa, P. 2024, ApJ, 975, 266 [Google Scholar]
1

We note that the core relaxation time may be slightly different depending on the stellar density in the core and the velocity distribution (see, e.g. Binney & Tremaine 2008; Merritt 2013; Pavlík et al. 2024).

2

The half-light and half-mass radii coincide in an equal-mass cluster.

3

Ten realisations of each smaller model and five of the larger models.

4

We note that the primordial binaries were not counted towards either population, but they did contribute to the potential term in Eq. (13).

5

We note that these are typical BH masses and that some detected gravitational-wave sources even exceed 30 M⊙. The upper limit of all massive objects per cluster also did not exceed the constraints coming from more realistic systems. Specifically, applying stellar evolution (from Hurley et al. 2000, 2002) on the Kroupa (2001) initial mass function (from 0.08 to 150 M⊙) results in ≈6% of the initial total GC mass being contained in BH progenitors.

6

See Eq. (3) where the theoretical hard/soft limit is a ≈ 200 au for a (5+5) M⊙ binary and ≈8000 au for a (30+30) M⊙ binary.

7

We note that while in the work of Pavlík et al. (2025), we showed almost full radial mixing of the clusters with 104 stars and ten massive primordial binaries, those initial conditions corresponded to a very idealised toy model. In the models presented here, which have more realistic velocity distributions, the mixing caused by the binary is slightly suppressed by the random radial motions of the single stars.

All Figures

thumbnail Fig. 1

Estimated semi-major axis (top) and the encounter time scale (bottom) of a binary star that can eject single stars from the core in a cluster of a given mass and half-mass radius. In this plot, we used Eqs. (8) and (9), for clusters containing stars of equal masses (m3 = 1 M⊙), one binary star with m1 = m2 = 30 M⊙ at its centre, and v∞ = 10 km s−1. Time in the bottom panel is normalised by the half-mass relaxation time, Eq. (12).
In the text
thumbnail Fig. 2

Evolution of the half-mass relaxation time, τrh(t), calculated from Eq. (12) in terms of the initial half-mass relaxation time, τrh(0). The left-hand set of panels is for the 10k-* clusters, the right-hand set is for 50k-*. The models are separated by the primordial binary star numbers and masses (columns), the binary star initial semi-major axes (colours), and the initial cluster half-mass radius (line styles). Each line is averaged over all realisations of the corresponding model. (We note that the models 10k-5M* and 50k-30M* do not show any differences from the models without binaries, which we also integrated for comparison but do not display them here.
In the text
thumbnail Fig. 3

Evolution of the core relaxation time, τrc(t) in the units of τrh(t), calculated from Eq. (11). The horizontal axis shows time in terms of the initial half-mass relaxation time, τrh(0). The left-hand set of panels is for the 10k-* clusters, the right-hand set is for 50k-*. The models are separated by the primordial binary star numbers and masses (columns), the binary star initial semi-major axes (colours), and the initial cluster half-mass radius (line styles). Each line is averaged over all realisations of the corresponding model.
In the text
thumbnail Fig. 4

Number of massive binaries in each 10k-* model in time (averaged over each model’s realisations). The models are separated by their initial rh (columns), primordial binary systems (rows), and their semi-major axes (colours). In each plot, the top horizontal axis shows the evolutionary time in the initial half-mass relaxation times, and the lower horizontal axis is in multiples of τrh at the end of the simulation (to mimic the ‘present-day’ value). While we permit the dynamical formation of the low-mass (m3+m3) or mixed (m1,2+m3) binaries in the simulations, the numbers here only show the massive (m1+m2) binaries, including component switching. However, if the primordial binary goes through a state of an unstable triple, it is not counted towards Nbin (which is why we sometimes see the number of binaries decreasing and then increasing again).
In the text
thumbnail Fig. 5

Same as Fig. 4 but for the 50k-* models.
In the text
thumbnail Fig. 6

Probability density function (PDF) of the binding energy change, α, per interaction, in the 10k-* models. Binary–single and binary–binary interactions are counted together. For reference, we also plot with dashed vertical lines the theoretical values of α, for which a single m3-star should theoretically escape, as derived using Eq. (8). The figure is separated by each model’s initial half-mass radius (columns), number and masses of the primordial binaries (rows), and the binary semi-major axes (colours).
In the text
thumbnail Fig. 7

Same as Fig. 6 but for the 50k-* models.
In the text
thumbnail Fig. 8

Selected 10k-* models with the highest A+ at 4 τrh for each primordial binary set-up (as labelled in the top left corners). Top panels: normalised projected cumulative radial distributions of stars in two populations (P1 and P2 – distinguished by colours) in the clusters, calculated from Eq. (16). All realisations of the corresponding model are plotted. Columns correspond to different values of Rlim. Bottom panels: the ratio of P2 stars in several radial bins (similarly to Leitinger et al. 2023). Each model realisation is plotted with a black line. The parameter ARlim+$\[A_{R\text {lim}}^{+}\]$ (its maximum range and median with standard deviation) is also displayed for each region.
In the text
thumbnail Fig. 9

Same as Fig. 8 but for the 50k-* models.
In the text
thumbnail Fig. 10

Time evolution of the parameter ARlim+$\[A_{R\text {lim}}^{+}\]$ in the 10k-* clusters – the top row is for Rlim = 1 rh, the bottom one is for Rlim = 4 rh, see Eq. (16) for the definition. Only the models with the initial binary semi-major axis of 5 au (dashed lines) and 50 au (solid lines) are shown, averaged over all realisations of the models. The shaded area corresponds to their combined 1-sigma uncertainties. The A4+$\[A_{4}^{+}\]$ data points of NGC clusters are taken from Leitinger et al. (2025) for comparison.
In the text
thumbnail Fig. 11

Same as Fig. 10 but for the 50k-* models.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.