Black holes and core expansion in massive star clusters

1 INTRODUCTION

• Young clusters typically have core radii of 1–2 pc, whereas the spread in core radius increases with age, reaching approximately 0–8 pc among the oldest clusters in earlier HST measurements.The same radius–age trend was found in both the LMC and SMC.Source: §1 · p. 2
• The oldest Magellanic Cloud clusters span core radii from approximately 0 to 17 pc, while the upper envelope is well defined only up to a few Gyr.Several old clusters are represented by upper limits because of severe crowding; all affected clusters have measured rc < 1 pc.Source: §1 · p. 2
• Old Magellanic Cloud clusters have core-radius distributions similar in range and overall shape to Galactic globular clusters, especially to the distant young-halo subsystem.The full Galactic sample has a sharper peak at small core radii, consistent with inner-Galaxy tidal forces destroying loosely bound clusters more rapidly.Source: §1 · pp. 2–3, Figure 2 · p. 3
• Standard models do not predict the observed order-of-magnitude core-radius expansion over a Hubble time, motivating a search for its origin.The authors note that changes in cluster-formation conditions over roughly 10 Gyr are another possibility, but available constraints are limited and they do not investigate it further.Source: §1 · p. 3
• Proposed explanations include early stellar-evolution mass loss from flat-IMF clusters, but the required mass loss would disrupt clusters within tens of Myr and cannot explain expansion across the observed age range.The mass-loss phase lasts only about the first 100 Myr, and observational evidence indicates that young-cluster IMFs are more or less invariant.Source: §1 · p. 3
• Time-varying tides and hard primordial binaries produced too little core expansion in earlier small-cluster simulations to account for the observed trend.The models also found similar core-radius evolution on circular and elliptical orbits.Source: §1 · p. 3
• The authors conduct large-scale, realistic N-body simulations to test whether black-hole heating and primordial mass segregation can provide an internal dynamical origin for the trend.The motivation is strengthened by the similar trend across the LMC and SMC and the consistency of its endpoints with globular clusters in several galaxies with different tidal fields.Source: §1 · p. 4

2.1 N-body code and initial conditions

• Initial clusters use EFF profiles with γ = 3, a Kroupa (2001) IMF, and a 0.1–100 M⊙ stellar mass range, yielding initial masses near log Mtot ∼ 4.75.The chosen profile and initial conditions were designed to match young Magellanic Cloud clusters.Source: §2.1 · pp. 4–5
• With N ≈ 10^5, direct N-body models follow stellar and binary evolution alongside cluster dynamics, enabling simulated observations without scaling results from smaller-N systems.The models use nbody4; large particle counts also reduce fluctuations in global evolution.Source: §2.1 · pp. 4–6
• Models orbit at 6 kpc in a point-mass LMC tidal field, with stars removed at 2rt; the authors caution that escape prescriptions can change evaporation rates and cluster lifetimes.The field is imposed to simulate gradual evaporation and compare models under a common external potential and escape criterion.Source: §2.1 · p. 6
• An initial length scale of 8.26 pc gives unsegregated models log ρ0 = 2.31 and rc = 1.90 pc, matching values observed for many young Magellanic Cloud clusters.The models assume initial tidal-radius filling; EFF profiles are truncated at the estimated tidal radius when generating stellar positions.Source: §2.1 · p. 6
• Each model forms 198 black holes with masses uniformly distributed from 8 to 12 M⊙, while natal kicks control the retained fraction.All stars initially above 20 M⊙ form black holes in these simulations; the code allows varying kick prescriptions and retention fractions.Source: §2.1 · pp. 6–7

2.2 Primordial mass segregation

• Primordial mass segregation is set by pre-evolving clusters with stellar masses truncated above 8 M⊙, then restoring original masses before full runs; longer pre-evolution produces smaller cores and higher central densities.The 450 Myr models have rc = 0.25 pc and log ρ0 = 4.58, similar to R136; 115 and 225 Myr models have intermediate structural properties. The models omit primordial binaries and assume solar metallicity.Source: §2.2 · pp. 7–8, §2.1 · p. 7

