Dynamical Evolution

5 Dynamical Evolution

The evolution of the compact binary population in a globular cluster depends on the interplay between the global evolution of the cluster as described in Section 2.3, the few-body interactions described in Section 4.3, and the binary stellar evolution described in Section 4.1. Simulations that attempt to describe the compact binary population of a globular cluster must be able to couple these processes and describe their mutual interaction.

The evolution of globular clusters is primarily governed by point-mass gravitational attraction between individual stars and thus their overall structure can be described in terms of classical Newtonian N-body dynamics. Other processes, however, play a major role. Stellar evolution obviously affects the number and properties of compact objects and compact binaries so describing it accurately is important for producing the correct population of compact objects. Stellar evolution can also result in significant mass-loss from individual stars and this mass-loss changes the mass of the cluster enough to affect its dynamical evolution. Consequently the global dynamical evolution of a globular cluster is coupled to nuclear burning in individual stars. This is particularly true in the first 1 – 2 Gyrs of cluster evolution where massive stars are still evolving and where significant mass-loss can occur on timescales of 10 – 100 Myrs. Indeed, since this timescale is significantly shorter that the relaxation timescale, simulating the dynamical response of the cluster to this mass loss can be very numerically challenging. Stellar evolution also affects the orbital parameters of the binary population and, as we have seen in Section 2.3, binaries can be an energy source for the cluster and can slow down or halt core collapse. Thus the details of both binary evolution and few-body interactions can affect the global dynamics of a globular cluster. Finally, the loop is closed by the fact that the global dynamics of a globular cluster can affect the properties of the binary population by bringing stars and binaries close enough that they can interact and perturb each other. We are interested in the binary population that is the end result of this complex interplay of different processes operating on very different spatial and temporal scales. Simulating star cluster dynamics is a problem that has a long history and has promoted contact between various branches of astrophysics. The MODEST (MOdeling DEnse STellar systems) collaboration, a collection of various working groups from different backgrounds, maintains a website that provides some up-to date information about efforts to model the dynamical and stellar evolution of star clusters [340 ].

5.1 Star cluster simulation methods

Broadly speaking, there are three approaches that can be used to simulate the dynamical evolution of star clusters. The first, direct N-body integration, treats the dynamical evolution of a star cluster by numerically solving Newton’s equations. The second, the distribution function method, represents the the star cluster using a 6N dimensional distribution function and then describes its evolution using a collisional form of the Boltzmann equation. Both of these methods have their strengths and weaknesses. Several authors have compared these techniques and in general find that they agree well [198 , 262

, 201 , 463

]. The summary of the MODEST-2 meeting gives a fairly recent review of efforts to implement these techniques and use them in simulations [439 ]. The third method, the encounter rate technique, does not self-consistently simulate the evolution of star clusters but uses a static cluster model combined with individual scattering experiments to estimate the effect of dense stellar systems on a specific population. In the next three sections we outline these techniques. We then conclude the chapter by summarizing some of the most recent results about compact binary populations in star clusters derived using these methods.

5.1.1 Direct N-body integration

Direct N-body simulations are conceptually the most straightforward method of simulating star clusters. Each star is explicitly represented by a massive (non-test) particle and the gravitational interactions are calculated by numerically integrating the 3N coupled differential equations of motion in classical Newtonian gravity:

