The "gravitational million-body problem," to model the dynamical evolution of a self-gravitating, collisional N-body system with ~10^6 particles over many relaxation times, remains a major challenge in computational astrophysics. Unfortunately, current techniques to model such systems suffer from severe limitations. A direct N-body simulation with more than 10^5 particles can require months or even years to complete, while an orbit-sampling Monte Carlo approach cannot adequately model the dynamics in a dense cluster core, particularly in the presence of many black holes. We have developed a new technique combining the precision of a direct N-body integration with the speed of a Monte Carlo approach. Our Rapid And Precisely Integrated Dynamics code, the RAPID code, statistically models interactions between neighboring stars and stellar binaries while integrating directly the orbits of stars or black holes in the cluster core. This allows us to accurately simulate the dynamics of the black holes in a realistic globular cluster environment without the burdensome N^2 scaling of a full N-body integration. We compare RAPID models of idealized globular clusters to identical models from the direct N-body and Monte Carlo methods. Our tests show that RAPID can reproduce the half-mass radii, core radii, black hole ejection rates, and binary properties of the direct N-body models far more accurately than a standard Monte Carlo integration while remaining significantly faster than a full N-body integration. With this technique, it will be possible to create more realistic models of Milky Way globular clusters with sufficient rapidity to explore the full parameter space of dense stellar clusters.

We present the first detailed comparison between million-body globular cluster simulations computed with a H\'enon-type Monte Carlo code, CMC, and a direct $N$-body code, NBODY6++GPU. Both simulations start from an identical cluster model with $10^6$ particles, and include all of the relevant physics needed to treat the system in a highly realistic way. With the two codes "frozen" (no fine-tuning of any free parameters or internal algorithms of the codes) we find excellent agreement in the overall evolution of the two models. Furthermore, we find that in both models, large numbers of stellar-mass black holes (> 1000) are retained for 12 Gyr. Thus, the very accurate direct $N$-body approach confirms recent predictions that black holes can be retained in present-day, old globular clusters. We find only minor disagreements between the two models and attribute these to the small-$N$ dynamics driving the evolution of the cluster core for which the Monte Carlo assumptions are less ideal. Based on the overwhelming general agreement between the two models computed using these vastly different techniques, we conclude that our Monte Carlo approach, which is more approximate, but dramatically faster compared to the direct $N$-body, is capable of producing a very accurate description of the long-term evolution of massive globular clusters even when the clusters contain large populations of stellar-mass black holes.

The predicted rate of binary black hole mergers from galactic fields can vary over several orders of magnitude and is extremely sensitive to the assumptions of stellar evolution. But in dense stellar environments such as globular clusters, binary black holes form by well-understood gravitational interactions. In this letter, we study the formation of black hole binaries in an extensive collection of realistic globular cluster models. By comparing these models to observed Milky Way and extragalactic globular clusters, we find that the mergers of dynamically-formed binaries could be detected at a rate of ~100 per year, potentially dominating the binary black hole merger rate. We also find that a majority of cluster-formed binaries are more massive than their field-formed counterparts, suggesting that Advanced LIGO could identify certain binaries as originating from dense stellar environments.

Hierarchical triple-star systems are expected to form frequently via close binary-binary encounters in the dense cores of globular clusters. In a sufficiently inclined triple, gravitational interactions between the inner and outer binary can cause large-amplitude oscillations in the eccentricity of the inner orbit ("Lidov-Kozai cycles"), which can lead to a collision and merger of the two inner components. In this paper we use Monte Carlo models of dense star clusters to identify all triple systems formed dynamically and we compute their evolution using a highly accurate three-body integrator which incorporates relativistic and tidal effects. We find that a large fraction of these triples evolve through a non-secular dynamical phase which can drive the inner binary to higher eccentricities than predicted by the standard secular perturbation theory (even including octupole-order terms). We place constraints on the importance of Lidov-Kozai-induced mergers for producing: (i) gravitational wave sources detectable by Advanced LIGO (aLIGO), for triples with an inner pair of stellar black holes; and (ii) blue straggler stars, for triples with main-sequence-star components. We find a realistic aLIGO detection rate of black hole mergers due to the Lidov-Kozai mechanism of 1yr^-1, with about 20% of these having a finite eccentricity when they first chirp into the aLIGO frequency band. While rare, these events are likely to dominate among eccentric compact object inspirals that are potentially detectable by aLIGO. For blue stragglers, we find that the Lidov-Kozai mechanism can contribute only up to ~10% of their total numbers in globular clusters.

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 $\sim 20 - 100\, M_\odot$. 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\times10^6$ stars. In almost all models we find that significant numbers of black holes (up to $\sim10^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 MSun, 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.