3 “OBSERVING” THE SIMULATIONS

• Simulated observations derive observational core radii at each output time, making model measurements directly comparable to those for real Magellanic Cloud clusters.The procedure incorporates subtleties of the HST measurements used to define the observed radius-age trend.Source: §3 · p. 8
• The sampled stellar-mass range changes with cluster age, so the simulations impose the bright and faint detection limits appropriate to each output time.The limits reflect saturation and approximate 50 per cent completeness; Figure 4 separates clusters into four age bins and distinguishes WFPC2 and ACS imaging.Source: §3 · p. 8, Figure 4 · p. 8
• The pipeline converts stellar luminosities and temperatures to standard magnitudes, applies detection limits, builds projected surface-brightness profiles, and fits EFF models within the camera field of view.Measurements average results from three orthogonal projections; the profiles are truncated to mimic the limited WFPC2 and ACS fields.Source: §3 · pp. 8–9
• The observational core radius rc is the projected half-central-brightness radius, whereas the N-body density radius rd has no general relationship to rc and may evolve differently.EFF fitting remains appropriate because the profiles are truncated to match the observations and simulated escapers can populate radii beyond rt.Source: §3 · p. 9

4 SIMULATIONS AND RESULTS

• The modeled clusters resemble young Magellanic Cloud clusters: unsegregated models have rc ∼1.9 pc and log ρ0 ∼2.3, while strongly segregated models have rc ∼0.25 pc and log ρ0 ∼4.6.The strongly segregated models also have γ ∼2.3 and resemble R136; the unsegregated models resemble several clusters around 20 Myr old. The authors conclude the simulations directly model massive Magellanic Cloud cluster evolution.Source: §4 · p. 9

4.1 Runs 1 and 2: No mass segregation

• Run 1 loses about 25% of its initial mass by 100 Myr without a corresponding core-radius response, then contracts through relaxation to core collapse at 17.4 Gyr.At greatest collapse, its central density reaches log ρ0 ≈4.5; the early mass loss is distributed across the cluster, while the core radius later shrinks as luminous stars segregate inward.Source: §4.1 · p. 10
• At ≈650 Myr, Run 2’s core radius begins expanding and continues through its 10.67 Gyr endpoint, while black-hole-free Run 1 contracts toward core collapse.The runs differ only in black-hole natal kicks and retention; Run 2 retains all 198 black holes. Its radius-age evolution runs parallel to, but begins too late to trace exactly, the observed upper envelope.Source: §4.1 · pp. 11–13, Figure 5 · p. 11
• Black holes sink to the centre, form binaries, and scatter or eject other black holes, heating the stellar core and suppressing further stellar mass segregation.The first black-hole binary forms at about 510 Myr, and several are present by about 650 Myr, when structural evolution diverges from the control; close interactions between stars and black-hole binaries contribute negligible heating.Source: §4.1 · pp. 12–15
• By the end of Run 2, rc/rh,l reaches ≈0.8, matching several extended Magellanic Cloud clusters; stellar-mass black holes can therefore produce large ratios also associated with proposed central IMBHs.Run 1 stays near rc/rh,l ≈0.45 for much of its evolution before the ratio falls during core collapse.Source: §4.1 · p. 14, Figure 9 · p. 14
• In Run 2, declining black-hole subsystem density slows encounters, ejections, and stellar-core heating, prolonging the subsystem’s survival and producing roughly logarithmic core-radius growth.At 10.67 Gyr, 65 single and 2 binary black holes remained in the cluster.Source: §4.1 · p. 15
• Run 2 ejected 15 black-hole binaries, typically at 4–6 AU before 2.5 Gyr and 10–30 AU after 5 Gyr, as cluster expansion lowered its escape velocity.Their mean black-hole mass was 10.98 M⊙; one triple black-hole system was also ejected at about 4100 Myr.Source: §4.1 · p. 15, Figure 10 · p. 15
• None of Run 2’s ejected black-hole binaries would merge within 12 Gyr, although modestly tighter or more eccentric orbits could meet that timescale.The most tightly bound had a_b = 4.2 AU and e = 0.839; the most eccentric had a_b = 8.2 AU and e = 0.972.Source: §4.1 · p. 16
• Run 2 showed no significant stellar density or velocity cusps at 5 or 10 Gyr; continual passage of scattered and ejected black holes may disrupt cusp formation.Only a very marginal density excess above the constant-density-core models was possible, with low significance.Source: §4.1 · pp. 16–17, Figure 11 · p. 16
• Run 2’s stellar velocity-dispersion profiles imply unseen central mass: fitting them required central densities about 3.4 and 2.1 times those inferred from density profiles at 5 and 10 Gyr.Projected 5-Gyr data also favored the higher-density model, and samples of upper-main-sequence and red-giant stars favored the larger mass-to-light ratio, though measurement is challenging.Source: §4.1 · pp. 17–18, Figure 12 · p. 17