Here, $mj$ is the mass of particle $j$ and $⃗ri,j$ are the positions of particles $i$ and $j$ respectively. A great deal of general information on direct N-body methods can be found in the review article by Spurzem [444

] and Aarseth’s book on N-body simulations [4

There are currently two major families of codes used for direct N-body simulations: the Kira/STARLAB/AMUSE environment and the NBODYX series of codes. Kira integrates the N-body equations of motion using a 4th-order Hermite predictor-corrector scheme [317 ]. Kira is available as part of the star cluster evolution code STARLAB [238 ], a code that also includes a stellar and binary evolution package, SeBa, based on the analytic tracks of Hurley et al. [230 , 233 ], and various other routines for generating initial conditions and dealing with tidal processes. Kira is also available in the more recent star cluster simulation environment AMUSE (Astrophysical Multipurpose Software Environment) [384 ]. AMUSE is designed to be highly modular and allows the user a choice between various integrators, stellar evolution routines and prescriptions for dealing with stellar collisions. The NBODYX series of codes are a long-standing feature of the stellar dynamics community and have been in constant development by Aarseth since the late 1960s [1 ]. These codes employ the same basic 4th-order Hermite predictor-corrector scheme as Kira and include (in some versions) stellar evolution based on the same Hurley et al. analytic tracks. The NBODYX codes, however, treat binary evolution rather differently and are highly optimized for computational efficiency. Aarseth has produced two excellent reviews of the general properties and development of the NBODYX codes [3 , 2 ] and includes further details in his book on N-body simulations [4 ]. A fairly complete summary of common N-body codes and related applications can also be found on the NEMO website [353 ].

The primary advantage of direct N-body simulations is the small number of simplifying assumptions made in treating the dynamical interactions. All gravitational interactions are explicitly integrated so the direct N-body method can, in principle, resolve all of the microphysics involved in the evolution of globular clusters. In this sense direct N-body simulations represent star clusters “as they really are”. Furthermore, since all trajectories involved in any interaction are known, it is easy to study each interaction in detail. Finally, because direct N-body simulations produce star-by-star models of globular clusters, it is relatively straightforward to couple them with additional stellar physics such as stellar evolution, stellar equations of state, stellar collision models and so on. Therefore, within the limits of numerical error, the results of direct N-body simulations are very reliable and represent the gold standard in globular cluster modelling.

Direct N-body simulations are, however, very computationally expensive. Over a relaxation time, the computational cost of a direct N-body simulation intrinsically scales as $N 3−4$ [199 ]. Two factors of N come from calculating the pairwise interactions, another factor comes from the increasing number of interactions per crossing time as N increases, and a final factor of $N ∕log N$ comes from the increasing length of $trlx$ compared to the crossing timescale (see Eq. 5 ). Another problem is the very large range of scales present in globular clusters [236 , 199 ]. Distances, for instance, can range from a few tens of kilometers (neutron star binaries) to several tens of parsecs (the size of the cluster itself). The timescales can have an even larger range, from several milliseconds (millisecond pulsar binaries) to several Gigayears (the relaxation times of very large globular clusters). In the most extreme cases it may be necessary to know positions to 1 part in 10¹⁴ while a cluster may need to be advanced on timescales 10²⁰ times shorter than its lifetime. Furthermore, N-body systems are known to be chaotic [335 , 336 , 169 ]. This means that the integration schemes must be very accurate (4th-order or better). Consequently they are very computationally expensive. Overall this means that with increasing N direct N-body simulations rapidly become computationally intractable. Pure N-body simulations are scale-free and hence it is, in principle, possible to scale up the results of a small-N simulation [35 ]. However, this is not possible once other physical processes, such as stellar evolution, are included. Thus the need for large-N simulations is unavoidable.

There are several methods that can be used to mitigate these problems. The Hermite scheme, for instance, requires only the first and second derivatives of acceleration, even though it is a 4th-order method. Since dynamical evolution proceeds more quickly in regions of high stellar density, the stars in cluster halos do not need to have their phase space parameters updated as frequently as those in cluster cores. Both STARLAB and NBODYX include individual adaptive timesteps to ensure that during each force calculation, only the necessary particles have their parameters updated. Some of the NBODYX codes also include a neighbour scheme so that not all pairwise forces need to be explicitly calculated at each timestep. Binaries tight enough that their internal evolution is largely unaffected by external forces from the cluster can also be removed from the global dynamical evolution and replaced with a single center-of-mass particle. The evolution of the binary itself can then be followed using an efficient special-purpose code. STARLAB calculates the evolution of binaries using a semi-analytic Keplarian two-body code that admits perturbations from the cluster environment. The NBODYX codes use a regularization technique to achieve the same effect with the additional advantage of removing the singularity created by very close binaries. Both of these methods avoid the global timestep of the star cluster being reduced to that of its tightest binary.

Attempts have been made to create parallel direct N-body codes with NBODY6++ [444 ] being the most successful example. A major problem with this type of parallization is that each particle in an N-body calculation formally requires information about every other particle for each force calculation. This means that parallel N-body codes are communication dominated and do not scale well with the number of processors. Another approach has been to use hardware acceleration. The GRAPE (GRavity PipE), invented by Makino [318 ], is a special-purpose computer designed to rapidly calculate the $1∕r2$ force necessary for direct N-body simulations. The current generation, the GRAPE-6, was reported to have a peak speed of 64 T-flops in 2002 [460 ]. A single GRAPE 6 board, the GRAPE-6A, is also available as PCI card [142 ]. The next generation GRAPE, the GRAPE-DR, is predicted to have a peak speed of $∼$ 2 P-flops [316 ]. It is not yet clear when (or if) the GRAPE-DR will become commercially available. Another very promising possibility is hardware acceleration using GPUs (Graphical Processing Units) [446 ]. GPUs provide a similar service to the GRAPE while having the advantages of being available off-the-shelf and commercially supported. AMUSE, NBODY6 and NBODY6++ have been modified to make use of GPUs.

Despite developments in both hardware and software, current direct N-body simulations are practically relegated to $N ≲ 105$ , about an order of magnitude lower than the number of stars in a typical globular cluster. Furthermore, such simulations can take several months of real-time computation to reach the $∼$ 10 Gyr age estimated for Galactic globular clusters and this makes large parameter-space studies impossible. For this reason more approximate methods employing distribution functions are still widely used in stellar dynamics.

5.1.2 Distribution function methods

Some of the computational problems associated with direct N-body simulations can be avoided by describing globular clusters in terms of distribution functions rather than as ensembles of particles. The distribution function for stars of mass $m$ at time $t$ , $fm(⃗x, ⃗v,t)$ , is 6N-dimensional in position-velocity phase space. The number of stars in the phase space volume $3 3 d x d v$ is given by $3 3 dN = fmd x d v$ . This description requires that either the phase space volume $3 3 d x d v$ is small enough to be infinitesimal yet large enough to be statistically meaningful or, more commonly, that $fm$ be interpreted as the probability of finding a star of mass $m$ in the volume $d3x d3v$ at time $t$ . The time evolution of the cluster is modelled by calculating the time evolution of $fm$ . The effect of gravity is divided into two components. The first is a smoothed potential, $ϕ$ , that describes the effect of the overall gravitational field of the cluster. $ϕ(⃗x)$ is found by solving the Poisson equation, which can be written in terms of the distribution function as:

The second component is a collision term, $Γ [f]$ , that can be thought of as the effect of random close encounters between individual stars. In this picture $ϕ$ represents the smooth background gravitational field of the cluster while $Γ [f ]$ represents the granularity induced by the mass concentration in individual stars. The time evolution of $fm$ is calculated using a modified form of the Boltzmann Equation:

In stellar dynamics, Eq. 43

is usually called the Fokker–Planck equation.

Eq. 43 can be directly solved (numerically) if a specific expression for $Γ [f]$ can be found. One common choice is to consider a cluster in the weak scattering limit where distant two-body interactions drive dynamical evolution. This is essentially equivalent to assuming the cluster evolves only by two-body relaxation as described in Section 2.3. In this case $Γ [f]$ can be described by an expansion in the phase space coordinates ( $⃗x,⃗v ∈ ⃗w)$ and takes on the form [57 ]:

where $⟨ΔX ⟩$ is the expectation value of the phase space variable over some time $Δt$ . In practise it is not efficient to solve Eqs. 43

and 44

in Cartesian coordinates but an appropriate choice of $Δt$ can alleviate this problem. The key observation is that there are two relevant timescales for evolution in the Fokker–Planck approximation: the crossing time, $tcr$ , that governs changes in position, and the relaxation time, $trlx$ , that governs changes in energy. In the case where $tcr ≪ trlx$ changes in position are essentially periodic and at any given time the orbit of a star can be written in terms of its energy, $E$ , and the magnitude of its angular momentum vector, $J$ . By making a careful choice of $Δt$ and suitable assumptions about the symmetry of the potential and the velocity dispersion [77 , 443

] it is possible to describe the evolution of the phase space structure of the cluster solely in terms of the orbital parameters and Eqs. 43

and 44

can be written in terms of $E$ and $J$ alone. In this representation Eq. 43

is called the orbit-averaged Fokker–Planck equation.

Numerically solving the orbit-averaged Fokker–Planck equation is much less computationally intensive than a direct N-body simulation and can be used to model a full-sized globular cluster. Since these simulations are much faster (taking hours or days rather than months) they are also useful for performing parameter space studies. It is also possible, in principle, to improve the accuracy of the Fokker–Planck method by including higher-order terms in Eq. 44 [420 ]. The Fokker–Planck method can be generalized to include velocity anisotropies, tidal stripping, non-spherical potentials and cluster rotation [118 , 260 , 262 ].

The Fokker–Planck method also has some serious limitations. Each stellar mass in the cluster requires its own distribution function in Eq. 43 . For more than a few distribution functions the Fokker–Planck method becomes numerically intractable so the ability of the method to resolve a continuous spectrum of masses is limited. Since the Fokker–Planck method uses a distribution function, it does not automatically provide a star-by-star representation of a globular cluster. This means that it can be difficult to compare the results of a Fokker–Planck simulation with observations. It is also difficult to include stellar evolution in Fokker–Planck simulations because there are no individual stars to evolve. For the same reason Fokker–Planck methods cannot easily model strong few-body interactions, interactions that are essential for describing the evolution of the compact binary fraction in star clusters [144 ]. Therefore Fokker–Planck methods are of limited use in studying binary populations in star clusters.

Another approach to solving Eq. 43 that avoids some of the problems of the Fokker–Planck approach is the Monte Carlo method. This method was first developed by Hénon [210 , 209 , 211 ] and significantly improved by Stodółkiewicz [449 , 450 ]. It does not use the distribution function explicitly but rather represents the cluster as a ensemble of particles, just like a direct N-body simulation. Each particle is characterized by a mass, a distance from the cluster center, and a tangential and radial velocity. The underlying treatment of relaxation is the same as in the Fokker–Planck approach but is calculated on a particle-by-particle basis. At the beginning of a Monte Carlo timestep the code calculates the global potential based on the current mass and radial position of each particle. This potential is then used to generate a plane-rosette orbit for each particle, defined by $E$ and $J$ (the orbit-averaged approximation). Next each particle has its orbit perturbed by a weak encounter, the parameters of which are calculated between it an a neighbouring particle. The key to the method is that assuming weak scattering, this perturbation is statistically representative of all weak encounters the particle will experience over a time $Δt$ . Thus the results of the perturbation can be multiplied by an appropriately chosen factor to represent the total perturbation the particle receives over $Δt$ . Once the new orbits have been found, new positions can be randomly chosen from a time-weighted average along the orbital path, a new potential calculated, and the process can start again.

The Monte Carlo method is mathematically identical to the Fokker–Planck method but, since each star is individually represented, it is much easier to couple with additional physics. Each star in the simulation can have a different mass with no loss of efficiency so the method can resolve a continuous mass spectrum. Stellar evolution can be included in exactly the same way as a direct N-body code. It is also relatively easy to include strong few-body interactions. For each star or binary at each timestep a cross-section for a strong interaction with another star or binary can be calculated and the encounter resolved using either analytic prescriptions or explicit few-body integration. Thus the Monte Carlo method has many of the advantages of the direct N-body method while being much faster. Since there are a fixed number of operations per particle, the Monte Carlo method scales almost directly with N rather than the $N 3−4$ dependence of a direct N-body code. Thus a simulation of $6 N ∼ 10$ can be completed in a matter of days, rather than the months required by the direct N-body method. Monte Carlo codes have been shown to re-produce observations [197 , 158 , 159 ] and the results of direct N-body codes [155 , 128 , 72 ] quite well, although some discrepancies remain [463 ].

The Monte Carlo method does have some shortcomings. It is limited by the weak scattering approximation and is not always reliable in the core region where interactions can be very strong. It is also strongly limited by what is assumed for the cross-sections for various types of few body interactions. If these are incorrect than the number and relative importance of single-single, binary-single and binary-binary interactions will be wrong. The only way to calibrate these interactions is by careful comparison with direct N-body simulations. The method is further limited by how accurately the outcomes of few-body interactions themselves are calculated. Analytic prescriptions for the outcomes are fast but cannot cover all initial conditions accurately while direct few-body integration is more general but imposes a computational cost. While not a fundamental requirement of the method, current Monte Carlo codes are one-dimensional (in that the position of each star is defined only by its radial position) and thus they are limited to spherical symmetry. In particular this means Monte Carlo codes currently cannot model rotating clusters. A rotating Monte Carlo code is possible in principle but would require a significant amount of development.

Despite these limitations Monte Carlo methods have proved to be quite effective for efficiently modelling star clusters and there are several version currently available, in particular those of Freitag [137 ], the Northwestern Group [253 , 252 , 130 , 72 ] and Giersz [153 , 154 , 155 , 159 , 160 ]. The Freitag code is focused on simulations of star clusters with central massive objects but the Northwestern and Giersz codes are suitable for general cluster dynamics and are roughly comparable in capabilities.

Another option for solving Eq. 43 is to use an analogy between between a globular cluster and a self-gravitating gaseous sphere [298 , 161 , 445 ]. These models are solved by taking moments of the Boltzmann equation and, like the Fokker–Planck method, can have their accuracy increased by including higher-order moments [420 ]. The term “anisotropic gas model” is often used to describe these methods due to the analogy between hydrodynamical models of stellar evolution and due to the fact that they can admit an anisotropic velocity distribution. An accessible description of some of the mathematical details of the method are given by Amaro-Seoane [15 , 12 ]). The anisotropic gas models still operate in the continuum limit so they have the same difficulty dealing with binary populations, strong few-body interactions and stellar evolution as the Fokker–Planck method. For this reason they are of limited use for simulating binary populations in globular clusters. However the method is applicable to the stellar dynamics of very large-N systems such as galactic centers and has been particularly useful for modelling systems containing massive central compact objects such as supermassive black holes (e.g. [15 ]). An attempt to improve the treatment of binaries in the anisotropic gas model has resulted in the so-called “hybrid” code where single stars are treated by the gaseous moment equations while strong few-body encounters are treated using a Monte Carlo method [162 , 163 ]. So far this method has not been used to examine compact binary populations in star clusters.

5.1.3 Encounter rate techniques

Encounter rate techniques are another method that has proven very useful for understanding the effect of dense stellar environments on binary populations. These methods do not deal with the global evolution of globular clusters but rather deal with the interactions between an environment and a particular stellar species. Normally a stellar density distribution is established for a globular cluster (or its core) and then the number of interactions per unit time is calculated for a specific object using cross-sections. The density distribution, interactions and outcomes can be calculated either analytically or numerically, depending on the complexity of the system in question. These methods are limited but are very useful for estimating which processes will be the most significant for specific stellar populations. There is no “general method” for encounter rate techniques but we will discuss some specific implementations in Section 5.2.

5.2 Results from models of globular clusters

Over the past two decades, there has been extensive work modeling the dynamical evolution of globular clusters. Published studies have employed encounter rate techniques in static cluster backgrounds, Fokker–Planck and Monte Carlo models and direct N-body simulations – e.g., [51 , 54 , 334 , 372 , 383 , 402 , 405 , 421 , 429 , 438 , 458 , 95 , 93 , 94 , 96 , 289 , 163 , 155 , 159 , 253 , 252 , 130 , 244 , 132 , 131 , 72 , 387 , 390 , 386 , 315 , 232 , 231 , 229 ]. These models can, if properly interpreted, provide a wealth of information about stellar populations in globular clusters. These models, however, have limitations that must be understood. The dynamical approximations employed coupled with an imperfect understanding of stellar and binary evolution make it impossible to compare the results of these simulations directly to observations. The problem is exacerbated by the fact that many authors are more concerned with the development of numerical methods than on producing astrophysical results. Furthermore, few of these studies focus on (and some do not even model) the compact binary population. Those that do often focus on the effect that the compact population will have on the cluster dynamics rather than the effect the cluster dynamics will have on the compact binary population. Nevertheless, many published simulations provide useful insights into the behaviour of compact binary populations in dense stellar systems. Here we report on these results for binaries where the primary is a WD, NS or BH, trying to give both a historical overview as well as details of the latest results. When reading past papers, attention must be paid to terminology. A neutron star binary, for instance, can imply a binary composed of two neutron stars or simply a binary where a neutron star is paired with any type of secondary depending on the authors preference. Here we attempt to clearly indicate the nature of the secondary in order to resolve any ambiguities.

5.2.1 Binaries with white dwarf primaries

There are two kinds of binaries with white dwarf primaries that have attracted particular attention from modelers: cataclysmic variables (CVs) and double white-dwarf binaries (DWDs). As described in Section 3.1, CVs consist of white dwarfs accreting matter from a companion, normally either a dwarf star or another white dwarf. CVs with a dwarf star companion are not, strictly speaking, relativistic binaries as defined in the introduction. Unlike many double-degenerate binaries, however, CVs can be observed in the electromagnetic spectrum and they are very useful for understanding the behaviour of compact objects in globular clusters. CVs with a helium white-dwarf companion are known as AM CVn systems and are a type of DWD. DWDs, which as the name suggests are binaries of two white dwarfs, are relativistic binaries but not all are AM CVns and are not always visible in electromagnetic radiation. The merger of DWDs may, however, be observed as a SNe Ia and at wider separations they may be observed as low-frequency gravitational wave sources. Thus they too have attracted attention from globular cluster modelers.

