Our current understanding of the stellar initial mass function and massive star evolution suggests that young globular clusters may have formed hundreds to thousands of stellar-mass black holes, the remnants of stars with initial masses from ∼ 20 − 100 M . Birth kicks from supernova explosions may eject some black holes from their birth clusters, but most should be retained. Using a Monte Carlo method we investigate the long-term dynamical evolution of globular clusters containing large numbers of stellar black holes. We describe numerical results for 42 models, covering a range of realistic initial conditions, including up to 1.6 × 10 6 stars. In almost all models we find that significant numbers of black holes (up to ∼ 10 3 ) are retained all the way to the present. This is in contrast to previous theoretical expectations that most black holes should be ejected dynamically within a few Gyr. The main reason for this difference is that core collapse driven by black holes (through the Spitzer "mass segregation instability") is easily reverted through three-body processes, and involves only a small number of the most massive black holes, while lower-mass black holes remain well-mixed with ordinary stars far from the central cusp. Thus the rapid segregation of stellar black holes does not lead to a long-term physical separation of most black holes into a dynamically decoupled inner core, as often assumed previously. Combined with the recent detections of several black hole X-ray binary candidates in Galactic globular clusters, our results suggest that stellar black holes could still be present in large numbers in many globular clusters today, and that they may play a significant role in shaping the long-term dynamical evolution and the present-day dynamical structure of many clusters.

Globular clusters should be born with significant numbers of stellar-mass black holes (BHs). It has been thought for two decades that very few of these BHs could be retained through the cluster lifetime. With masses ∼ 10 M ⊙ , BHs are ∼ 20 times more massive than an average cluster star. They segregate into the cluster core, where they may eventually decouple from the remainder of the cluster. The small-N core then evaporates on a short timescale. This is the so-called Spitzer instability. Here we present the results of a full dynamical simulation of a globular cluster containing many stellar-mass BHs with a realistic mass spectrum. Our Monte Carlo simulation code includes detailed treatments of all relevant stellar evolution and dynamical processes. Our main finding is that old globular clusters could still contain many BHs at present. In our simulation, we find no evidence for the Spitzer instability. Instead, most of the BHs remain well-mixed with the rest of the cluster, with only the innermost few tens of BHs segregating significantly. Over the 12 Gyr evolution, fewer than half of the BHs are dynamically ejected through strong binary interactions in the cluster core. The presence of BHs leads to long-term heating of the cluster, ultimately producing a core radius on the high end of the distribution for Milky Way globular clusters (and those of other galaxies). A crude extrapolation from our model suggests that the BH-BH merger rate from globular clusters could be comparable to the rate in the field.

We study the dynamical evolution of globular clusters containing primordial binaries, including full single and binary stellar evolution using our Monte Carlo cluster evolution code updated with an adaptation of the single and binary stellar evolution codes SSE/BSE from Hurley et al. (2000Hurley et al. ( , 2002. We describe the modifications we have made to the code. We present several test calculations and comparisons with existing studies to illustrate the validity of the code. We show that our code finds very good agreement with direct N-body simulations including primordial binaries and stellar evolution. We find significant differences in the evolution of the global properties of the simulated clusters using stellar evolution compared to simulations without any stellar evolution. In particular, we find that the mass loss from stellar evolution acts as a significant energy production channel simply by reducing the total gravitational binding energy and can significantly prolong the initial core contraction phase before reaching the binaryburning quasi steady state of the cluster evolution as noticed in Paper IV. We simulate a large grid of clusters varying the initial cluster mass, binary fraction, and concentration and compare properties of the simulated clusters with those of the observed Galactic globular clusters (GGCs). We find that our simulated cluster properties agree well with the observed GGC properties. We explore in some detail qualitatively different clusters in different phases of their evolution, and construct synthetic Hertzprung-Russell diagrams for these clusters.

We present a new parallel code for computing the dynamical evolution of collisional N-body systems with up to N ∼ 10 7 particles. Our code is based on the Hénon Monte Carlo method for solving the Fokker-Planck equation, and makes assumptions of spherical symmetry and dynamical equilibrium. The principal algorithmic developments involve optimizing data structures, and the introduction of a parallel random number generation scheme, as well as a parallel sorting algorithm, required to find nearest neighbors for interactions and to compute the gravitational potential. The new algorithms we introduce along with our choice of decomposition scheme minimize communication costs and ensure optimal distribution of data and workload among the processing units. Our implementation uses the Message Passing Interface (MPI) library for communication, which makes it portable to many different supercomputing architectures. We validate the code by calculating the evolution of clusters with initial Plummer distribution functions up to core collapse with the number of stars, N, spanning three orders of magnitude, from 10 5 to 10 7 . We find that our results are in good agreement with self-similar core-collapse solutions, and the core collapse times generally agree with expectations from the literature. Also, we observe good total energy conservation, within 0.04% throughout all simulations. We analyze the performance of the code, and demonstrate near-linear scaling of the runtime with the number of processors up to 64 processors for N = 10 5 , 128 for N = 10 6 and 256 for N = 10 7 . The runtime reaches saturation with the addition of processors beyond these limits, which is a characteristic of the parallel sorting algorithm.The resulting maximum speedups we achieve are approximately 60×, 100×, and 220×, respectively.

We investigate the dynamical decay of non-hierarchical accreting triple systems and its implications on the ejection model as Brown Dwarf formation scenario. A modified chain-regularization scheme is used to integrate the equations of motion, that also allows for mass changes over time as well as for momentum transfer from the accreted gas mass onto the bodies. We integrate an ensemble of triple systems within a certain volume with different accretion rates, assuming several prescriptions of how momentum is transferred onto the bodies. We follow their evolution until the systems have decayed. We analyze the end states and decay times of these systems and determine the fraction of Brown Dwarfs formed, their escape speeds as well as the semi-major axis distribution of the formed Brown Dwarf binaries. We find that the formation probability of Brown Dwarfs depends strongly on the assumed momentum transfer which is related to the motion of the gas. Due to ongoing accretion and consequent shrinkage of the systems, the median escape velocity is increased by a factor of 2 and the binary separations are decreased by a factor of 5 compared with non-accreting systems. Furthermore, the obtained semi-major axis distribution drops off sharply to either side of the median, which is also supported by observations. We conclude that accretion and momentum transfer of accreted gas during the dynamical decay of triple systems is able to produce the observed distribution of close binary Brown Dwarfs, making the ejection model a viable option as Brown Dwarf formation scenario.Comment: 31 pages, 8 figures, accepted for publication in Ap