4.2 Runs 3 and 4: Strong mass segregation

• In strongly segregated Run 3, central density fell rapidly from an R136-comparable initial value to match other young Magellanic Cloud clusters by about 10 Myr, indicating early core expansion.Run 1’s lower central density changed little during early stellar evolution, and Runs 1 and 3 together spanned observed densities of the youngest clusters.Source: §4.2 · p. 18, Figure 13 · p. 18
• Run 3’s simulated radial mass- and luminosity-function measurements generally matched observations of R136, NGC 1805, NGC 1818, and NGC 330, supporting the constructed primordial segregation.The model evolved from an R136-like profile to resemble NGC 1805 and NGC 1818 at 15 Myr and NGC 330 at 30 Myr through internal dynamics and rapid stellar-evolution mass loss.Source: §4.2 · p. 19, Figure 14 · p. 19
• Retained black holes restart Run 4’s core expansion around 750 Myr, reaching rc ≈11 pc by 10 Gyr; models exceeding ∼6 pc after a Hubble time require prolonged black-hole-driven expansion.Early mass loss alone acts for too short a time to produce that long-term expansion.Source: Figure 15 · p. 20, §4.2 · pp. 21–22
• Centrally concentrated early stellar-evolution mass loss heats mass-segregated cores and drives dramatic expansion; Run 3 expands more than Run 4 because its 198 black holes escape and it loses more mass.The supernovae, and thus the mass loss, are more centrally concentrated in Run 3 than in Run 1.Source: §4.2 · p. 20, Figure 16 · p. 21
• Early expansion pushes Run 4’s centrally formed black holes outward before they sink inward; black-hole binary formation therefore begins at similar times in segregated and non-segregated models.In Run 4, the first black-hole binary forms at ∼570 Myr.Source: §4.2 · pp. 21–22

4.3 Runs 4a and 4b: Variable mass segregation

• Runs 4a and 4b have two core-expansion phases, with greater early expansion at higher primordial mass segregation; their tracks lie between Runs 2 and 4 and later run approximately parallel.The first phase follows rapid stellar evolution and the second begins around 600–800 Myr with retained black-hole activity.Source: §4.3 · pp. 23–24, Figure 21 · p. 23
• First black-hole binaries formed at 510, 570, 540, and 460 Myr in Runs 2, 4, 4a, and 4b, respectively, despite their different early structural evolution.The models’ early expansion delays central black-hole accumulation, so primordial mass segregation does not make the black-hole subsystem dynamically active significantly earlier.Source: §4.3 · p. 24

4.4 Run 5: Intermediate BH retention