White dwarfs are the lowest-mass compact stellar remnants and are not necessarily significantly more massive than the main sequence turn-off of present-day globular clusters. Because of this they may not experience significant mass segregation or be as centrally concentrated as higher-mass stellar remnants. Indeed for much of the life of a star cluster, white dwarfs may be lighter than the average stellar mass in the cluster and could be pushed to the outskirts. This means that white dwarfs are not necessarily found in regions of high stellar density so they may actually have a low probability of interaction with other stars and binaries. Furthermore, when interactions do occur their lower masses mean they may not be clearly preferred for exchange into binaries. Due to these uncertainties modeling is crucial to understanding their behaviour in globular clusters.

Lee (1987) [289 ] is an example of an early numerical study of the compact binary population in globular clusters, focusing on the effect of the binaries on the cluster dynamics. It is based on Fokker–Planck models using two different stellar mass components to represent main sequence stars and compact objects. It demonstrated that compact objects have a strong effect on the cluster core, implying they experience many interactions and highlighting the importance of the cluster environment for compact binaries.

One of the earliest studies of CVs in globular clusters is found in Di Stefano and Rappaport (1994) [103 ] who performed a large series of simulations using an encounter-rate technique. Specifically, they performed numerical two-body scattering experiments in static cluster backgrounds coupled with binary stellar evolution prescriptions in order to investigate the tidal capture and subsequent evolution of CVs. The cluster parameters were chosen to match the current states of 47 Tuc and $ω$ Cen. Since these clusters are frequently compared in studies of compact binaries, some of their properties are given in Table 4 (more details can be found in the Harris Catalogue of globular clusters [193 ]). The primary differences between these clusters are that $ω$ Cen is significantly more massive that 47 Tuc while 47 Tuc is much more concentrated and dynamically older than $ω$ Cen.

Table 4: Structural parameters for 47 Tuc and $ω$ Cen taken from Davies and Benz [95

]. $Mtot$ is the total mass, $ρ 0$ the central density, $r c$ the core radius, $r t$ the tidal radius and $t rh$ the current estimated half-mass relaxation time.

Cluster	$Mtot$ [ $M ⊙$ ]	$ρ0$ [ $M ⊙$ /pc]	$rc$ [pc]	$rt$ [pc]	$trh$ [Gyr]
$ω$ Cen	$3.34 × 106$	$4.0 × 103$	3.4	74.0	16.0
47 Tuc	$1.06 × 106$	$5.1 × 104$	0.59	69.0	8.7

The simulations of Di Stefano and Rappaport produced 150 – 200 CVs in 47 Tuc, $∼$ 45 of which had an accretion luminosity of $> 1033 erg s−1$ and would be visible at the current epoch. For $ω$ Cen they predicted 100 – 150 and $∼$ 20 respectively. They claimed a factor of $≲$ 10 enhancement in the number of CVs over the field population. The authors point out both that this number of CVs was barely within the upper limits on the number of CVs in clusters known observationally at the time and that the number could be even higher when binary-single and binary-binary interactions were taken into account. Observations of CVs in globular clusters were in their infancy at the time and many more have since been discovered. However the over-production of CVs remained a problem in later works.

The next studies to touch on CVs in globular clusters were a series of papers by Davies and collaborators in 1995 – 1997 [95 , 93 , 94 ] concerned with exotic binaries produced by dynamical interactions. These were also encounter rate models but included binary-single interactions as well as tidal capture. Davies and collaborators [95 , 93 ] reported many single-single interactions between white dwarfs and main sequence stars but found that the majority of strong encounters led to mergers rather than CVs. Overall Davies predicted more CVs in 47 Tuc (smaller, concentrated, dynamically older) than than in $ω$ Cen (larger, less concentrated, dynamically younger) but more merged systems than CVs in both. They also predicted a reasonable number of DWD mergers in both clusters, most of which were above a Chandrasekhar mass and could lead to SNe Ia.

In 1997 Davies [94 ] produced another encounter rate study concentrating on the effect a cluster environment would have on a population of binaries that would become CVs due to binary stellar evolution in isolation. This study observed that disruption of binaries in a dense stellar environment is important and can reduce the number of CVs compared to the field. Davies found that most CV progenitors in cluster cores with number densities $> 103$ are disrupted during dynamical interactions. Thus, if dynamical interactions form no new CVs, the cores of galactic globular clusters will be depleted in CVs. CV progenitors in the low-density halos of globular clusters may survive to become CVs if the relaxation time of the cluster is long enough that they do not mass-segregate to the core. Davies proposed that comparing the CV population in the cores and halos of globular clusters could be used as a probe of their dynamical properties. A cluster with a core depleted in CVs compared to the halo will not have experienced much dynamical binary formation while a cluster with a large number of CVs in the core must have experienced dynamical binary formation at some point in its history.

The first detailed study of DWDs in dense stellar environments became available with the direct N-body simulations of Shara and Hurley (2002) [429 ]. They were specifically interested in the galactic rate of SNe Ia and if globular clusters could significantly enhance it. They used the NBODY4 code to perform four 22,000 particle simulations with 10% primordial binaries, two different metallicities and detailed stellar evolution prescriptions. These simulations were in the open rather than globular cluster regime and were terminated after only 4.5 Gyr when they had lost more than 75% of their initial mass. They were run on 16 chip GRAPE-6 boards and took approximately five days each to complete.

Shara and Hurley found a factor of 2 – 3 enhancement in the number of DWDs in a clustered environment compared to a similar field population. Exchange interactions seemed to be particularly important since of the 135 DWDs produced between all four simulations, 93 formed by exchanges. The enhancement in SNe Ia progenitors, DWDs with a total mass above a Chandrasekhar mass and a gravitational wave inspiral timescale shorter than a Hubble time, was even greater. They found $∼$ 4 SNe Ia progenitors per 2000 binaries, a factor of $∼$ 15 more than in the field. Both exchange interactions and hardening during fly-bys were responsible for the enhancement. Any DWD that went through a common envelope phase was assumed to be circularized and they found that later encounters did not induce significant eccentricity in the binaries. Thus they concluded that most DWDs in globular clusters will be circular. It is worth noting, however, that most of these DWDs formed promptly, at less than 1 Gyr of cluster age, and few had inspiral timescales longer than 1 Gyr. Therefore the applicability of these results to globular clusters is unclear. They found that neither normal CVs nor AM CVns were significantly enhanced over the field population. However, those in the clusters were often the result of exchange interactions and their properties were different from the AM CVns in the field. Overall this implies that most CVs in star clusters are the product of dynamics and that pure binary stellar evolution channels for CV and AM CVn creation are shut down in clustered environments.

Shara and Hurley (2006) [430 ] continued this work using a single 100,000 simulation with 5000 primordial binaries. This simulation took six months to run to an age of 20 Gyr using a 32-chip GRAPE-6 board. The focus of this simulation was on normal CVs, 15 of which were produced during the course of the simulation and three of which would be visible at the current epoch. A comparison between these simulations and the more statistical methods of Davies et al. illustrate the relative strengths of these methods. After six months the direct N-body simulations of Shara and Hurley produced only 15 objects for analysis, too few for solid statistics. Additional simulations for more objects or for verification purposes were too computationally expensive. By contrast the simulations of Davies et al. produced many more objects and multiple realizations of each cluster model were possible. This allowed for greater confidence in the statistical significance of the results. However, the N-body simulations of Shara and Hurley made no assumptions about the cross-sections for various types of interactions and they naturally included all dynamical processes. Without comparison to N-body simulations there is no way to tell if there were systematic errors in the cross-sections used by Davies et al. (or any other statistical method) or if important interactions were excluded. Both types of simulations are necessary for investigating the compact binary population in globular clusters. Direct N-body simulations provide a complete understanding of the physical processes involved while encounter rate techniques or distribution function methods can be used to provide the multiple realizations for the statistical verification and parameter space coverage that is too computationally expensive to produce using direct N-body methods.

Shara and Hurley found that dynamical effects pushed several of the progenitors to the active CV state in a shorter time than would have occurred in the field. Indeed one of these binaries would have merged without becoming a CV if not for an orbital perturbation. Furthermore the simulation produced two binaries formed by exchange interactions that clearly have no analogs in the field. The exchange CVs were, on average, more massive than CVs in the field. They confirmed that many CV progenitors were destroyed by dynamical interactions and showed that a field population of stars with similar properties would produce 17 CVs over 20 Gyr, nine of which would be visible at the current epoch. Thus Shara and Hurley predicted $∼$ 2 – 3 times fewer CVs in globular clusters than in the field.

This simulation, although primarily focused on CVs, included some discussion of the DWD population and generally confirmed the picture presented in Shara and Hurley 2002. In particular short-period ( $P ≤ 1$ day, comparable to the largest period of SNe Ia progenitors in Shara and Hurley 2002) DWDs were enhanced by a factor of 5 – 6 over the field population. Although this was somewhat less than the factor of $∼$ 15 previously reported it is still much larger than the enhancement found for in CVs. Shara and Hurley proposed several reasons for more efficient DWD production in star clusters. In a DWD both components must be sufficiently massive to become white dwarfs within a Hubble time while only one does in a standard CV. Thus CVs are a much more common outcome of pure binary stellar evolution than are DWDs. DWDs are also, on average, more massive than CVs and thus tend to mass-segregate further towards the core and undergo more interactions. Finally, since each component of a DWD is likely to be more massive than the main sequence star in a CV and since exchanges tend to place massive objects in binaries, a DWD is a more likely outcome of an exchange interaction between a main sequence star and two white dwarfs than is a CV.

One of the most complete projects on the compact binary population in the recent literature is found in a series of papers by Ivanova and collaborators [405 , 244 , 248 , 247 , 246 , 245 ]. These papers simulate the binary – and particularly the compact binary – population of star clusters using a combination of Monte Carlo simulations and well-developed encounter rate techniques coupled with stellar and binary evolution algorithms tailored towards compact objects [230 , 233 , 47 ]. These models generally assume a simplified two-zone core-halo model for a globular cluster. The central number density of stars, the central 1-D RMS velocity dispersion and $trh$ are all taken to be fixed quantities. The stars are then distributed with the correct core density with velocities drawn from a Maxwellian choose to give an appropriate 3-D velocity dispersion in the core. Stars move from the core to the halo stochastically based on probabilities calculated from estimates of the mass-segregation timescale. This type of model can be used to produce a number density, mass distribution and velocity dispersion for the core and halo at given intervals. From these it is possible to calculate cross sections for stars and binaries to interact. Whether and what type of interaction occurs is determined by Monte Carlo sampling while the interactions themselves are resolved by explicit few-body integration.