• Run 5 was evolved to 10.06 Gyr with strong primordial mass segregation and approximately 0.5 black-hole retention (96 black holes after 100 Myr), testing whether expansion requires near-total retention.Its initial conditions otherwise matched Runs 3 and 4, which had zero and full black-hole retention, respectively.Source: §4.4 · p. 24, Figure 24 · p. 25
• Run 5 retains an estimated fBH = 0.485, or as little as 0.440 under an alternative escape-counting definition; natal-kick ejections make the fraction stochastic.All 198 black holes form by τ = 10 Myr; 102 have escaped by τ = 100 Myr. Natal kicks also give retained black holes extra kinetic energy, delaying their inward sinking.Source: §4.4 · pp. 25–26
• Run 5’s black-hole heating begins late and proceeds more slowly: its first black-hole binary forms at τ = 1350 Myr, and subsequent hardening, scattering, and ejection are reduced.The first post-binary single-black-hole ejection occurs at τ = 2050 Myr, and the first binary ejection at τ = 4000 Myr.Source: §4.4 · p. 26
• Despite retaining only about half as many black holes as Run 4, Run 5 reaches rc/rh,l ∼0.8 after ∼10 Gyr, demonstrating that complete black-hole retention is unnecessary for prolonged core expansion.The late expansion begins at τ ≈1400 Myr; after it begins, Run 5’s rc/rh,l locus lies on top of Run 4’s, although its core-radius expansion is slower.Source: §4.4 · pp. 25–26

4.5 Run 6: Can neutron stars replace BHs?

• Retaining 425 neutron stars in Run 6 produces no core expansion: without runaway subsystem collapse or central neutron-star binaries, the cluster contracts toward core collapse.Neutron stars sink toward the centre on a ∼500 Myr timescale, but are not sufficiently more massive than other cluster members to make their subsystem unstable.Source: §4.5 · p. 27, Figure 27 · p. 27
• Run 6 reaches deepest collapse at τ ≈12.8 Gyr, earlier than Run 1 at τ ≈17.4 Gyr; the respective points occur after 4.40 and 8.37 median relaxation times.The difference is linked to shorter central relaxation times in Run 6, where the centralmost objects are predominantly neutron stars. Binary formation, including NS-NS binaries, accompanies deepest collapse.Source: §4.5 · p. 27
• During collapse, Run 6’s measured core radius is larger than Run 1’s because its central neutron stars contribute no light; rc/rh,l oscillates around ∼0.2 rather than ≲0.05.The luminosity-based core-radius calculation therefore traces a more widely distributed luminous core in Run 6 than in Run 1.Source: §4.5 · p. 27

• Early stellar-evolution mass loss and later heating by retained stellar-mass black holes together produce varied evolutionary paths spanning the observed Magellanic Cloud radius-age distribution.The simulations identify these as physically motivated drivers of the observed trend.Source: §5 · p. 28
• Early expansion depends on centrally concentrated mass loss: greater primordial mass segregation or a larger mass loss produces faster expansion, while unsegregated models show essentially none.The rapid phase is mostly complete by 100–200 Myr; its effects can persist when relaxation times become long.Source: §5 · pp. 28–29
• The models omit residual-gas expulsion, which may affect clusters younger than about 50 Myr; jointly modelling it with early stellar mass loss requires more sophisticated codes.The authors also suggest that expansion near 100 Myr should be considered when interpreting young-cluster dynamical masses, although its integrated-spectra signal is unclear.Source: §5 · pp. 29, 31
• Black-hole-driven expansion begins only after a dense central subsystem forms binaries, whose interactions scatter and eject black holes and heat the cluster core.In otherwise similar models, the first black-hole binaries form at roughly 500–600 Myr; natal kicks can delay their formation and reduce the retained population.Source: §5 · p. 29
• Retained black holes can sustain late expansion: long-duration models retain sizeable populations after a Hubble time, and expansion continues even with roughly 50% retention.Expansion proceeds more slowly with fewer black holes or a lower-density stellar core; by about 10 Gyr, the core-to-half-light-radius ratio approaches 0.8 in expanding models.Source: §5 · pp. 29–30
• Early stellar-evolution mass loss alone cannot explain the oldest, most extended clusters; the models require retained black holes for that expansion, which neutron stars do not reproduce.Without sufficient black-hole retention, a cluster can contract again as two-body relaxation takes over.Source: §5 · p. 30
• The proposed radius-age explanation requires differences in black-hole retention among similar clusters, potentially arising from sampling, progenitor evolution, cluster mass, or natal kicks.The authors associate very large old cores with significant retained populations and late core collapse with few retained black holes.Source: §5 · p. 30
• Core expansion may alter cluster-survival and formation studies because prolonged expansion can increase the number of diffuse, tidally vulnerable clusters and may help form unusually extended systems.The authors also caution that prolonged expansion makes the half-mass radius unstable, complicating its use to infer cluster formation.Source: §5 · p. 31

• EFF models may provide more appropriate initial conditions than Plummer spheres or King models for some scenarios, and the appendix supplies formulae to generate them.The authors note that such models become more feasible as computing power and software sophistication increase.Source: Appendix A · p. 35

A1 General properties

• EFF models were empirically fitted to young massive LMC clusters, whose outer regions lack tidal cut-offs.They are represented by projected three-parameter profiles.Source: A1 General properties · p. 35
• Deprojection preserves the density profile’s functional form but changes its index from γ to γ + 1.ρ0 denotes central luminosity density; multiply by the global mass-to-light ratio Υ to obtain central mass density.Source: A1 General properties · p. 35
• Integrating the density profile gives expressions for the gravitational potential and enclosed mass, with the potential involving a Gauss hypergeometric function.The enclosed quantity is mass when ρ0 is a mass density, or luminosity when it is a luminosity density.Source: A1 General properties · p. 35
• For γ > 2, the asymptotic mass converges to M∞ = 2πµ0Υa2/(γ −2).For γ ≤ 2, the expression does not converge as r approaches infinity.Source: A1 General properties · p. 35
• The Jeans equations yield radial isotropic velocity dispersion, and integer γ values give elementary analytic expressions for M(r) and σ2(r).These expressions can be used in code to generate initial conditions; γ = 4 is the Plummer sphere, while γ = 3, 5, and 6 are also considered.Source: A1 General properties · p. 35

A2 The γ = 3 case

• The γ = 3 EFF case has analytic M(r) and σ2(r) expressions and, in projection, most closely matches observed young massive star cluster profiles in the LMC and SMC.The 18 young LMC and SMC clusters span γ = 2.05–3.79, with a median γ = 2.67.Source: A2 The γ = 3 case · pp. 35–36
• For γ = 3, the hypergeometric function reduces to an explicit form, yielding analytic enclosed-mass and isotropic-velocity-dispersion expressions used for the initial conditions.The section presents the substitutions and resulting expressions.Source: A2 The γ = 3 case · p. 36

A3 Steeper cases: γ = 5 and γ = 6

• The γ = 5 and γ = 6 models are motivated by observed globular clusters and a Fornax diffuse cluster with steep brightness-profile fall-offs.Reported γ values include 4.55 ± 0.61 for NGC 1841, 5.21 ± 0.99 for NGC 339, and 7.52 ± 0.64 for Fornax diffuse cluster 1.Source: A3 Steeper cases: γ = 5 and γ = 6 · p. 36

γ = 5

• For γ = 5, the hypergeometric function is reduced and expressions are derived for enclosed mass and isotropic velocity dispersion.The text identifies the mass expression after substitution into Eq. A4 and the velocity-dispersion expression after evaluating Eq. A5.Source: A3 Steeper cases: γ = 5 and γ = 6 · p. 36

γ = 6

• The γ = 6 case is highlighted because its properties are comparable in simplicity to those of the Plummer sphere (γ = 4).The hypergeometric function in Eqs. A3 and A4 can be written for γ = 6.Source: A3 Steeper cases: γ = 5 and γ = 6 · p. 36

A3 Steeper cases: γ = 5 and γ = 6

• Although similar analytic expressions can be derived for increasing integer γ, they become correspondingly more complicated.Source: A3 Steeper cases: γ = 5 and γ = 6 · p. 36