Initially, Ivanova et al. (2005) [244 ] considered the effect of cluster dynamics on the number of all types of binaries globular clusters and showed that dynamical disruption will significantly reduce the binary fraction for reasonable initial conditions. Indeed they showed that a globular cluster where 100% of stars are initially binaries will normally have a core binary fraction of only $∼$ 10% after several Gyrs. This implies that globular clusters may have had high binary fractions when young and thus binary-single and binary-binary interactions may be more important than previously thought. They also demonstrated both the binary burning that supports the cluster against core collapse and a build-up of short-period white-dwarf binaries in the core, as shown in Figure 12 .

Figure 12: Binary period distributions from the Monte Carlo simulation of binary fraction evolution in 47 Tuc. The bottom panel indicates the period distribution for binaries containing at least one white dwarf. $Nb$ is the total number of binaries and $Nb,p$ is the number of binaries per bin. Image reproduced by permission from Ivanova et al. [244

The models were improved by Ivanova et al. (2006) [247 ] in order to investigate the white-dwarf binary population in globular clusters. Stellar evolution and tidal capture were updated and single-single, binary-single and binary-binary interactions were all included. They confirmed that CVs are not highly over-produced in globular clusters relative to the field population but found that CV production is moderately efficient in cluster cores. In the cores they predicted an enhancement over the field population of a factor of 2 – 3. As with Shara and Hurley, the globular cluster CVs had different properties than field CVs. In particular the globular cluster CVs were slightly brighter and had a larger X-ray-to-optical flux. Ivanova et al. attributed this to slightly higher-mass white dwarfs in the cluster CVs, just as found by Shara and Hurley. The more massive white dwarfs led to globular cluster CVs that extracted more energy at the same mass transfer rate than field CVs and had slightly higher magnetic fields. This was enough to account for the observed X-ray excess.

Ivanova et al. shed more light on the dynamical processes that form CVs than previous encounter rate studies because they included most of the likely formation mechanisms in the same simulation. Their results showed that some 60 – 70% of CVs in the core form because of dynamical interactions while most of the CVs in the halo originate from primordial binaries. Overall they found that only $∼$ 25% of CVs formed in the clusters would have become CVs in the field. They concluded that tidal capture of main sequence stars is normally unimportant for CVs but that collisions between red giants and main sequence stars – a situation where the envelope of the red giant can be stripped away leaving its core as white dwarf bound to the main sequence star – could provide up to 15% the total number of CVs in the clusters. Exchange interactions were the most efficient channel, forming $∼$ 32% of all CVs in their simulations, while collisions during binary-single or binary-binary encounters produced another $∼$ 13%. The remaining CVs were produced by primordial binaries but $∼$ 20% of these experienced hardening encounters in fly-by interactions and would not have become CVs in the field. Finally, Ivanova et al. found that 60% of their CVs did not form solely by common envelope evolution, the most common channel for field CV production and thus highlighted the importance of cluster dynamics for the CV population. Ultimately Ivanova et al. predicted 1 detectable CV per $1000 – 2000 M ⊙$ of mass in Galactic globular clusters. They specifically predicted 35 – 40 CVs in 47 Tuc, of the same order as the 22 observed there at the time. The discrepancy could be due to stochasticity, some inherent deficiencies in the modeling or because not all CVs in 47 Tuc have been observed. Observations of globular clusters in the 2000s, particularly in the UV [205 , 105 , 306 ] have indicated a larger number of CVs than previously thought. The combination of fewer CVs in current models and a larger number of detected CVs has gone a long way towards solving the “CV problem” in globular clusters.

Ivanova et al. did not find a strong enhancement in DWDs compared to their reference field population, unlike Shara and Hurley who found an order-of-magnitude enhancement. This could be “‘real” in the sense that their simulations model larger systems than Shara and Hurely’s work and could be taken as evidence that, while DWDs are enhanced in open clusters, they are not in globular clusters. It could also be a result of the different numerical methods. The two-zone model of Ivanova et al., for example, may not be able to represent mass segregation as accurately as the direct N-body models of Shara and Hurley and thus the effect of interactions in the dense core could be underestimated. It is also possible that the encounter prescriptions were missing details that are important for DWD creation. Although not order-of-magnitude, Ivanova et al. did find a modest enhancement in the number of SNe Ia progenitors over the field population. They predicted that this will not contribute significantly to SNe Ia from Milky Way-type galaxies because only a very small fraction of their mass is in globular clusters. Elliptical galaxies, however, tend both to have a larger fraction of their mass in globular clusters and to have significantly older stellar populations. Since SNe Ia are usually associated with younger stellar populations the contribution of globular clusters to the SNe Ia rate in elliptical galaxies could be important. Ivanova et al. pointed out that DWD mergers and accretion induced collapse of WDs can produce neutrons stars and may do so without imparting a natal kick. Given their assumptions about common envelope evolution Ivanova et al. found that the merger of the DWDs produced a number of neutron stars comparable to the number formed by massive stellar evolution with natal kicks and then retained. Indeed they found that merger-induced collapse of DWDs led to an over-production of neutron stars in globular clusters. They argued that either merger-induced SNe Ia must not lead to neutron star production in the majority of cases or that they over-estimated the energy dissipation efficiency of common envelope evolution.

Willems et al. 2007 [487 ] performed an extended analysis of the simulations of Ivanova et al. 2006 in order to estimate the number of eccentric DWDs in the low frequency gravitational wave band accessible to space-based gravitational wave detectors. Taking the DWD eccentricities from the 10th to 13th Gyr of the Ivanova et al. simulations, Willems et al. generated the probability distribution for eccentric DWDs in Galactic globular clusters shown in Figure 13 . These results indicate that there may be a small but non-negligible population of eccentric DWDs in the Galactic globular cluster system (see Table 5). They will have no analogs in the field because field DWDs are expected to be circularized by common envelope evolution. Thus a detection of an eccentric DWD by an appropriate detector is both probable and a clear indication of a dynamical origin.

Figure 13: Statistical properties of eccentric DWDs ( $e > 0.01$ ) in the LISA band ( $P < 5000 s$ ) for three different cluster models. Hatched histograms are for the entire eccentric DWD population while shaded histograms arefore the DWDs with at least two harmonics individually resolvable by LISA with an $S∕N > 8$ (5 year mission, 10 kpc distance). The PDFs are normalized to unity. From left to right: frequency, eccentricity, total mass, mass ratio and lifetime before merger. Plot reproduced by permission from Willems et al. 2007 [487

], copyright by IOP.

The consensus that emerges from the simulations is that dynamical destruction of potential CVs and dynamical formation of new CVs are equally important in globular clusters. Therefore the number of CVs per unit mass is unlikely to be much higher in globular clusters than it is in the field. Globular cluster CVs may, however, be slightly more massive than field CVs and thus have slightly higher X-ray fluxes. Overall current simulations seem to be able to re-produce the observed number of CVs in clusters reasonably well but still tend to slightly over-predict the population. However, the observational trend is towards finding more CVs in globular clusters so this problem may soon disappear.

Table 5: The mean number of eccentric binaries per cluster ( $⟨N ⟩$ ) and the Poisson probability ( $P N$ ) that at least $N$ eccentric DWDs are present at any given time. Data taken from Willems et al. 2007 [487

GC Model	$⟨N ⟩$	$P1$	$P2$	$P3$
Typical	0.39	0.32	0.06	0.01
47 Tuc	1.90	0.84	0.55	0.29
Terzan 5	1.60	0.80	0.48	0.22

It seems likely that the DWD population in globular clusters is larger per unit stellar mass than that in the field however it is not clear if the enhancement is by factors of ten or only factors of a few. The disagreement between various simulations could be due to differing numerical methods or because scaling the DWD population from open to globular clusters is not straightforward. Systematic comparison of the different simulation methods will be necessary in order to resolve this unambiguously and ultimately large N-body simulations may be the only way to understand white-dwarf binaries in globular clusters in detail.

5.2.2 Binaries with neutron star primaries

Neutron star binaries were the first binaries with degenerate companions to be specifically investigated in globular clusters. LMXBs, neutron stars accreting from a companion as described in Section 3.2, have attracted particular interest since results from 1975 [256 , 74 ] gave indications that they are over-abundant in globular clusters. These papers immediately suggested this was a result of dynamical encounters and by the end of the 1970s a debate over which type of encounter, whether tidal capture by MS stars [121 ], collisions with RG stars [455 ], exchange interactions [217 ] or three-body interactions [433 ], was already underway. Like CVs, LMXBs are not necessarily relativistic by our definition but are electromagnetically active and thus useful for calibration purposes. The ultracompact systems discussed in Section 3.2, often called UCXBs, are almost certainly double-degenerate and may consist of a neutron star accreting from a white dwarf.

Neutron stars are, on average, more massive than the main sequence turn-off of old globular clusters and also more massive than white dwarfs. As such they should both mass-segregate to cluster cores and be preferentially exchanged into binaries during interactions. Therefore globular clusters may produce a large number of double neutron star binaries (DNSs), the mergers of which will be excellent sources of gravitational radiation for ground-based detectors and may also be a source of short gamma-ray bursts [181 , 226 ]. However, neutron star production in globular clusters is by no means certain. Neutron stars are rare and thus primordial DNSs are not common. Furthermore large natal kicks can both break up primordial DNSs and remove neutron stars from globular clusters. Also, neutron stars are only a factor of two to three times more massive than the average cluster mass, leaving the extent of their mass segregation and their probability of exchange unclear.

Neutron stars are low-luminosity and, if not members of an X-ray binary, are normally detected as MSPs. Their presence was long-expected in globular clusters due to the large population of LMXBs that are thought to be MSP progenitors but were not confirmed until 1987 [304 ]. Since all MSPs are thought to need a spin-up phase during accretion in a binary system, even single MSPs can be used as a probe of binary evolution. As discussed in Section 3.3, these are now known to be quite common in globular clusters. Of more interest for gravitational wave detection are the many binaries containing MSPs and particularly the double MSP binaries (DMSPs) that have been detected in some globular clusters.

Originally tidal captures or collisions were thought to be the main paths for producing X-ray binaries. When the LMXB enhancement was first discovered, globular clusters were thought to have very few primordial binaries and thus would be have few exchange interactions. Since three-body interactions tend to produce soft binaries they were largely ignored. During the 1980s studies employing static cluster models and encounter rate estimates based either on analytic estimates (e.g., [240 , 469 , 23 , 472 ]) or few-body scattering experiments (e.g., [277 , 278 ]) were able to roughly re-produce the X-ray binary population in globular clusters using only tidal capture and collisions. They also indicated that hardening encounters might be important in bringing dynamically formed binaries into the mass-transferring, and thus the X-ray, regime. Verbunt and Meylan [472 ] showed that mass segregation has a significant effect on the production of CVs and X-ray binaries, an early indication that globular cluster structure needs to be treated in some detail rather than as a homogeneous static background.

As stellar modeling improved, the central role of tidal capture began to be called into question. Several authors demonstrated that close encounters between MS or RG stars lead to mergers rather than bound pairs [407 , 324 , 27 ] and thus will not produce X-ray binaries. It was also realized that, since dynamical interactions in globular cluster cores can destroy binaries [244 ], a low binary fraction in old Milky Way globular clusters did not necessarily imply a low primordial binary fraction. Thus, particularly in young clusters, interactions between stars and binaries are likely to be important.

One of the earliest systematic attempts to simulate double-degenerate binaries in globular clusters was the Rappaport et al. (1989) [403 ] study of the effect of a dense stellar environment on a wide binary containing an MSP with a low-mass white dwarf companion. The study utilized an encounter rate technique where a large number of numerical binary-single scattering experiments were conducted in a static cluster background with parameters chosen to represent $ω$ Cen and 47 Tuc. They found that binary-single interactions were capable of introducing significant eccentricity into some of the binaries, enough to affect their gravitational wave inspiral timescales. It also appeared that disruption interactions were unlikely, at least for binaries with $P ≤ 300$ days so hard double-degenerate binaries should be able to survive in globular cluster cores.

Sigurdsson and Phinney [437 , 438 ] conducted similar scattering experiments to those performed by Rappaport et al. [403 ] but with tighter binaries and a much extended mass spectrum. They discovered that exchanges were even more important than previously thought, and indeed could be the dominant type of interaction between massive stars and moderately hard binaries. A secondary effect of exchange is that when a light star is swapped for a massive one, the binary will have a larger semi-major axis for the same binding energy, which increases the cross section of the binary to undergo further interactions. They also found that three-body interactions with very hard binaries tended to produce a merger between two of the members and accretion onto the neutron star during such an event could plausibly create a MSP. Indeed Sigurdsson and Phinney argued that such interactions are sufficient to explain the entire MSP population of globular clusters. They also showed that a small number of low mass double millisecond pulsars (DMSPs) can be formed in this way.

Sigurdsson and Hernquist [435 ] considered the role of binary-binary interactions in forming DMSPs. They proposed a mechanism where two neutron star-main sequence binaries interacted and both main sequence stars were disrupted. The result was a DNS surrounded by an envelope that could be accreted to spin up both members to produce a DMSP. The study employed both point-mass scattering experiments and an SPH representation of the main sequence stars, when appropriate. The model appeared plausible but has not factored into later studies and would be a rather rare, if interesting, outcome of binary-binary interactions.

Davies and collaborators [95 , 93 ] specifically investigated the population of LMXBs in globular cluster in the same simulations described in Section 5.2.1. Davies and Benz (1995) [95 ] considered the evolution of 1000 binaries injected into various cluster cores modelled as King models with a concentration parameter $W0$ between 3 and 12 – a range that covers the concentrations of galactic globular clusters. In this study they found that the number of LMXBs produced depends on many factors including the IMF, the specific model of common envelope evolution, the assumptions made about stellar evolution and the exact concentration of the cluster. Any of these parameters can affect the numbers and types of binaries by factors of several.

Davies (1995) [93 ] concentrated on the blue straggler and compact binary population of 47 Tuc and $ω$ Cen. Particular interest was shown in comparing the efficiency of single-single and binary-single encounters in producing various kinds of binaries. Here the effect of different cluster models became clear. For the model of $ω$ Cen (massive, low concentration) blue stragglers were produced significantly more efficiently by binary-single than by single-single interactions. By contrast in 47 Tuc (lower mass, higher concentration) binary-single and single-single interactions were of comparable importance for blue straggler formation. In both systems binaries containing neutron stars seemed more likely to be formed by single-single interactions with a main sequence star than by binary-single interactions. Like the CVs, however, such interactions tended to lead to the disruption of the main sequence star and a “smothering” of the neutron star with the main sequence star’s envelope. This could lead to the formation of a Thorne–Zytkow object, a red giant star which contains a neutron star in its center [461 ]. The neutron star eventually accretes the matter, which could lead to black hole production, if enough matter is accreted, or an MSP but will not lead to an X-ray binary. A smothered system could also be the result of a binary-single interactions but was a less common outcome. For the most optimistic set of parameters Davies predicted a 4:1 ratio of smothered to non-smothered neutron star binaries in $ω$ Cen and a 5:1 ratio in 47 Tuc.

Davies predicted that $ω$ Cen would form more X-ray binaries than in 47 Tuc since the IMF in $ω$ Cen model was skewed towards higher masses and thus produced a larger number of neutron stars. By contrast more CVs were predicted in 47 Tuc. Both clusters were predicted to produce more CVs than LMXBs. This was difficult to reconcile with observations at the time but actually agrees with more recent observations. The number of predicted LMXBs was consistent with contemporary observations of 47 Tuc and $ω$ Cen (the possibility of one or two in each at the current epoch).

A new model for neutron star-white dwarf binaries was introduced by Rasio et al. (2000) [405 ] who proposed that short-period MSP-white dwarf binaries can be formed by exchanges with primordial binaries followed by a common envelope phase. In this model a $∼ 1 –3 M ⊙$ main sequence star in a binary acquires an neutron star companion through an exchange interaction and forms a new binary at a separation of 0.1 – 1.0 AU. The main sequence star undergoes Roche lobe overflow and has its envelope stripped by the neutron star. Spin-up during the accretion phase leaves an MSP in a close orbit with a white dwarf. Rasio et al. used Monte Carlo scattering simulations to predict the MSP-white dwarf binary population produced by this scenario and compared it to the population of short-period MSP-white dwarf binaries observed in 47 Tuc. The results, given in Figure 14 , show very good agreement with the observations and confirm the importance of interactions for the compact binary population in globular clusters.

Figure 14: Results of the Monte Carlo simulation of neutron star–white dwarf binary generation and evolution in 47 Tuc. Each small dot represents a binary system. The circles and error bars are the 10 binary pulsars in 47 Tuc with well measured orbits. Systems in A have evolved through mass transfer from the white dwarf to the neutron star. Systems in B have not yet evolved through gravitational radiation to begin RLOF from the white dwarf to the neutron star. Systems in C will not undergo a common envelope phase. Image reproduced by permission from Rasio et al. [405 ], copyright by IOP.

The tail end of the systems in group B of Figure 14 represents neutron star–white dwarf binaries in very short period orbits that are undergoing a slow inspiral due to gravitational radiation. These few binaries can be used to make an order of magnitude estimate on the population of such objects in the galactic globular cluster system. If we consider that there are two binaries with orbital period less than 2000 s out of $∼ 106 M ⊙$ in 47 Tuc, and assume that this rate is consistent throughout the globular cluster system as a whole, we find a total of $∼$ 60 such binaries. This is a crude estimate and may be optimistic considering the results of Ivanova et al. [247 , 246 ] but is an example of how simulations can be extended to produce gravitational wave detection estimates.

Ivanova et al. (2005) [248 ] Re-introduced the plausibility of a collision between a red giant and a neutron star forming an eccentric neutron star-white dwarf UCXB. The result was based on the hydrodynamical simulations of Lombardi et al. [295 ]. These simulations showed that while a common envelope phase is likely after such a collision, much of the envelope is ejected to infinity and an eccentric UCXB is a common outcome of the interaction. Ivanova et al. used the results from 40 SPH simulations and simple, analytic cluster modeling to estimate that it is possible for all UCXBs in globular clusters to be formed by this mechanism.

Ivanova et al. (2008) [246 ] extended the results of the CV study described in Section 5.2.1 to cover neutron star binaries as well. The dynamical situation for neutron star binaries was found to be slightly different than for CVs. About half of the neutron star binaries in the core formed as the result of an exchange interaction with only a few percent more formed by either tidal capture or collision. The rest of the binaries were primordial. Very few neutron star binaries became X-ray binaries and Ivanova et al. estimated there is only a $∼$ 3% chance that any given neutron star will ever be in an LMXB. Furthermore, most primordial binaries that became X-ray binaries did so quite early in the life of the cluster and were no longer luminous in the current epoch. Thus Ivanova et al. predicted that X-ray binaries observed in Galactic globular clusters today are likely to be products of exchange, tidal capture, or collisions. Although exchange interactions produced most of the dynamically formed neutron star binaries, X-ray binaries in their simulations were produced in equal numbers by all three channels. This was due to each channel having a different probability of producing a mass-transferring system. Finally Ivanova et al. predicted $∼$ 7.5 UCXBs in the entire Galactic globular cluster system and $∼$ 180 qLMXB systems. They also made specific predictions that certain clusters should have at least one LMXB and all of these predictions have been met.

In addition to X-ray binaries, Ivanova et al. 2008 [246 ] also considered the population of MSPs and double neutron star binaries. Overall they had difficulty reconciling the number of X-ray binaries and MSPs in their simulations with the observed ratios. If they assumed that all mass-gaining events lead to MSPs then the number produced was too high unless only $∼$ 10% of MSPs in clusters have been detected. They speculated that either not all accretion leads to spin up or that their model for mass segregation did not work properly for these systems and more should have sunk to the core and be disrupted. This posed a problem, however, because all of the channels that lead to MSP production were necessary in order to produce the correct number of LMXBs. As a solution they proposed that physical collisions or accretion onto post-accretion-included collapse systems produces MSPs with strong magnetic fields. Such MSPs are short-lived and would not appear as MSPs at current GC ages. They also proposed that mass transfer in NS-WD LMXBs does not lead to the formation of binary MSPs. By using these modifications and adjusting their assumptions about spin-up during various kinds of accretion they were able to bring the number of MSPs down to a level that is consistent with observations of 47 Tuc and Terzan 5. They also noted that the fraction of single MSPs to MSPs in binaries tends to increase with increasing cluster density. This implies that destructive interactions dominate for MSP binaries at high density.

They found that the production of double neutron star binaries was very inefficient. There were only three primordial DNSs in the $7 ∼ 2 × 10 M ⊙$ of stars that make up their entire simulation set while only 14 double neutron stars were formed by dynamical interactions. Given their chosen IMF and cluster ages, this implies that it takes $∼$ 500 single NSs to produce one merging double neutron star binary in 11 Gyrs. These results were very density dependent and the vast majority of the DNSs formed in the highest density simulations. They estimated that merging DNSs in globular clusters could provide at most 10 – 30% of short gamma ray bursts but noted that this agrees fairly well with the predictions of Grindlay et al. (2006) [181 ].

This result also agrees with results of Downing et al. (2010) [111 ] who used a similar-sized Monte Carlo dataset, primarily to study black hole-black hole binaries, and found no DNSs in any of their simulations. That DNSs are not produced in similar numbers to X-ray binaries in globular clusters may be a combination of the lower number of potential progenitors, kicks for both progenitors, and the likelihood of merging or producing a WD during interactions.

Ivanova et al. also reported a significant number of LISA sources in their simulations. They predicted that there should be an average of 10 LISA sources per globular cluster at any given instant and that on average 180 LISA sources are produced per globular cluster per Gyr. Scaling this to the Milky Way, they found that there may be up to 1500 LISA binaries, both DWD and white dwarf-neutron star, in the Galactic globular cluster system at the current time. They cautioned, however, that they may have overestimated the white dwarf-neutron star binary population since it was not fully consistent with UCXB observations.

Banerjee and Ghosh [29 , 30 ] employed a novel method to investigate the compact binary population of globular clusters. They used the collisional Boltzmann equation as given by Spitzer [443 ] to calculate the evolution of the number distribution of binaries in a uniform, non-evolving background of stars. They did not employ the Fokker–Planck approximation but rather considered mean encounter rates with a randomly fluctuating component. This is called a Wiener process and is a mathematical description of Brownian motion. The Wiener process means that the encounter rates, normally ordinary differential equations in the Boltzmann equation, become stochastic differential equations and must be treated with a mathematical tool known as Itō calculus. Both, Wiener processes and Itō calculus are described in Banerjee and Ghosh (2008) [30 ].

Banerjee and Ghosh specifically considered a bimodal population of $0.6 M ⊙$ main sequence stars interacting with $1.4 M ⊙$ neutron stars. They consider binary formation by tidal capture and exchange, destruction by 3-body disruption of the binary or by exchange of the compact component for a main sequence star, and hardening by magnetic breaking and gravitational radiation. They were able to re-produce the build-up of short-period binaries found by other authors but, more interestingly they could compare the mean rates predicted by the Boltzmann equation to the random variations introduced by the stochasticity of the Wiener process. Thus they were able to determine if it is possible to predict the number of X-ray sources expected in a given globular cluster or if stochastic effects will dominate. Figure 15 shows both the mean number of X-ray binaries and the stochastic fluctuation expected for various cluster parameters as well as observed numbers from several Galactic globular clusters. Although the number of X-ray binaries in the models and observations follow the same trend, the stochastic fluctuations are large and indicate it is probably impossible to predict exact numbers of compact binaries in any given globular cluster. Thus, although the method of Banerjee and Gosh is in some ways less detailed than that of Ivanova et al., it is still valuable since it constrains the effects of stochasticity in a quantitative way.

Figure 15: Contours of constant number of X-ray binaries as a function of the “Verbunt parameters”. Here $Γ = ρ2r3c∕vc$ , $γ = ρ∕vc$ , $ρ$ is the average stellar density, $rc$ is the core radius and $vc$ the core velocity dispersion. $Γ$ describes the two-body stellar encounter rate and $γ$ describes the rate of encounters between a single binary with the stellar background. The squares represent the location of various Galactic globular clusters in the $γ − Γ$ plane and the number in parentheses is the number of X-ray binaries observed in the cluster. Plot reproduced by permission from Banerjee and Ghosh [30 ], copyright by IOP.

Recently Bagchi and Ray [21 , 22 ] used scattering experiments to re-investigate the eccentricities of MSP binaries in light of improved observations. They noted that observationally there are three, statistically distinct, eccentricity groupings in MSP binaries, those with $e ∼ 0$ , those with $− 6 2 × 10 < e ≤ 0.01$ and those with $0.01 ≤ e$ . They performed scattering experiments, utilizing the STARLAB package, to attempt to determine the origin of the different groupings. They found that most of the higher eccentricity binaries could be explained by exchanges or mergers during binary-single or binary-binary interactions. The intermediate eccentricities could be explained by fly-by encounters while the lowest eccentricities were binaries that had been circularized without experiencing fly-bys. They also found one binary that they could only explain by invoking resonant interactions within a hierarchical triple [276 ].

Of a more speculative nature was the work of Trenti et al. (2008) [464 ] considering the number triple systems that might contain at least one pulsar. They studied triple formation using direct N-body simulations with between 512 and 19 661 bodies, both in equal mass systems and systems with a mass function, for various initial binary fractions. They found that triple systems were about two orders of magnitude less common than binary systems within globular cluster cores but, given the number of MSP binaries currently known in globular clusters, this implies that a triple containing an MSP should be detected soon. Indeed, the LMXB 4U 1820–30 is an excellent candidate for such a system. 4U 1820–30 is a WD-NS binary in the globular cluster NGC 6624 with a period of $∼$ 685 seconds [447 ] and a factor-of-two luminosity variation with a period of $∼$ 171 days [394 , 73 ]. It has been shown that this system is probably a hierarchical triple [175 ] where the third body is in a Kosai resonance with the inner binary [491 , 395 ]. The most recent work of Prodan and Murray [395 ] gives a detailed investigation of the nature of the Kosai resonance in 4U 1820-30 and a plausible dynamical history. However these authors state that a better understanding of both the long-term evolution of compact X-ray binaries and the potential of NGC 6624 are necessary to fully understand this system. 4U 1820–30 promises to be a very important system for understanding some of the more complicated dynamics that can affect compact binaries in globular clusters.

There are several overall conclusions that can be drawn about the population of neutron star binaries in globular clusters. Dynamical interactions can certainly enhance the number of X-ray binaries but the current consensus is that the enhancement will by factors of a few, not factors of hundreds. Recent simulations show that detailed descriptions of dynamics, binary stellar evolution, and stellar collisions are all necessary in order for the results to be reliable. The same interactions that produce X-ray binaries are also effective at producing MSPs, although many of the resulting binaries will be subsequently disrupted and the MSPs will continue their lives as single objects. It also seems that neutron star-white dwarf binaries may be fairly common in globular clusters. By contrast detailed simulations show that DNS binaries are not common in globular clusters since few exist primordially and they are not produced efficiently by stellar dynamics. The individual results of these studies are highly dependent on the detailed treatment of various phases of binary stellar evolution as well as the kick distribution. It is worth noting that many authors, using a wide variety of approximations, have successfully predicted the number of X-ray binaries in globular clusters. This indicates that it is actually reasonably easy for simulations to re-produce the X-ray binary population and could be taken as a vindication of the dynamical models. However it also indicates that X-ray binaries are a rather generic outcome of globular cluster dynamics. Therefore successfully re-producing the observed number of X-ray binaries in globular clusters may not be a very helpful diagnostic when trying to constrain our assumptions about globular cluster dynamics.

5.2.3 Binaries with black hole primaries

Black holes are non-luminous and can only be detected in the electromagnetic spectrum if they are accreting matter, normally from a binary companion. Mergers of black holes and neutron stars may be detectable in both gravitational waves and as short gamma-ray bursts [182 ]. Black hole-black hole binaries (DBHs) can only be detected in gravitational radiation.

Black holes are the most massive compact stellar remnants, up to $35 –80 M ⊙$ [141 , 43 ], and can, on average, be several times more massive than the average mass of stars in even a fairly young globular cluster. Thus black holes are strongly affected by mass segregation, may form Spitzer-unstable subsystems in the dense core of the cluster and will be preferentially exchanged into binaries during interactions. As such, the black hole binary population in globular clusters will be particularly strongly affected by dynamics.

DBHs are of great interest for gravitational wave detection since they are quite massive compared to neutron star mergers and may be detectable by the next generation of ground based detectors at Gigaparsec distances. Strong dynamical interactions should produce many hard DBHs in globular clusters. It is not clear, however, that globular clusters themselves will host a large number of DBH mergers because the binaries could be ejected from the cluster cores before they reach the gravitational wave regime. Kulkarni et al. (1993) [281 ] and Sigurdsson and Hernquist (1993) [436 ] presented the first detailed studies of black hole-black hole binaries in globular clusters. They used purely analytic estimates to show that black holes would indeed mass segregate and form a Spitzer unstable subsystem in the cluster core. In this subsystem black holes interact with each other to form binary and triple systems. Interactions between these systems and the remaining black holes remove both from the cluster very efficiently and Sigurdsson and Hernquist in particular predicted that there should be only zero to four black holes left in the cores of old galactic globular clusters, with perhaps a few more in the cluster halos. Remarkably, this basic picture of mass segregation, dynamical interaction, binary formation and ejection has not changed much since these pioneering publications.

Studies of black hole binaries remained dormant until the direct N-body simulations of Portegies Zwart and McMillan (2000) [388 ]. Using small simulations of between 2048 and 4096 particles with a bimodal mass distribution, 0.5 – 1% of bodies ten times more massive than the rest, they essentially confirmed the results of Kulkarni et al. and Sigurdsson and Hernquist. They found that the black holes mass segregate rapidly, form binaries through three-body interactions, and are then ejected from the clusters. Overall they found only $∼$ 8% of black holes were retained by their parent clusters while $∼$ 61% were ejected as singles and $∼$ 31% were ejected as binaries. However, many of the ejected, dynamically-formed DBHs were tight enough to merge in the galactic field outside the clusters within a Hubble time. Indeed, extrapolating their results to globular clusters they estimated up to one advanced LIGO detection of a dynamically-formed DBH binary merger per day. Thus Portegies Zwart and McMillan concluded that globular cluster would be an important source of DBH binaries. It must be noted, however, that these simulations did not include stellar evolution or natal kicks for black holes.

Further studies of the behaviour of black holes in globular clusters were prompted by the suggestion of Miller and Hamilton [334 ] that IMBHs of up to several hundred solar masses could be formed by repeated mergers of black holes in dense stellar environments (see Section 5.3 for further information). Gültekin et al. (2004) [186 ] performed repeated three-body scattering experiments between single black holes and DBH binaries of various masses in order to investigate such growth. They used the static cluster background density to calculate an encounter timescale between interactions and, if the quadrapole gravitational wave radiation timescale was less than the encounter timescale, a merger occurred. They confirmed that the DBHs are ejected before an IMBH is formed but that many of the ejected binaries pass thorough the LISA band and then merge in the VIRGO/LIGO band within a Hubble time. They also found that many of the binaries retain significant eccentricity in the LISA band.

O’Leary et al. (2006) [360 ] considered the same problem using a slightly different approximation. They assumed all black holes were completely mass-segregated and formed a Spitzer-unstable subsystem interacting only with itself. The they followed the evolution of this subsystem, embedded in a globular cluster potential, for a range of initial conditions using a Monte Carlo selection of encounters and numerical integration of binary-single and binary-binary interactions. In some of their simulations they included a prescription for gravitational wave recoil after a DBH merger (see [345 ] for further discussion of this effect). O’Leary et al. also found little evidence of IMBH growth but their simulations did form a large number of DBH binaries, many of which merged within a Hubble time. They quantified the detection rate for advanced LIGO, using simplified cosmological assumptions and assuming optimal orientation, and predicted 10s of black hole merger detections per year (the exact number depending on the initial conditions of the cluster model). They found that between 50% and 70% of these mergers will take place between DBHs that were ejected from the clusters. Again, the clusters produce many DBHs but do not retain them.

Both Gültekin et al. and O’Leary et al. included only self-interacting black hole systems in their simulations. By contrast, Sadowski et al. 2008 [415 ] modeled entire globular clusters using the same Monte Carlo method as Ivanova et al. [247 ] and assumed that the black holes remained in thermal equilibrium with the rest of the cluster. This is the opposite dynamical assumption to that made by O’Leary et al. since it assumes that the black holes are minimally mass segregated. It maximizes the number of interactions between black holes and other cluster stars while minimizing the number of interactions between black holes and other black holes. The main result of this is that there are fewer DBHs formed but those that do form are far less likely to experience a further encounter with another black hole. Interactions with lower-mass cluster members are less energetic and thus less likely to disrupt or eject DBHs. In practice the lower destruction and ejection rate wins out and Sadowski et al. predicted a very large number of DBH binaries in globular clusters. Indeed they estimated an advanced LIGO detection rate of 25 – 3000 per year, 90% of which will occur within the clusters. When compared with the results of O’Leary et al. this proves that the treatment of global star cluster dynamics has a strong effect on the DBH population.

Two recent sets of studies have tried to address the problem of more accurate global dynamics. The Monte Carlo simulations of Downing et al. 2010 and 2011 [111 , 112 ] and the direct N-body simulations of Banerjee et al. 2010 [28 ]. Downing et al. used the Monte Carlo code of Giersz [153 , 155 ] to perform 160 fully self-consistent simulations of globular clusters with 500,000 particles, a realistic mass function, stellar evolution, 10 – 50% primordial binaries, and various combinations of initial concentration and metallicity. The results of the simulations tend to support the approximation of O’Leary et al., namely the black holes mass segregate very quickly and form a dense system in the core of the cluster. Downing et al. confirmed that many DBH binaries were formed, particularly by three-body interactions, but many of these binaries were either destroyed in subsequent interactions or were ejected from the clusters. In dense clusters DBH production was strongly time-dependent. Many DBHs were formed early but this number dropped substantially after a few Gyrs as the clusters were depleted of black holes. By contrast the clusters with a low initial concentration did not suffer as strong a depletion in their black hole population and produced DBHs slowly but constantly throughout their lives. The number of black hole binaries produced was highly stochastic with large simulation-to-simulation variations, even between models with the same initial conditions, as shown in Table 6. Depending on the initial conditions Downing et al. predicted between 1 – 100 DBH merger detections per year for advanced LIGO with an average of $∼$ 15 per year. This is consistent with, although slightly higher than, O’Leary et al. The detection rate also showed simulation-to-simulation variation and was similar for some simulations with different initial conditions. This implies that the DBH merger detection rate alone will probably not carry much astrophysical information. Downing et al. also predicted the possibility of eccentric DBH binaries in the low-frequency gravitational wave radiation band.

Table 6: The cumulative number of BH–BH binaries after 3, 9, and 25 $trh$ , and also after 1 $T H$ . The first column gives the initial conditions. 10 or 50 in the first place indicates 10% or 50% primordial binaries. sol or low indicates a metallicity of z = 0.02 or z = 0.001. The number in the third place indicates $rt∕rh$ = 21, 37, 75 or 180. Each column is averaged over ten independent realizations and includes the $1σ$ uncertainty. A dash in a column indicates that the cluster did not reach that number of half-mass relaxation times within one Hubble time. Reproduced by permission from Downing et al. 2010 [111

]

Simulation	$t = 3 trh$	$t = 9trh$	$t = 25 trh$	$t = 14 Gyr$
10sol21	1 ± 1	–	–	1 ± 1
10sol37	1 ± 1	14 ± 11	–	14 ± 11
10sol75	0 ± 0	8 ± 6	49 ± 19	52 ± 19
10sol180	0 ± 0	12 ± 6	54 ± 21	123 ± 27
50sol21	1 ± 1	–	–	3 ± 2
50sol37	3 ± 2	36 ± 10	–	50 ± 11
50sol75	1 ± 1	26 ± 8	115 ± 24	147 ± 28
50sol180	0 ± 1	11 ± 5	111 ± 23	354 ± 33
10low21	22 ± 10	–	–	27 ± 10
10low37	23 ± 4	44 ± 6	–	44 ± 6
10low75	18 ± 10	32 ± 13	54 ± 6	54 ± 16
10low180	26 ± 8	51 ± 8	79 ± 20	112 ± 24
50low21	104 ± 16	–	–	127 ± 16
50low37	93 ± 22	175 ± 26	–	184 ± 29
50low75	64 ± 10	155 ± 22	173 ± 22	202 ± 38
50low180	103 ± 19	205 ± 35	294 ± 50	453 ± 109

Banerjee et al. used NBODY6 to perform direct N-body simulations of young and intermediate age massive clusters with between 5000 and 100 000 stars, including both stellar evolution and a realistic IMF. The also performed a series of 3000 to 4000 body simulations surrounded by a reflective boundary to model the core of a globular cluster. They found that the black holes mass-segregated quickly, on timescales of 50 – 100 Myr, almost independently of the mass of the cluster. Once mass-segregated the black holes interacted strongly, making many DBHs during three-body encounters. Like most previous authors, Banerjee et al. found that the black holes and black hole binaries are efficiently disrupted or ejected and that their simulated clusters were depleted of black holes on a 4.5 Gyr timescale. Thus they predicted that current galactic globular clusters will not contain black holes. They also estimated that intermediate age massive clusters will produce $∼$ 31 black hole binary merger detections per year for advanced LIGO. This is of the same order as found by Downing et al. and slightly higher than the rate found by O’Leary et al. while it is about an order of magnitude lower than that found in the previous N-body simulations of Portegies Zwart and McMillan. They cautioned, however, that this rate is not for the old globular clusters that were the focus of the previous studies and thus care must be taken in the comparison.

There have also been efforts made to include Post-Newtonian (PN) dynamics into direct N-body simulations, both for stellar-mass black holes [284 ] and for black holes around a central massive object [5 ]. The algorithm has been tested successfully and used to discover significant new dynamical behaviour around central massive objects, particularly in the context of extreme mass ratio inspirals (e.g., the “butterfly effect” [13 ]). However, in simulations of globular clusters, where typical velocities are low and the likelihood of relativistic encounters is small (see, e.g., [14 ]), it is not clear that the inclusion of PN terms has a significant effect on the DBH population of a cluster [6 ]. For stellar-mass binaries, and particularly for estimating relativistic merger rates, a Newtonian treatment of the cluster dynamics with relativistic effects added in post-processing may be sufficient.

In 2010, Ivanova et al. [245 ] considered the possibility of black hole-white dwarf X-ray binaries as UCXBs in globular clusters. This study employed scattering experiments and a simplified treatment of cluster dynamics. Ivanova et al. determined that stellar collisions, exchange interactions and three-body binary formation can all produce black hole-white dwarf X-ray binaries, however in most cases the binary will be too wide to become an X-ray source within a Hubble time. Therefore most black hole-white dwarf X-ray binaries must have experienced dynamical hardening, either through fly-by encounters or as part of a hierarchical triple where the inner black hole-white dwarf binary has its separation reduced by the Kozai mechanism [276 ]. They found that each formation mechanism, coupled with dynamical hardening, could explain the observed number of ULXBs in globular clusters but that all required $≥$ 1% of black holes remain in the cluster upon formation. This places limits on the possible natal kick distribution of black holes. Interestingly, recent observations of the black hole X-ray binary XMMU 122939.7+075333 in the globular cluster RZ 2109 (itself in the Virgo cluster galaxy NGC 4472) shows variability that may be the result of the Kozai mechanism [309 ]. This strengthens the case for black hole X-ray binary hardening in hierarchecal triples.

Unlike some of the other species mentioned, the basic behaviour of black hole-black hole binaries in globular clusters is well-understood. The black holes mass segregate, interact, form binaries and are swiftly ejected. The focus of future investigations will be on event rates and details of mass distributions. It seems certain, however, that globular clusters will be a major source of DBH binary mergers. When compared with the results for field DBH binaries of Belczynski et al. 2007 [48 ] all of the studies cited indicate that dynamically formed binaries will dominate the DBH binary merger detection rate for advanced LIGO by factors of several. We note that the most recent estimates of field DBH rates for advanced LIGO include a range of rates from pessimistic to highly optimistic that spans over three orders of magnitude [7 ], and so substantial uncertainty remains. Some more recent results [44 ] indicate that the number of field black hole binaries may be larger than previously estimated but dynamically formed binaries are still expected to be a major source of detections. Including younger clusters before they have been depleted in black holes may be important for estimating event rates and will probably only cause them to increase. Should the next generation of ground-based detectors perform according to expectations, detection of a DBH binary inspiral is almost certain.

5.3 Intermediate-mass black holes

Intermediate-mass black holes (IMBHs) are a hypothetical class of black holes with masses of $102– 105M ⊙$ . Ultimately IMBH is a catch-all term for BHs in the mass range between the stellar mass BHs ( $2 < 10 M ⊙$ ) that are known to exist from stellar evolution and and the super-massive BHs ( $> 105M ⊙$ ) that are known to exist at the centers of galaxies. It is as yet unclear whether such objects exist and, if they do, whether they exist as binaries in globular clusters. If they do exist, a binary of two IMBHs or an IMBH and a stellar mass BH in a globular cluster could be a very interesting source of gravitational wave radiation.

There are several methods that can be used to detect IMBHs. The most commonly used are X-ray or radio accretion luminosity, burst emission from tidal disruption of stars and the structure and kinematics of the inner regions of globular clusters.

If IMBHs exist as close binaries with any companion other than a neutron star or black hole they will accrete matter and the accretion will produce X-ray emission and can be detected, as is the case for stellar-mass black holes in binaries. The mass of the IMBH can be estimated from the total X-ray luminosity. There are several extra-galactic ultra-luminous X-ray sources (ULXs) where the luminosities are consistent with an IMBH and that are provisionally associated with globular clusters [294 , 122 , 396 , 151 , 272 , 274 , 380 ]. While suggestive these observations are not conclusive since X-ray observations often lack sufficient resolution to definitively associate them with the proposed host clusters. Some are also consistent with stellar-mass X-ray binaries [380 , 274 ] or have been resolved into multiple point sources which are then more consistent with stellar-mass X-ray binaries [359 ]. Another possibility for detection would be radio emission from Bondi accretion of ambient material from the ISM [305 ]. Several studies used this method to place upper limits on the possible mass of IMBHs in Galactic globular clusters [98 , 465 , 32 , 300 ] but there are no definitive detections. Indeed many studies are consistent with there being no IMBHs in the Galactic globular clusters (e.g., [311 , 312 ].

Another possible way of detecting an IMBH is through tidal disruption of a passing star. Simulations have shown that if a globular cluster contains an IMBH then the IMBH is quite likely tidally capture stars that pass close to it [36 ]. If these stars are tidally disrupted they will produce a short soft X-ray flare followed by a several-hundred year X-ray afterglow from the accretion disk [398 ] that may be detectable. There is some evidence for the disruption of a white dwarf in two extra-galactic globular clusters [243 , 313 ]. If an IMBH is in a binary with another cluster star it may also be possible for it to exchange it for a millisecond pulsar and form an IMBH-millisecond pulsar binary. Although the lifetime of such a binary would be short due to gravitational wave radiation, it might be detectable with future radio observations [101 ]. The merger of such a system might also be an interesting gravitational wave source.

Perhaps the strongest evidence for IMBHs in globular clusters comes from an analysis of the structure and kinematics of their core regions. A central massive object, such as an IMBH, affects both the density profile and velocities of the stars surrounding it. If a globular cluster contains an IMBH it will sink to the center due to dynamical friction and will form binaries with the compact objects in the center, particularly the black holes. This will be a source of dynamical heating for the cluster [38 , 463 , 275 ] and will cause the core of the cluster to expand. This means that clusters with an IMBH will have relatively large cores with constant density and without a significant peak in surface brightness at the center [39 , 356 ]. This may lead to a rather large ratio of $rc∕rh$ in clusters with IMBHs and this could be used to narrow the search for potential IMBH hosts (although various other processes in globular clusters can mimic this effect [228 ]). The activity of the IMBH in the core of a globular cluster will also increase both the velocity dispersion and the mass-to-light ratio in the cluster core. Structural and kinematic data suggests possible IMBHs in the Galactic globular clusters $ω$ Cen [357 , 358 , 467 , 250 ], M15 [148 , 149 , 466 ], NGC 6388 [285 , 303 ], 47 Tuc [322 ] and the extra-galactic globular cluster G1 [40 , 146 ]. In all cases the data are merely consistent with an IMBH being the best-fitting model given the dynamical assumptions and give only upper limits on possible IMBH masses. Of all the candidates G1 has perhaps the strongest dynamical evidence and is supported by X-ray observations [274 ]. $ω$ Cen is also a reasonably good candidate, although the current upper mass limit on the IMBH is disputed by a factor of 3 – 4 [358 , 467 ]. By contrast the kinematic signal in 47 Tuc is weaker and the signal in M15 is disputed [314 ].

All of these promising but uncertain detections beg the question as to whether we should expect IMBHs in globular clusters to begin with. There are several mechanisms proposed to form IMBHs in globular clusters, the most straightforward of which is the collapse of zero-metallicity (Pop III) stars. Pop III stars may have characteristic masses of $≥ 100 M ⊙$ [9 ] and can collapse in the same manner as today’s stars but leave more massive remnants [141 ]. Such a collapse could produce BHs of up to a few hundred solar masses. Another possibility is the direct collapse of a $∼ 105M ⊙$ cloud of gas to a BH through post-Newtonian instabilities [41 ]. Both of these mechanisms are proposed for forming the seed black holes for SMBHs in young dark matter halos in the early universe. However globular clusters are not thought to be a result of cosmological structure formation as they form later and not necessarily in the same location as Pop III stars. Therefore there is no particular reason to expect that either of these mechanisms will lead to IMBHs in globular clusters.

Figure 16: Total number of black holes ejected in building up to a $1000 M ⊙$ IMBH as a function of seed mass. The four numbered curves are made assuming different cluster escape velocities. The dotted line is 10³ black holes, the assumed population for a young cluster. As can be seen, a very large seed mass is necessary before IMBH formation by merging black holes becomes a possibility. Image reproduced by permission from Gültekin et al. [186

], copyright by IOP.

There are, however, ways of forming IMBHs that occur only in globular clusters. These include repeated collisions and mergers of stellar mass black holes and the formation of very massive stars (VMSs) through the runaway mergers of young massive stars. The collision and merger of black holes model, first proposed by Miller and Hamilton [334 ] is simple, DBH binaries form by the processes described in Section 5.2.3 and then merge due to gravitational wave inspiral. The resulting merger product becomes more massive, is more likely to be exchanged into BH binaries and by repeated mergers grows into and IMBH. This model originally seemed promising however, as discussed in Section 5.2.3, the three-body numerical scattering experiments of Gültekin et al. (2004) [186 ] showed that it was not possible to create a BH with a mass greater than $∼ 240M ⊙$ for any reasonable seed mass. Either the proto-IMBH–BH binary was ejected from the cluster or the core depleted of BHs due to scattering interactions before this point. Their final mass limits for different seed masses are illustrated in Figure 16 . The results of O’Leary et al. (2006) [360 ] tend to confirm this result, finding similar (although slightly higher) limiting masses in their simulations where merger products do not receive gravitational wave recoils and almost no growth at all when recoils were included. It should be noted that O’Leary et al. assume total mass segregation and use fairly massive black holes. This means that their interactions will be frequent and strong and thus the ejection rate maximal. However, many simulations using more realistic models of cluster dynamics have now confirmed that DBH binaries efficiently eject each other from globular clusters [388 , 111 , 112 , 28 ] and it is unlikely many would remain long enough to merge and form a seed IMBH before they are ejected during binary-binary interactions.

A more promising path for the formation of IMBHs is the merger of young massive stars in globular clusters while the clusters themselves are still young. In this model, first suggested in the context of IMBHs by Portegies Zwart and McMillan (2002) [389 ], young star clusters with short relaxation timescales ( $≲$ 25 Myr) can experience mass-segregation enhanced core collapse on timescales shorter than the lifetimes of massive stars ( $≲$ 5 Myr for stars of $∼ 50 M ⊙$ ). Such stars have large radii and the combination of mass-segregation and core collapse can provide a region of sufficiently high stellar density for collisions between them to become likely. If sufficiently high central densities are reached the product of one collision can rapidly experience another collision in a runaway process, resulting in the production of a VMS. Several studies suggest that a star of up to several times $103M ⊙$ can be formed by runaway collisions in young, dense star clusters [389 , 385 , 189 , 139 , 138 , 173 ] for a variety of plausible initial conditions. For massive, dense clusters an increase in central potential by gas accretion onto the core could enhance this effect [75 , 341 , 37 ]. Collisions during binary interactions may also play a role and there is a suggestion that they could lead to the formation of a pair of VMSs – a possible path to an IMBH binary [188 ].

Figure 17: The effect of stellar wind mass-loss on the evolution of a runaway merger product during core hydrogen and core helium burning. The blue curve corresponds to the mass-loss used by Portegies Zwart et al. (2004) [385 ]. The black curve corresponds to the mass-loss used by Belkus et al. (2007) [49 ] with Z = 0.02 while the red line is for the same models but Z = 0.001 Image reproduced by permission from Vanbeveren et al. [468

], copyright by Springer.

The runaway collision scenario is more promising that the BH merging scenario for a variety of reasons. Young massive stars are much larger than BHs, making direct collisions between them more common. The fact that they are not compact objects means that tides can be raised during interactions which leads to easier capture and less chance of scattering. Finally, the relativistic recoils expected after the merger of two BHs will not occur in the merger of two massive stars and the merger product is more likely to remain in the cluster. However, there are some problems with the runaway merger scenario. The first is the short relaxation times and high stellar densities necessary to achieve a runaway merger. There seem to be no star-forming regions or young open cluster in the Milky Way with a sufficiently high concentration to produce a runaway merger [19 , 37 , 341 ]. The formation conditions for globular clusters are not know and it is certainly possible that they form at much higher densities. However, if they form from the hierarchical merging of many sub-clusters they may not reach a sufficient central density for runaway collapse before the massive stars have evolved away from the main sequence. It is also possible that there is significant mass-loss during the collisions and mergers (up to 25% [143 ]). Such mass-loss could slow the growth of the merger product. Furthermore this mass may be lost asymmetrically giving the merger product a velocity kick, removing it from the core of the cluster, further hindering its growth. Another problem is the evolution of the merger product itself. Such a massive object may be highly unstable, will have a short lifetime and may suffer significant mass-loss due to stellar winds. The relationship between the mass-loss and the collision timescale is key because if mass is lost to stellar evolution faster than it can be replaced by collisions then the merger product cannot grow. There are both studies that predict that the merger timescale is shorter than the mass-loss timescale [138 , 456 ] and studies that predict that the mass-loss timescale is shorter than the collision timescale [165 , 468 ], at least for clusters of solar metallicity. Some of the different possible outcomes are illustrated in Figure 17 . A high merger rate can also cause problems if further collisions occur before the merger product has regained thermal equilibrium. The outcome of such a merger and its effect on the evolution of the growing massive star is difficult to predict but might lead to the disruption of the VMS. Even the equilibrium evolution of such a star is poorly understood and it is not clear if the remnant of such evolution will actually be an IMBH [138 ].

Finally, it is not clear that, if formed, an IMBH binary will remain in a globular cluster for very long. Binary-binary scattering with DBH binaries could remove a low-mass IMBH–BH binary from the cluster in a similar way to a standard DBH binary while the gravitational wave recoil from the merger of a $1000 M ⊙$ IMBH with a $50 M ⊙$ BH has a greater than 30% probability of removing the merer product from the cluster [223 ]. Overall IMBHs – and by extension IMBHs in binaries – remain an intriguing but speculative component of globular clusters.

	Abstract
1	Introduction
2	Globular Clusters
	2.1	Stellar populations in globular clusters
	2.2	The structure of globular clusters
	2.3	The dynamical evolution of globular clusters
3	Observations
	3.1	Cataclysmic variables
	3.2	Low-mass X-ray binaries
	3.3	Millisecond pulsars
	3.4	Black holes
	3.5	Extragalactic globular clusters
4	Relativistic Binaries
	4.1	Binary evolution
	4.2	Mass transfer
	4.3	Globular cluster processes
5	Dynamical Evolution
	5.1	Star cluster simulation methods
	5.2	Results from models of globular clusters
	5.3	Intermediate-mass black holes
6	Prospects of Gravitational Radiation
7	Summary
	Acknowledgements
	References
	Updates
	Figures
	Tables