1 Introduction

For the purposes of this review, relativistic binaries are systems containing two stellar mass degenerate or collapsed objects that are in close orbits. In the Galactic field, where these systems evolve in relative isolation, their final properties are set solely by their initial conditions and are the result of mass transfer, common envelope evolution, or other interactions that may interrupt the course of single star evolution due to the presence of a nearby neighbor. When considered as a fraction of the total stellar mass, the number of relativistic binaries in Galactic globular clusters is overrepresented compared to the Galactic field. Thus, the dynamical interactions found in the environment of dense stellar clusters provide additional channels for the formation of these systems. Relativistic binaries reveal themselves observationally as UV or X-ray sources and are potential sources of gravitational radiation.

This review will concentrate on the Galactic globular cluster system for the bulk of the text. We shall touch on extra-galactic globular cluster systems briefly in the sections on observations (Section 3) and gravitational radiation prospects (Section 6).

We begin in Section 2 by looking at the physical structure and general history of the galactic globular cluster system that leads to the concentration of evolved stars, stellar remnants, and binary systems in the cores of these clusters. Current observations of globular clusters that have revealed numerous populations of relativistic binaries and their tracers are presented in Section 3. We also consider the prospects for future observations in this rapidly changing area. Many relativistic binaries are the product of stellar evolution in close binaries. In Section 4, we will look at how mass transfer between one star and a nearby companion can dramatically alter the evolution of both stars. The enhanced production of relativistic binaries in globular clusters results from dynamical processes that drive binaries toward tighter orbits and that preferentially exchange massive, degenerate objects into binary systems. Numerical simulations of globular cluster evolution, which can be used to predict the rate at which relativistic binaries are formed, are discussed in Section 5. These models can be compared with the observable members of the population of relativistic binaries in order to try and constrain the entire population. Finally, we conclude with a brief discussion of the prospects for observing these systems in gravitational radiation in Section 6.

Readers interested in further studies of the structure and evolution of globular clusters are invited to look at Binney and Tremaine [57], Spitzer [443], and Volumes I and II of Padmanabhan’s Theoretical Astrophysics [362, 363] for a good introduction to the physical processes involved. Review articles of Meylan and Heggie [326] and Meylan [325] also provide a comprehensive look at the internal dynamics of globular clusters. Although our focus is mainly on the Galactic globular cluster system, the physics of globular cluster systems associated with other galaxies is well covered in the review article by Harris [192] as well as his lecture notes from the Saas-Fee course on star clusters [65]. Carney has a thorough introduction to evolution of stars in globular clusters [66]. An observational perspective on the role of binaries in globular clusters is presented in an excellent review by Bailyn [25], while a good introduction to the details of observing binary systems in general can be found in An Introduction to Close Binary Stars [215]. Although slightly out of date, the review of binaries in globular clusters by Hut et al. [241] is an excellent introduction to the interaction between globular cluster dynamics and binary evolution, as is a short article on globular cluster binaries by McMillan, Pryor, and Phinney [323]. Rappaport et al. [402] and Rasio et al. [404] have written reviews of numerical simulations of binary populations in globular clusters. An excellent introduction to the astrophysics and numerical techniques relevant to globular cluster dynamics can be found in the book by Heggie and Hut [199]. Finally, a shorter and more observationally focused review of compact objects in globular clusters can be found in Maccarone and Knigge [306].

2 Globular Clusters

Globular clusters are gravitationally bound associations of 104 − 106 stars, distinct both from their smaller cousins, open clusters, and the larger, dark matter dominated dwarf galaxies that populate the low-mass end of the cosmological web of structure. Globular clusters are normally associated with a host galaxy and most galaxies, including the Milky Way, are surrounded and penetrated by a globular cluster system. A good estimate of the number of globular clusters in the Milky Way is the frequently updated catalogue by Harris [193], which has 157 entries as of 2010. Although fairly complete, a few new clusters have been discovered in recent years at low Galactic latitudes [234, 271] and there may be more hidden behind the galactic disc and bulge. The distribution of known globular clusters in the Galaxy is given in Figure 1. Other galaxies contain many more globular clusters and the giant elliptical M87 alone may have over 10000 [194]. The richness of the globular cluster system of a galaxy can be classified by the number of globular clusters associated with the galaxy normalized to its luminosity. One widely used measure of this is the specific frequency, \({S_N} = {N_{{\rm{GC}}}} \times {10^{0.4}}({M_{V + 15}})\) where NGC is the number of globular clusters and MV is the V-band magnitude of the galaxy [195]. SN can vary significantly between different galaxy types. For instance Sn ∼ 1.3 for the local spiral galaxy M31 while Sn ∼ 14.4 for M87. On the whole Sn seems to be higher in massive elliptical galaxies than in spiral galaxies. For more information on extragalactic globular cluster systems see the review by Brodie & Strader [61].

Figure 1
figure 1

Globular cluster distribution about the galaxy. Positions are from Harris [193] and are plotted as black circles on top of the COBE FIRAS 2.2 micron map of the Galaxy using a Mollweide projection. Image reproduced from Brian Chaboyer’s website [69].

Milky way globular clusters are old, having typical ages of 13 Gyr and an age spread of less than 5 Gyr [67]. This is on the order of the age of the Galaxy itself, thus Galactic globular clusters are thought to be left over from its formation. By contrast other galaxies such as the small and large Magellanic clouds (SMC and LMC) have intermediate age globular clusters (< 3 Gyr old, e.g., [346, 337]) and in some galaxy mergers, such as the Antennae, massive star-forming regions that may become globular clusters are observed [117]. Taken together, this implies that globular clusters of all ages are relatively common objects in the universe.

2.1 Stellar populations in globular clusters

Most of the detailed information on stellar populations in globular clusters comes from those in the Milky Way since only they are close enough for stars to be individually resolved. The stars in individual Galactic globular clusters all tend to have the same iron content [174] so globular clusters are thought to be internally chemically homogeneous. The colour-magnitude diagram (CMDs) for most Galactic globular clusters (e.g., M80, Figure 2) also indicate a single stellar population with a distinct main-sequence, main-sequence turn-off, horizontal and giant branch. The single main sequence turn-off in particular indicates that all stars in the cluster have the same age. This leads to a so-called “simple stellar population” model for globular clusters where all stars have the same composition and age and differ only by their masses, which are set by the initial mass function (IMF). This simple picture has been challenged in recent years as observations have shown systematic star-to-star light element variations in globular clusters [174, 376]. Specific effects include different populations in s-process abundances (e.g., [321, 338]), anti-correlations between Na and O (e.g., [320, 321]), variations in CNO elements (e.g., [321, 338]) and even differences in iron abundance (e.g., [320, 321]). The best way of explaining these anomalies so far has been to use self-enrichment models where a single globular cluster experiences several bursts of star formation, each enriched by pollution from the previous generation [81]. How multiple populations affect the CMD of a globular cluster is shown in Figure 3. The importance of these scenarios for relativistic binaries has not yet been explored but if the first and second generation have different IMFs this could affect the number of compact remnants. For this review, we will focus mainly on the case of a simple stellar population but we will discuss details of the multi-generation case further in Section 2.3.

Figure 2
figure 2

Colour-magnitude diagram for M80. Image reproduced from the catalogue of 52 globular clusters (see [413]). The entire catalogue is available at the Padova Globular Cluster Group website [166].

Figure 3
figure 3

CMD for NGC 288 showing evidence of two populations, and models incorporating different metallicities, helium fraction, and age. The best fit model has an age difference of Δ Age = 1.5 Gyr. Image reproduced by permission from Roh et al. [412], copyright by IOP.

The IMF is thought to be universal [34] and is usually taken to be a power-law of the form

$${{dN} \over {dM}} \propto {M^{- {\alpha _i}}},$$
(1)

where M is the mass of a star, N the number of stars, dN/dM is the number of stars in an infinitesimal mass range and αi can take different values for different mass ranges. For values above ∼ 1 M, αi is usually assumed to have a single value, the so-called Salpeter slope, of ∼ 2.35 [416]. There is much more debate about the value of αi in low-mass regime. One solution is to treat the IMF as a broken power-law with a break around 1M, such as that proposed by Kroupa & Weidner [279, 280]:

$$\begin{array}{*{20}c} {{\alpha _0} = 0.3 \pm 0.7,0.01 \leq M/{M_ \odot} < 0.08} \\ {{\alpha _1} = 1.3 \pm 0.5,0.08 \leq M/{M_ \odot} < 0.50} \\ {{\alpha _2} = 2.3 \pm 0.3,0.50 \leq M/{M_ \odot} < 1.00} \\ {{\alpha _3} = 2.3 \pm 0.7,1.00 \leq M/{M_ \odot}}.\quad\quad\\ \end{array}$$
(2)

Another possibility, introduced by Chabrier [70], uses a log-normal distribution of masses below 1 M and a Salpeter slope above. In both cases the power-law strongly favours low masses so stars massive enough to form neutron stars (NSs) and black holes (BHs) will be rare. It is worth noting that both of these IMFs were derived in the context of clustered star formation and young open clusters in the solar neighbourhood. It is often assumed that these IMFs hold for globular clusters as well but, because nearby globular clusters are old, most of the stars above ∼ 0.8 M have already moved off the main-sequence. Those just above 0.8 M will be on the red giant branch and are readily visible in most optical images of globular clusters, such as the image of M80 shown in Figure 4. Those that are more massive have evolved past the giant and horizontal branches and become low-luminosity compact remnants. The fact that the original high-mass population is no longer visible produces an intrinsic uncertainty in our knowledge of the high-mass end of the IMF of globular cluster stars. It is possible that the characteristic mass for star formation is higher in dense, optically thick regions and this would lead to an IMF more biased towards high-mass stars [348]. This in turn would increase the number of massive stellar remnants and could have an effect on the compact binary population.

Figure 4
figure 4

Hubble Space Telescope photograph of the dense globular cluster M80 (NGC 6093).

If we assume that the consequences of multiple stellar generations in globular clusers are negligible, then the picture of Galactic globular cluster stellar populations that emerges from this analysis is of a simple, nearly co-eval, chemically homogeneous set of luminous, low-mass population II stars combined with a low-luminosity population of high-mass stellar remnants. It is interactions with members of this remnant population that will be of particular interest for producing relativistic binaries.

2.2 The structure of globular clusters

Globular clusters are classically modelled as spherical N-body systems. This approximation is relatively good given that the mean ellipticity, e, of Galactic globular clusters is ∼ 0.08 (where e = 1 − (b/a), a is the semi-major axis and b is the semi-minor axis) and that e < 0.3 for all Galactic globular clusters [193]. Globular clusters have a core-halo structure where the core is highly concentrated, reaching densities of up to 106M⊙/pc3, and strongly self-gravitating. The surrounding halo is of much lower density and is less strongly self-gravitating. The structure of a globular cluster can be classified using three basic radii: the core radius (rc), the half-mass radius (rh), and the tidal radius (rt). One definition of the core radius relates it to the central velocity and density through the equation

$${{4\pi G} \over 3}{\rho _c}r_c^2 = v_c^2,$$
(3)

where \(\upsilon _c^2\) is the mean-squared central velocity and is the central density [199]. For an isothermal model this corresponds roughly to the radius at which the density drops to about one third, and thus the surface density to one half, of its central value. Observationally this corresponds to the radius at which the surface brightness drops to half of its central value. For multi-mass systems it is less clear how to arrive at an appropriate theoretical definition of the core radius and more empirical measures, such as the density-weighted average of the distance of each star from the density center are often used [68]. The half-mass radius is simply the radius thatcontains half of the mass of the system. The corresponding observational value is the half-light radius, which contains half the light of the system (the two radii do not necessarily agree). The tidal radius is the radius at which the gravitational field of the host galaxy becomes more important than the self gravity of the star cluster. A classical estimate of this for a circular orbit is given by Spitzer [443] as:

$$r_t^3 = {{{M_{{\rm{GC}}}}} \over {2{M_{{\rm{GAL}}}}}}R_{{\rm{GAL}}}^3,$$
(4)

where MGC is the mass of the globular cluster, MGAL the mass of the galaxy and RGAL the galactocentric radius of the circular orbit. In a time-dependent galactic field the escape process becomes significantly more complicated (see [119] for a more complete theory of tidal escape in a time-varyingexternal potential). For a given cluster, cluster orbit, and galaxy model the tidal radius of the cluster can, in principle, be clearly defined by comparing the effect of the galactic versus globular cluster gravitational field on a test mass. Observationally tidal radii can be difficult to determine due to the low stellar density of globular cluster halos and an imperfect knowledge of the gravitational potential of the host galaxy. Median values for rc, rh, and rt in the Galaxy are ∼ 1 pc, ∼ 3 pc and ∼ 35 pc respectively [199].

There are also two important timescales that characterize globular cluster evolution: the crossing time (tcr) and the relaxation time (trlx). The crossing time is simply the time required for a star traveling at a typical velocity to cross some characteristic cluster radius. Thus, tcrR/v where, for example, or might be typical radii of interest and v could be the velocity dispersion — normally taken to be the root-mean-square velocity and observed to be ∼ 10 km s−1 in Galactic globular clusters [193, 199]. tcr is also, roughly speaking, the orbital timescale for the cluster. For typical values of rh and v, tcr for Galactic globular clusters is on the order of 0.1–1 Myr but is longer at the tidal radius and much shorter in the core.

The relaxation time describes how long it takes for orbits to be significantly altered by stellar encounters. In particular trlx is often defined as the time necessary for the velocity of a star to change by an order of itself [57]. This can be thought of as the time necessary for a cluster to lose the memory of its initial conditions or, more exactly, the time necessary for stellar encounters to transform an arbitrary velocity distribution to a Maxwellian [443]. The relaxation time is related to the number and strength of encounters and, thus, to the number density and energy of a typical star in the cluster. It can be shown that the mean relaxation time in a globular cluster is [57, 443]

$${t_{{\rm{rlx}}}} \propto {N \over {\ln N}}{t_{{\rm{cr}}}}.$$
(5)

Taking a typical value of N = 105 and tcr as before, typical values of trlx are 0.1–1 Gyr. Thus, with ages typically greater than 10 Gyr, Galactic globular clusters are expected to be dynamically relaxed objects. In reality the value of trlx varies significantly within a globular cluster due to the highly inhomogeneous density distribution of the core-halo structure. It is possible for the core of a globular cluster to be fully relaxed while the halo remains un-relaxed after 13 Gyr. By making various approximations for tcr it is possible to relate trlx to local cluster properties. For instance the criterion of Meylan and Heggie [326, 288]:

$${t_{{\rm{rlx}}}} = {{0.065{{\langle {v^2}\rangle}^{3/2}}} \over {\rho \langle m\rangle \ln (\gamma N)}}$$
(6)

relates trlx to the local mass density (ρ), the mass-weighted mean squared velocity (〈v2〉), and the average mass (〈m〉). The criterion by Spitzer [443]:

$${t_{{\rm{rlx}}}} = {{0.138{N^{1/2}}{r^{3/2}}} \over {{{\langle m\rangle}^{1/2}}{G^{1/2}}\ln (\gamma N)}}$$
(7)

relates trlx to a characteristic radius (r), the average mass, and the number of stars in the system. In practise is normally used for the characteristic radius in the Spitzer criterion and trlx is renamed the half-mass relaxation time, trh. The factor ln γN that appears in both definitions is called the Coulomb logarithm and describes the relative effectiveness of small and large angle collisions. It is calculated as an integral over the impact parameters for two-body scattering encounters, b:

$$\ln (\gamma N) \equiv \int\nolimits_{{b_{\min}}}^{{b_{\max}}} {{b \over {db}} = \ln \left({{{{b_{\max}}} \over {{b_{\min}}}}} \right)},$$
(8)

where γ is a constant of order unity. The exact value of γ is a matter of some debate and must be determined empirically. Values in the literature range from 0.02–0.4 depending on the mass distribution of the system [156, 157]. Both definitions of trlx are used extensively in stellar dynamics.

On timescales shorter than a relaxation time interactions between individual stars do not govern the overall evolution of a stellar system and the granularity of the gravitational potential can be ignored. On these timescales the background structure of the cluster can be modelled using a static distribution function, f, that describes the probability of finding a star at a particular location in a 6-D position-velocity phase space. Formally, depends on position, velocity, mass, and time so we have \(f(\vec x,\vec \upsilon, m,t)\). For times less than trlx, however, the evolution of f is described by the collisionless Boltzmann equation:

$${{\partial f} \over {\partial t}} + \vec v \cdot \vec \nabla f - \vec \nabla \phi \cdot {{\partial f} \over {\partial \vec v}} = 0$$
(9)

and the explicit time dependence can be removed: ∂f/∂t = 0. The gravitational potential, ϕ, is given by Poisson’s equation:

$${\vec \nabla ^2}\phi = 4\pi G\rho$$
(10)

and can be calculated at any position by integrating the distribution function over mass and velocity:

$${\vec \nabla ^2}\phi = 4\pi G\int {f(\vec x,\vec v,m){d^3}\vec vdm}.$$
(11)

Solutions to Eq. 9 are often described in terms of the relative energy per unit mass, ɛ ≡ Ψ − v2/2, rather than in terms of the phase-space variables. Here, Ψ = −ϕ + ϕ0 is the relative potential and ϕ0 is defined such that no star has an energy less that zero (f > 0 for ɛ > 0 and f = 0 for ɛ < 0). A simple class of solutions to Eq. 9 are Plummer models [378]:

$$f(\mathcal{E}) = F{\mathcal{E}^{7/2}},$$
(12)

the stellar-dynamical equivalent of an n = 5 ploytrope [199]. Another class of models that admit anisotropy and a distribution in angular momentum, L, are known as King-Mitchie models [264, 327]. Their basic distribution function is:

$$f(\mathcal{E},L) = {\rho _1}{(2\pi {\sigma ^2})^{- 3/2}}\exp ({{- {L^2}} \over {2r_a^2{\sigma ^2}}})\left[ {{e^{\mathcal{E}/{\sigma ^2}}} - 1} \right],\quad \mathcal{E} > 0,$$
(13)

where σ is the velocity dispersion, the anisotropy radius where the velocity distribution changes from nearly isotropic to nearly radial, and ρ1 is a constant related to the density. Although not as well theoretically supported as the single-mass case, King-Mitchie models have been extended to include a spectrum of stellar masses [187] and even external gravitational field [200]. Multi-mass King models in particular are often fit to observed globular cluster cluster surface brightness profiles in order to determine their masses. A good example of the construction of a multi-mass King-Mitchie model is found in the appendix of Miocchi [339].

2.3 The dynamical evolution of globular clusters

Although static models can be used to describe the instantaneous structure of a globular cluster, there is no stable equilibrium for self-gravitating systems [57] and, therefore, their structure changes dramatically over time. Accessible descriptions of globular cluster evolution are given in Hut et al. [241], Meylan and Heggie [326] and Meylan [325] and have also been the subject of several texts (e.g., Spitzer [443] and Heggie and Hut [199]). Here we merely outline some of the more important aspects of globular cluster evolution.

The initial conditions of globular clusters are not well constrained since they seem to form only very early in the history of galaxy formation or in major mergers [150, 61] — both situations quite different from the environment of our Galaxy today. However, it is possible to make some general statements. Like all stars, the stars in globular clusters collapse out of molecular gas. Due to their small age spread and chemical homogeneity, it seems likely that all the stars in a globular cluster formed by the collapse and fragmentation of a single giant molecular cloud [61]. However, the exact details of how this collapse proceeds are unclear. In the classical picture a globular cluster forms in a single collapse event — all stars form rapidly at almost exactly the same time from a globally collapsing cloud of gas and, thus, have almost exactly the same age and chemical composition. Not all of the gas in the collapsing cloud is necessarily converted into stars and the fraction of gas that becomes stars is described through the star formation efficiency (SFE). An SFE of less than 100% implies that not all of the primordial gas forms stars and the resulting globular cluster will be less massive than its parent cloud. Star formation is then terminated and the left-over gas expelled by a combination of ionizing radiation from young stars and energy injection from supernovae. Although more applicable to open clusters, Dale et al. [85] provides some the the most recent results on the details of how this process works in young star clusters. If the star formation efficiency is low, then the amount of mass loss through gas expulsion can leave the cluster out of virial equilibrium and may lead to its immediate dissolution [171, 172, 484], a process called “infant mortality”. It is estimated that more than half of young star clusters in the local universe are destroyed in this manner [172] and even the surviving clusters will lose a large fraction of their stars. Even if the cluster survives the gas expulsion, the rapid change in potential will cause the energy of individual stars to change in a mass-independent manner [57]. This process, called violent relaxation, means that the positions, velocities, and masses of cluster stars will be initially uncorrelated.

As mentioned in Section 2.1, observations are now beginning to challenge some of the details of this simplistic picture of globular cluster formation. Rather than a single burst of star formation the observed abundance anomalies suggest a more drawn-out formation scenario where a first generation of stars forms, evolves and enrichs the cloud with ejecta processed through nuclear burning. Later generations form from the enriched gas and carry the chemical tracers of the pollution. A general description of how such a scenario might proceed can be found in Conroy and Spergel [81] however there is little consensus on the exact details. The initial cloud could suffer a global collapse, experience star formation, have the collapse halted by feedback and then re-collapse to form the second generation. Conversely, different regions of the cloud could collapse into sub-clumps at different times and enrich neighbouring regions with their ejecta. The globular cluster would then assemble from the merger of these sub-clumps. It is not clear what effect, if any, these different formation scenarios will have on the compact binary population of globular clusters. There is some suggestion that the early generations must have an IMF biased towards massive stars in order to explain the observed abundance anomalies (e.g., [91, 92]) and if so this would certainly have implications for the number of compact objects and compact binaries. However, the question is far from resolved. Due to the uncertainties in the multi-generation scenarios and in particular the fact that the details are not necessarily important for compact binary production, we will assume a single star formation event in globular clusters for the purposes of this review.

Equipartition of energy dictates that the most massive stars should have the lowest kinetic energies [57, 443] but violent relaxation leaves the velocities and masses uncorrelated. Therefore, as soon as the residual gas is expelled from a young globular cluster, massive stars with large kinetic energies will start to transfer this energy to low-mass stars through stellar encounters. As the massive stars lose kinetic energy they will sink to the center of the cluster while the low-mass stars gain kinetic energy and move to the halo. This process is known as mass segregation and proceeds on a timescale tmsmi/〈m〉 × trlx [442, 479, 257]. Mass segregation can also occur in other situations, such as around super-massive black holes in galactic centres [136]. Due to mass-loss from stellar evolution, compact remnants rapidly become the most massive objects in globular clusters as the star population ages [419]. Thus they, along with binaries, are strongly affected by mass segregation and are particularly likely to be found in cluster cores [475]. Mass segregation continues until energy equipartition has been achieved. There are, however, initial conditions for which it is impossible to achieve energy equipartition [479, 257] and it is formally impossible to halt mass segregation, leading to a singularity in the core of the cluster. This phenomenon was first noted by Spitzer [442] and is thus called the “Spitzer instability”. In reality, the massive objects in such systems form a strongly interacting subsystem, dynamically decoupled from the rest of the cluster. Due to their high masses, black holes in star clusters are particularly likely to experience the Spitzer instability and this has been the starting point for several investigations of BH binaries in star clusters (e.g. [360]).

The longer-term evolution of star clusters is driven by two-body relaxation, where the orbits of stars are perturbed by encounters with their neighbours. The theory of two-body relaxation was first quantified by Chandrasekhar in 1942 [71]. Two-body relaxation becomes important on timescales longer than the local relaxation time. The evolution of a globular cluster over these timescales can still be described (at least formally) by the Boltzmann equation but with a collisional term added to the right-hand side. Eq. 9 then takes on the form:

$${{\partial f} \over {\partial t}} + \vec v \cdot \vec \nabla f - \vec \nabla \phi \cdot {{\partial f} \over {\partial \vec v}} = \Gamma [f].$$
(14)

Eq. 14 is sometimes called the collisional Boltzmann equation and the term Γ[f] describes the effect of two-body (and in principle higher-order) interactions on the distribution function. Practically speaking, it is not possible to evaluate Γ[f] analytically and various numerical approximations must be employed. Approaches include the Fokker-Planck method, where Γ[f] is approximated in the weak scattering limit by an expansion in powers of the phase-space parameters; the Monte Carlo method, where Γ[f] is approximated by a Monte Carlo selection of weak encounters over a time shorter than the relaxation time; the anisotropic gas model where the cluster is approximated as a self-gravitating gaseous sphere; or direct N-body integration, where rather than solving Eq. 14 the orbits of each star in the cluster are explicitly integrated. Each method has its strengths and weaknesses and will be discussed further in Section 5.

Thermodynamically speaking, strongly self-gravitating systems have a negative heat capacity. This can be seen by relating the kinetic energy of the system to a dynamical temperature [57]:

$${1 \over 2}mv_{{\rm{ave}}}^2 = {3 \over 2}{k_B}T.$$
(15)

This, in combination with the virial theorem, can be used to define a total internal energy for the cluster: 5

$$E = - {3 \over 2}N{k_B}T,$$
(16)

where N is the number of bodies in the system. Finally, this can be used to calculate a heat capacity:

$$C = {{dE} \over {dT}} = - {3 \over 2}N{k_B}.$$
(17)

Since all the constants on the right hand side of Eq. 17 are positive, C is always negative. A negative heat capacity means that heating a self-gravitating system actually causes it to lose energy. For a core-halo star cluster, the core is strongly self-gravitating while the halo is not, so the halo acts as a heat bath for the core. Any perturbation in which the core becomes dynamically hotter than the halo causes energy to flow into the halo. The negative heat capacity means that the core becomes even hotter, increasing the flow of energy to the halo in a runaway process. This causes the core to contract, formally to a singularity. The runaway process is called the gravothermal catastrophe and the consequent shrinking of the core is called core collapse. It affects all self-gravitating systems and was first noted in the context of star clusters by Antonov [18]. In equal-mass systems core collapse will occur after 12–20 trh [443] but may be accelerated in systems with a spectrum of masses due to mass segregation. Core collapse not only appears in analytic models, but has also been found in a variety of numerical simulations such as the model shown in Figure 5 [253]. Furthermore, the Harris catalogue lists several Galactic globular clusters that from their surface brightness profiles are thought to have experienced a core collapse event [193].

Figure 5
figure 5

Lagrange radii indicating the evolution of a Plummer model globular cluster for an N-body simulation and a Monte Carlo simulation. The radii correspond to radii containing 0.35, 1, 3.5, 5, 7, 10, 14, 20, 30, 40, 50, 60, 70, and 80% of the total mass. Image reproduced by permission from Joshi et al. [253], copyright by IOP.

Core collapse can be halted, at least temporarily, by an energy source in the core. For stars (also self-gravitating systems) this energy source is nuclear burning. In star clusters tightly bound binaries perform a similar role. Stars in the core scatter off these binaries and gain kinetic energy at the expense of the orbital energy of the binary. This process reverses the temperature gradient and consequently the gravothermal instability which in turn causes the core to re-expand. The core will then cool again due to the expansion, the temperature gradient will again reverse and the process of core collapse will repeat. These repeated core expansions and contractions are know as gravothermal oscillations [56]. The heating from the binaries is the trigger for the process and, in analogy to nuclear burning in stars, is called “binary burning”. The binaries taking part in binary burning may be either primordial (binaries where the stars were born bound to each other) or dynamically formed by a variety of interactions that will be discussed in Section 4.3. It is worth noting that while tight binaries serve as energy sources for the cluster, loosely bound binaries with orbital velocities below the local velocity dispersion can actually act as energy sinks and may significantly hasten the onset of core collapse [131]. The importance of this effect on the evolution of star clusters remains largely unexplored.

Stars can escape from a star cluster if they gain a velocity greater than the cluster’s escape velocity, \(\left\langle {\upsilon _e^2} \right\rangle = - 4U/M\) where U is the potential energy of the star cluster and M its total mass. Using the virial theorem it is possible to show that \(\left\langle {\upsilon _e^2} \right\rangle = 4\left\langle {{\upsilon ^2}} \right\rangle\) where 〈v2 is the RMS velocity in the cluster [57]. There are two means through which a star can reach the escape velocity. The first is ejection where a single strong interaction, such as occurs during binary burning, gives the star a sufficient velocity impulse to exceed ve. This process is highly stochastic. The second is evaporation where a star reaches escape velocity due to a large number of weak encounters during the relaxation process. Relaxation tends to maintain a local Maxwellian in the velocity distribution and, since a Maxwellian distribution always has a fraction γ = 7.38 × 10−3 stars with v > 2vrms, there are always stars in the cluster with a velocity above the escape velocity. Thus, it is the fate of all star clusters to evaporate. The evaporation time can be estimated as:

$${t_{{\rm{evap}}}} \approx {{{t_{{\rm{rlx}}}}} \over \gamma} = 136{t_{{\rm{rlx}}}}.$$
(18)

This is much longer than a Hubble time so few globular clusters are directly affected by evaporation. Evaporation can, however, be accelerated by the presence of a tidal field. In this case stars whose orbits extend beyond rt are stripped from the cluster and lost. Tidal dynamics are more complicated than a simple radial cutoff would imply and detailed prescriptions taking into account orbital energy and angular momentum as well as a time-varying field for star clusters in elliptical orbits are necessary to capture all of the important processes [457, 119]. It seems likely that the ultimate fate of most globular clusters is destruction due to tidal effects [167, 152]. Since most compact objects are likely to live deep in the core of globular clusters where ejection will normally be due to violent interactions, the details of tidal stripping are unlikely to be critical for the treatment of relativistic binaries in star clusters.

3 Observations

By their very nature, relativistic binaries are compact and faint, so that observations of these systems are difficult unless they are in an interacting phase. Furthermore, observations of binaries that have segregated to the centers of clusters can be complicated by crowding issues. Nonetheless, observations across a broad spectrum using a variety of ground- and space-based telescopes have revealed populations of binaries and their tracers in many Milky Way and extra-galactic globular clusters.

In a CMD of a globular cluster, some stars can be found on the main sequence above and to the left of the turn-off (see Figure 6). These stars are known as blue stragglers. As their name suggests, if these stars are coeval with the stellar population of the cluster, they are too hot and too massive and should have evolved off the main sequence. They are thought to have been rejuvenated through mass transfer, merger, or direct collision of binaries [129, 97, 377, 268, 369]. As such, they can be interpreted as tracers of binary interactions within globular clusters. Recent observations of the blue straggler populations of 13 globular clusters indicate a correlation between the specific frequency of blue stragglers and the binary fraction in the globular cluster [441]. This supports observations which also seem to suggest that binary coalescences are the dominant formation mechanism for blue stragglers in globular clusters [291]. However, many of the populations that have been observed exhibit a bimodal radial distribution [127, 488, 123, 414, 478, 286, 87, 42, 86] for which the inner population can be interpreted as tracing dynamical interactions and the outer population as representative of the primordial binary population.

Figure 6
figure 6

CMD of NGC 6205 from the Hubble Space Telescope Advanced Detector for Surveys. Different populations are identified as horizontal branch (HB), red giant branch (RGB), main sequence turn off (MSTO), and blue stragglers (BS). Image reproduced by permission from Leigh, Sills, & Knigge [292].

Blue stragglers are some of the most visible and populous evidence of the dynamical interactions that can also give rise to relativistic binaries. Current observations that have been revealing the blue straggler population and their variable counterparts (SX Phe stars) are the ACS Survey of Galactic Globular Clusters [418] and the Cluster AgeS Experiment (CASE) [76]. For a good description of surveys that use far-ultraviolet in detecting these objects, see Knigge [266]. For somewhat older but still valuable reviews on the implications of blue stragglers on the dynamics of globular clusters, see Hut [236] and Bailyn [25].

The globular cluster population of white dwarfs can be used to determine the ages of globular clusters [343], and so they have been the focus of targeted searches despite the fact that they are arguably the faintest electromagnetically detectable objects in globular clusters. These searches have yielded large numbers of globular cluster white dwarfs. For example, a recent search of ω Centauri has revealed over 2000 white dwarfs [344], while Hansen et al. [191] have detected 222 white dwarfs in M4. Deep ACS observations of NGC 6397 [411] have identified a substantial population of approximately 150 white dwarfs [453]. Dieball et al. [107] have found approximately 30 white dwarfs and about 60 probably cataclysmic variables and white dwarf binaries in M80. In general, however, these searches uncover single white dwarfs. Optical detection of white dwarfs in binary systems tends to rely on properties of the accretion process related to the binary type. Therefore, searches for cataclysmic variables (CVs, see Section 3.1) generally focus on low-luminosity X-ray sources [251, 179, 470] and on ultraviolet-excess stars [104, 105, 177, 269, 319], but these systems are usually a white dwarf accreting from a low mass star. The class of “non-flickerers” which have been detected recently [82, 459] have been explained as He white dwarfs in binaries containing dark CO white dwarfs [116, 180, 190, 454].

Pulsars, although easily seen in radio, are difficult to detect when they occur in hard binaries, due to the Doppler shift of the pulse intervals. Thanks to an improved technique known as an “acceleration search” [328], which assumes a constant acceleration of the pulsar during the observation period, more short-orbital-period binary pulsars are being discovered [62, 64, 88, 90, 133, 140, 400]. For a good review and description of this technique, see Lorimer [296]. The progenitors of the ultracompact millisecond pulsars (MSPs) are thought to pass through a low-mass X-ray binary phase [100, 179, 246, 401, 405] (LMXBs, see Section 3.2). These systems are very bright and when they are in an active state, they can be seen anywhere in the Galactic globular cluster system. There are, however, several additional LMXBs that are currently quiescent [179, 207, 471, 185, 184] (qLMXBs). As these systems turn on, they can be added to the list of known LMXBs, which is currently at 15 [202]. Additional evidence of a binary spin-up phase for MSPs comes from measurements of their masses, which indicate a substantial mass-transfer phase during the spin-up. Several observed globular cluster MSPs in binary systems are seen to have masses above the canonical mass of 1.4 M [135].

Although there are many theoretical predictions of the existence of black holes in globular clusters (see, e.g., [334, 388, 333, 91]), these have generally predicted that only a small number would be retained in the cluster. Prior to the recent discovery of two black hole binaries in M22 [451], there have been only a few observational hints of their existence. Measurements of the kinematics of the cores of M15 [145, 183], NGC 6752 [113], and ω Centauri [357] provide some suggestions of a large, compact mass. One proposed explanation would be a large black hole of ∼ 103 M in the cluster core. Black holes in this mass range are much larger than stellar mass black holes (10–100 M) but much smaller than the supermassive black holes (≳ 106) found in galactic centers. As such they are called intermediate mass black holes (IMBH) and may potentially form through mergers of stars or stellar-mass black holes in globular clusters. However, these observations can also be explained without requiring an IMBH [314, 367] and the existence of such objects is questionable. Furthermore, ω Centauri may be the stripped core of a dwarf galaxy. Therefore, any large black hole in its center may have arisen by other, non-dynamical means. VLA observations of M80, M62, and M15 do not indicate any significant evidence of radio emission, which can be used to place some limits on the likelihood of an accreting black hole in these clusters. While radio observations have provided the strongest limit on the mass of any IMBH in M15, M19, and M22 [452]. However, uncertainties in the expected gas density makes it difficult to place any upper limits on a black hole mass [32]. Recent observations of the kinematics of NGC 2808 [302] and NGC 6388 [303] have placed upper bounds of ∼ 104 M on any IMBH in their cores. The unusual millisecond pulsar in the outskirts of NGC 6752 has also been argued to be the result of a dynamical interaction with a possible binary intermediate mass black hole in the core [80]. If the velocity dispersion in globular clusters follows the same correlation to black hole mass as in galactic bulges, then there may be black holes with masses in the range 50–103 M in many globular clusters [493, 297]. Recent Hubble Space Telescope observations of the stellar dynamics in the core of 47 Tuc have placed an upper bound of 1500 M for any intermediate mass black hole in this cluster [322]. Stellar mass black hole binaries may also be visible as low luminosity X-ray sources, but if they are formed in exchange interactions, they will have very low duty cycles and hence are unlikely to be seen [254]. For a good recent review on neutron stars and black holes in globular clusters, see Rasio et al. [406]. Heinke [202] has an excellent review of X-ray observations of Galactic globular clusters, including CVs, LMXBs and qLMXBs, MSPs, and other active binaries.

Recent observations and catalogs of known binaries are presented in the following sections.

3.1 Cataclysmic variables

Cataclysmic variables (CVs) are white dwarfs accreting matter from a companion that is usually a dwarf star or another white dwarf. They have been detected in globular clusters through identification of the white dwarf itself or through evidence of the accretion process. White dwarfs managed to avoid detection until observations with the Hubble Space Telescope revealed photometric sequences in several globular clusters [83, 82, 366, 408, 409, 410, 459, 191]. Spectral identification of white dwarfs in globular clusters has begun both from the ground with the VLT [342, 343] and in space with the Hubble Space Telescope [82, 116, 459, 344]. With spectral identification, it will be possible to identify those white dwarfs in hard binaries through Doppler shifts in the Hβ line. This approach has promise for detecting a large number of the expected double white dwarf binaries in globular clusters. Photometry has also begun to reveal orbital periods [349, 115, 255] of CVs in globular clusters.

Accretion onto the white dwarf builds up a layer of hydrogen-rich material that may eventually lead to a dwarf nova outburst as the layer undergoes a thermonuclear runaway. These dwarf novae cause an increase in luminosity of a factor of a few to a few hundred. Identifications of globular cluster CVs have been made through such outbursts in the cores of M5 [319], 47 Tuc [365], NGC 6624 [431], M15 [428], and M22 [16, 58]. With the exception of V101 in M5 [319], original searches for dwarf novae performed with ground based telescopes proved unsuccessful. This is primarily due to the fact that crowding obscured potential dwarf novae up to several core radii outside the center of the cluster [425, 427]. Since binaries tend to settle into the core, it is not surprising that none were found significantly outside of the core. Subsequent searches using the improved resolution of the Hubble Space Telescope eventually revealed a few dwarf novae close to the cores of selected globular clusters [424, 426, 431, 428, 16]. To date, there have been 14 found [374, 423], using the Hubble Space Telescope and Las Campanas Observatory (CASE).

A more productive approach has been to look for direct evidence of the accretion around the white dwarf. This can be in the form of excess UV emission and strong Hα emission [124, 178, 269, 270, 104] from the accretion disk. This technique has resulted in the discovery of candidate CVs in 47 Tuc [124, 269], M92 [126], NGC 2808 [104], NGC 6397 [82, 116, 459], and NGC 6712 [125]. The accretion disk can also be discerned by very soft X-ray emissions. These low luminosity X-ray binaries are characterized by a luminosity LX < 1034.5 erg/s, which distinguishes them from the low-mass X-ray binaries with LX > 1036 erg/s. Initial explanations of these objects focused on accreting white dwarfs [24], and a significant fraction of them are probably CVs [471, 480]. There have been 10 identified candidate CVs in NGC 6752 [382], 15 in NGC 6397 [78], 19 in NGC 6440 [381], 2 in ω Cen [147], 5 in Terzan 5 [204], 22 in 47 Tuc [114], 5 in M80 [206], 7 in M54 [399], 2–5 in NGC 288 [273], 4 in M30 [301], 4 in NGC 2808 [422], 1 in M71 [227], and 1 in M4 [33]. However, some of the more energetic sources may be LMXBs in quiescence [471], or even candidate QSO sources [33].

The state of the field at this time is one of rapid change as Chandra results come in and optical counterparts are found for the new X-ray sources. A living catalog of CVs was created by Downes et al. [108], but it has been allowed to lapse and is now archival [109, 110]. It may still be the best source for confirmed CVs in globular clusters up to 2006.

3.2 Low-mass X-ray binaries

The X-ray luminosities of low-mass X-ray binaries are in the range Lx ∼ 1036–1038 erg/s. The upper limit is close to the Eddington limit for accretion onto a neutron star, so these systems must contain an accreting neutron star or black hole. All of the LMXBs in globular clusters contain an accreting neutron star as they also exhibit X-ray bursts, indicating thermonuclear flashes on the surface of the neutron star [251]. Compared with ∼ 100 such systems in the galaxy, there are 15 LMXBs known in globular clusters. The globular cluster system contains roughly 0.1% of the mass of the galaxy and roughly 10% of the LMXBs. Thus, LMXBs are substantially over-represented in globular clusters. Within the Galactic globular cluster system, high metallicity clusters are more likely to host an LMXB, as noted by Grindlay in 1993 [176] and Bellazzini in 1995 [50].

Because these systems are so bright in X-rays, the globular cluster population is completely known — we expect new LMXBs to be discovered in the globular cluster system only as quiescent (qLMXBs) and transient systems become active. The 15 sources are in 12 separate clusters. Five have orbital periods greater than a few hours, six ultracompact systems have measured orbital periods less than one hour, and four have undetermined orbital periods. A member of the ultracompact group, 4U 1820-30 (X1820-303) in the globular cluster NGC 6624, has an orbital period of 11 minutes [447]. This is one of the shortest known orbital periods of any binary and most certainly indicates a degenerate companion. The orbital period, X-ray luminosity, and host globular clusters for these systems are given in Table 1.

Table 1 Low-mass X-ray binaries in globular clusters. Host clusters and LMXB properties.

The improved resolution of Chandra allows for the possibility of identifying optical counterparts to LMXBs. If an optical counterpart can be found, a number of additional properties and constraints for these objects can be determined through observations in other wavelengths. In particular, the orbital parameters and the nature of the secondary can be determined. So far, optical counterparts have been found for X0512-401 in NGC 1851 [224], X1745−203 in NGC 6440 [473], X1746−370 in NGC 6441 [99], X1830-303 in NGC 6624 [265], X1832-330 in NGC 6652 [100, 203], X1850-087 in NGC 6712 [84, 26, 354], X1745−248 in Terzan 5 [204], and both LMXBs in NGC 7078 [20, 485]. Continued X-ray observations will also further elucidate the nature of these systems [347].

The 15 bright LMXBs are thought to be active members of a larger population of lower luminosity quiescent low mass X-ray binaries (qLMXBs) [486]. Early searches performed with ROSAT data (which had a detection limit of 1031 erg/s) revealed roughly 30 sources in 19 globular clusters [251]. A more recent census of the ROSAT low luminosity X-ray sources, published by Verbunt [470], lists 26 such sources that are probably related to globular clusters. Recent observations with the improved angular resolution of Chandra have begun to uncover numerous low luminosity X-ray candidates for CVs [179, 180, 203, 225, 204, 206, 114, 115, 147, 382, 381]. For a reasonably complete discussion of recent observations of qLMXBs in globular clusters, see Verbunt and Lewin [471] or Webb and Barret [480] and references therein. Table 2 lists the 27 qLMXBs known to date.

Table 2 Quiescent Low-mass X-ray binaries in globular clusters.

3.3 Millisecond pulsars

The population of known millisecond pulsars (MSPs) is one of the fastest growing populations of relativistic binaries in globular clusters. Several ongoing searches continue to reveal millisecond pulsars in a number of globular clusters. Previous searches have used deep multi-frequency imaging to estimate the population of pulsars in globular clusters [140]. In this approach, the expected number of pulsars beaming toward the earth, Npuls, is determined by the total radio luminosity observed when the radio beam width is comparable in diameter to the core of the cluster. If the minimum pulsar luminosity is Lmin and the total luminosity observed is Ltot, then, with simple assumptions on the neutron star luminosity function,

$${N_{{\rm{puls}}}} = {{{L_{{\rm{tot}}}}} \over {{L_{\min}}\ln ({L_{{\rm{tot}}}}/{L_{{\rm{min}}}})}}.$$
(19)

In their observations of 7 globular clusters, Fruchter and Goss have recovered previously known pulsars in NGC 6440, NGC 6539, NGC 6624, and 47 Tuc [140]. Their estimates based on Eq. (19) give evidence of a population of between 60 and 200 previously unknown pulsars in Terzan 5, and about 15 each in Liller 1 and NGC 6544 [140]. Additional Fermi LAT work uses γ-ray emission to estimate total numbers of ≃ 2600–4700 MSPs in the Galactic globular cluster system, which corresponds to estimates from encounter rates [8].

Current searches include the following: Arecibo, which is searching over 22 globular clusters [214]; Green Bank Telescope (GBT), which is working both alone and in conjunction with Arecibo [249, 214]; the Giant Metrewave Radio Telescope (GMRT), which is searching over about 10 globular clusters [134]; and Parkes, which is searching over 60 globular clusters [89]. Although these searches have been quite successful, they are still subject to certain selection effects such as distance, dispersion measure, and acceleration in compact binaries [63]. For an excellent review of the properties of all pulsars in globular clusters, see the review by Camilo and Rasio [63] and references therein. The properties of known pulsars in binary systems with orbital period less than one day are listed in Table 3, which has been extracted from the online catalog maintained by Freire [397].

Table 3 Short-orbital-period binary millisecond pulsars in globular clusters. Host clusters and orbital properties. See http://www.naic.edu/∼pfreire/GCpsr.html for references to each pulsar.

With the ongoing searches, it can be reasonably expected that the number of millisecond pulsars in binary systems in globular clusters will continue to grow in the coming years.

3.4 Black holes

Radio observations of M22 have revealed two stellar mass black holes with likely low-mass stellar companions [451]. These objects have flat radio spectra, but no observed X-ray emission. As such, they are likely to be undergoing very low level accretion from their companions. Mass estimates put these black holes at ∼ 10–20 M. There have been no confirmed observations of X-ray black hole binaries in Galactic globular clusters. All of the Galactic globular cluster high luminosity LMXBs exhibit the X-ray variability that is indicative of nuclear burning on the surface of a neutron star. It is possible that some of the recently discovered low luminosity LMXBs may house black holes instead of neutron stars [471], it is more likely that they are simply unusual neutron star LMXBs in quiescence [486]. Finally, there is very circumstantial evidence for the possible existence of an intermediate mass black hole (IMBH) binary in NGC 6752 based upon an analysis of the MSP binary PSR A [79, 80]. Further evidence for black holes in globular clusters will probably have to wait for the next generation of gravitational wave detectors.

3.5 Extragalactic globular clusters

Recently, the X-ray observatories XMM-Newton, Chandra, and Swift have been combined with HST observations to search for extragalactic globular cluster binaries. These have yielded a number of systems found in M31, M104, NGC 1399, NGC 3379, NGC 4278, NGC 4472, and NGC 4697.

Three super soft X-ray sources have been found in M31 globular clusters [212, 375, 213]. Super soft sources have little or no radiation at energies above 1 keV and have a blackbody spectrum corresponding to temperatures between 3 × 104 K to 2 × 105 K. These are assumed to be from classical novae, and two have been associated with optical novae [212, 213]. Stiele et al. [448] have found 36 LMXBs and 17 candidate LMXBs that are co-located with M31 globular clusters, while Peacock et al. [368] claim 41 LMXBs in 11% of the M31 globular clusters. Barnard et al. [31] have found 4 black hole X-ray binaries associated with M31 globulars. Finally, there is the suspected intermediate mass black hole in G1 [17, 102].

The X-ray luminosity function (XLF) of the population of LMXB’s in individual galaxies has been found to obey a power law with a slope in the range of 1.8 to 2.2 [164, 258]. When the XLF’s for many galaxies are combined, a broken power law provides a better fit [258]. At the break luminosity of ∼ 5 × 1038 erg, the XLF steepens to a slope consistent with ∼ 2.8. The break luminosity is near the Eddington luminosity for a 3 M accretor. The XLFs for several individual galaxies are shown in Figure 7 and the combined XLFs from Kim and Fabbiano [258] are shown in Figure 8. When the globular cluster XLF is separated from the field XLF, we see a dramatic drop in the number of low luminosity LMXBs for the globular cluster population, indicative of a different formation mechanism [477, 492].

Figure 7
figure 7

X-ray luminosity functions for several individual galaxies. Each can be fit will with a single power law. Image reproduced by permission from Kim & Fabbiano [258], copyright by IOP.

Figure 8
figure 8

X-ray luminosity function for several galaxies combined together. On the left is the best fit to a single power law, while a broken power law fit is shown in the right. Image reproduced by permission from Kim & Fabbiano [258], copyright by IOP.

Searches for bright X-ray sources in other nearby galaxies have yielded LMXB’s thought to contain black holes based on their luminosity. There are 2 in NGC 1399 [432, 243, 313], and Paolillo et al. [364] estimate that 65% of globular clusters in NGC 1399 host LMXBs. Brassington et al. [60] have found 3 LMXBs in globular clusters associated with NGC 3379, of which 1 is a presumed black hole. Fabbiano et al. [120] have found 4 LMXBs in globular clusters associated with NGC 4278, all of which are bright enough to be black holes. Maccarone and collaborators [308, 310] have found 2 LMXBs containing black holes in NGC 4472 globular clusters. Kim et al. [259] claim 75 globular clusters LMXBs in NGC 3379, NGC 4278, and NGC 4697. Finally, Li et al. [293] have found 41 X-ray sources that are presumably LMXBs in globular clusters in M104. The extragalactic LMXBs have been found to follow the same trend toward over-representation in red (or metal rich) globular clusters [440, 282, 307, 261, 283], as can be seen in Figure 9.

Figure 9
figure 9

GC magnitudes (Mz) vs. GC colors (g − z), with integrated histograms of the properties above and to the right. GCs unmatched to LMXBs are indicated by small black dots and open black histograms (scaled down by a factor of 10). Blue GCs with LMXBs are indicated by filled blue squares and histograms. Red GCs with LMXBs are indicated by filled red circles and histograms. The histograms of the GCs with LMXBs are stacked on each other. GCs that are redder and brighter are more likely to contain LMXBs. Image reproduced by permission from Sivakoff et al. [440], copyright by IOP.

4 Relativistic Binaries

We define relativistic binaries to be binary systems containing two degenerate or collapsed objects and an orbital period such that they will be brought into contact within a Hubble time. Note that this definition in intentionally vague and includes binaries which are already in contact as well as systems that will never have relativistic orbital velocities (such as double white dwarfs). Frequently, in the process of becoming a relativistic binary, a binary will exist with a single degenerate or collapse object and a normal star. These systems are tracers of relativistic binaries. Outside of dense stellar clusters, most relativistic binary systems arise from isolated primordial binary systems whose evolution drives them to tight, ultracompact orbits. Dynamical processes in globular clusters can drive wide binary systems toward short orbital periods and can also insert degenerate or collapsed stars into relativistic orbits with other stars. Before addressing specific evolutionary scenarios, we will present the generic features of binary evolution that lead to the formation of relativistic binaries.

4.1 Binary evolution

The evolution of an isolated binary system of two main-sequence stars can significantly affect the evolution of both component stars if the orbital separation is sufficiently small. For low mass stars (M ≲ 1.25 M), tidal dissipation in the convective envelope will circularize the binary on a time scale given by \({\tau _{{\rm{circ}}}} \propto P_{{\rm{orb}}}^{10/3}\) [168, 489, 490]. For such stars, this implies that nearly all binaries with orbital periods less than 10 days will be circularized by the early main sequence phase. For high mass stars with radiative envelopes, radiative damping is the dominant mechanism for circularization. The circularization time scale for high mass stars is τcirc ∝ (R/a)21/2 where R is the stellar radius and a is the relative semimajor axis [355]. All high mass binaries with R/a > 0.25 will be circularized. Both stars start in the main sequence with the mass of the primary Mp and the mass of the secondary Ms defined, such that MpMs. The binary system is described by the orbital separation r, and the mass ratio of the components qMs/Mp. The gravitational potential of the binary system is described by the Roche model where each star dominates the gravitational potential inside regions called Roche lobes. The two Roche lobes meet at the inner Lagrange point along the line joining the two stars. Figure 10 shows equipotential surfaces in the orbital plane for a binary with q = 0.4. If the volume of either star exceeds the effective volume of its Roche lobe, then it is said to fill its Roche lobe. Matter will stream from a Roche lobe filling star through the inner Lagrange point to the other star in a process known as Roche lobe overflow (RLOF). This mass transfer affects both the evolution of the components of the binary as well as the binary properties such as orbital period and eccentricity.

Figure 10
figure 10

Cross section of equipotential surfaces in the orbital plane of a binary with q = 0.4. The values of the potential surfaces are 5.0, 3.9075, 3.8, 3.559, 3.2, 3.0, and 2.8. The units have been normalized to the orbital separation, so a = 1.

Roche lobe overflow can be triggered by the evolution of the binary properties or by evolution of the component stars. On the one hand, the orbital separation of the binary can change so that the Roche lobe can shrink to within the surface of one of the stars. On the other hand, stellar evolution may eventually cause one of the stars to expand and fill its Roche lobe. When both stars in the binary are main-sequence stars, the latter process is more common. Since the more massive star will evolve first, it will be the first to expand and fill its Roche lobe. At this stage, the mass exchange can be conservative (no mass is lost from the binary) or non-conservative (mass is lost). Depending on the details of the mass exchange and the evolutionary stage of the mass-losing star there are several outcomes that will lead to the formation of a relativistic binary. The primary star can lose its envelope, revealing its degenerate core as either a helium, carbon-oxygen, or oxygenneon white dwarf; it can explode as a supernova, leaving behind a neutron star or a black hole; or it can simply lose mass to the secondary so that they change roles. Barring disruption of the binary, its evolution will then continue. In most outcomes, the secondary is now the more massive of the two stars and it may evolve off the main sequence to fill its Roche lobe. The secondary can then initiate mass transfer or mass loss with the result that the secondary can also become a white dwarf, neutron star, or black hole.

The relativistic binaries that result from this process fall into a number of observable categories. A WD-MS or WD-WD binary may eventually become a cataclysmic variable once the white dwarf begins to accrete material from its companion. If the companion is a main-sequence star, RLOF can be triggered by the evolution of the companion or more commonly by a process known as magnetic braking. If the companion is a low mass star with a convective envelope (M ≲ 1.5 M) then the binary can lose angular momentum through its wind in a process proposed by Verbunt and Zwaan [474] in the context of LMXBs. The general mechanism of magnetic braking involves the removal of large amounts of angular momentum through an ionized stellar wind that is forced to co-rotate with the magnetic field of the star out to large distances. Stars with convective envelopes can sustain the enhanced magnetic activity necessary to support this mechanism. Further details of the magnetic braking process can be found in Postnov and Yungelson [392] and Knigge, Baraffe, and Patterson [267]. If the companion is another white dwarf, then RLOF is triggered by the gradual shrinking of the orbit through the emission of gravitational radiation. Some WD-WD cataclysmic variables are also known as AM CVn stars if they exhibit strong He lines. If the total mass of the WD-WD binary is above the Chandrasekhar mass (∼ 1.4, M), the system may be a double degenerate progenitor to a Type Ia supernova.

The orbit of a NS-MS or NS-WD binary will shrink due to the emission of gravitational radiation. At the onset of RLOF, the binary will become either a low-mass X-ray binary (if the donor star has a lower mass than the accretor, typically a white dwarf or main sequence star with M < 2 M), or a high-mass X-ray binary (if the donor is the more massive component). These objects may further evolve to become millisecond pulsars if the neutron star is spun up during the X-ray binary phase [96, 405]. A NS-NS binary will remain virtually invisible unless one of the neutron stars is observable as a pulsar. A BH-MS or BH-WD binary may also become a low- or high-mass X-ray binary. If the neutron star is observable as a pulsar, a BH-NS binary will appear as a binary pulsar. BH-BH binaries will be invisible unless they accrete matter from the interstellar medium. A comprehensive table of close binary types that can be observed in electromagnetic radiation can be found in Hilditch [215].

The type of binary that emerges depends upon the orbital separation and the masses of the component stars. During the evolution of a 10 M star, the radius will slowly increase by a factor of about two as the star progresses from zero age main sequence to terminal age main sequence. The radius will then increase by about another factor of 50 as the star transitions to the red giant phase, and an additional factor of 10 during the transition to the red supergiant phase. These last two increases in size occur very quickly compared with the slow increase during the main-sequence evolution [372]. Depending upon the orbital separation, the onset of RLOF can occur any time during the evolution of the star. Mass transfer can be divided into three cases related to the timing of the onset of RLOF (see Hilditch [215] for more details):

  • Case A: If the orbital separation is small enough (such that the period is usually a few days), the star can fill its Roche lobe during its slow expansion through the main-sequence phase while still burning hydrogen in its core.

  • Case B: If the orbital period is less than about 100 days, but longer than a few days, the star will fill its Roche lobe during the rapid expansion to a red giant with a helium core. If the helium core ignites during this phase and the transfer is interrupted, the mass transfer is case BB.

  • Case C: If the orbital period is above 100 days, the star can evolve to the red supergiant phase before it fills its Roche lobe. In this case, the star may have a CO or ONe core.

The typical evolution of the radius for a low metallicity star is shown in Figure 11. Case A mass transfer occurs during the slow growth, Case B during the first rapid expansion, and Case C during the final expansion phase. The nature of the remnant depends upon the state of the primary during the onset of RLOF and the orbital properties of the resultant binary depend upon the details of the mass transfer.

Figure 11
figure 11

Evolution of the radius for a 10 M star with a metallicity of Z = 0.001. Image reproduced by permission from Pfahl et al. [372], copyright by IOP.

4.2 Mass transfer

Although there are still many unanswered theoretical questions about the nature of the mass transfer phase, the basic properties of the evolution of a binary due to mass transfer can easily be described. The rate at which a star can adjust to changes in its mass is governed by three time scales. The dynamical time scale results from the adiabatic response of the star to restore hydrostatic equilibrium, and can be approximated by the free fall time across the radius of the star,

$${t_{{\rm{dyn}}}} \simeq {\left({{{2{R^3}} \over {GM}}} \right)^{1/2}}\sim 40{\left[ {{{\left({{R \over {{R_ \odot}}}} \right)}^3}{{{M_ \odot}} \over M}} \right]^{1/2}}\min,$$
(20)

where M and R are the mass and radius of the star. The thermal equilibrium of the star is restored over a longer period given by the thermal time scale

$${t_{{\rm{th}}}} \simeq {{G{M^2}} \over {RL}}\sim 3 \times {10^7}{\left({{M \over {{M_ \odot}}}} \right)^2}{{{R_ \odot}} \over R}{{{L_ \odot}} \over L}{\rm{yr}},$$
(21)

where L is the luminosity of the star. Finally, the main-sequence lifetime of the star itself provides a third time scale, which is also known as the nuclear time scale:

$${t_{{\rm{nuc}}}}\sim 7 \times {10^9}{M \over {{M_ \odot}}}{{{L_ \odot}} \over L}{\rm{yr}}.$$
(22)

The rate of mass transfer/loss from the Roche lobe filling star is governed by how the star’s radius changes in response to changes in its mass. Hjellming and Webbink [221] describe these changes and the response of the Roche lobe to mass changes in the binary using the radius-mass exponents, ζd ln R/d ln M, for each of the three processes described in Eqs. (20, 21, 22) and defining

$${\zeta _{\rm{L}}} = (1 + q){{d\ln {R_{\rm{L}}}} \over {d\ln q}}$$
(23)

for the Roche lobe radius-mass exponent. If ζL > ζdyn, the star cannot adjust to the Roche lobe, then the mass transfer takes place on a dynamical time scale and is limited only by the rate at which material can stream through the inner Lagrange point. If ζdyn > ζl > ζth, then the mass transfer rate is governed by the slow expansion of the star as it relaxes toward thermal equilibrium, and it occurs on a thermal time scale. If both ζdyn and ζth are greater than ζL, then the mass loss is driven either by stellar evolution processes or by the gradual shrinkage of the orbit due to the emission of gravitational radiation. The time scale for both of these processes is comparable to the nuclear time scale. A good analysis of mass transfer in cataclysmic variables can be found in King et al. [263] and Knigge, Baraffe, and Patterson [267].

Conservative mass transfer occurs when there is no mass loss from the system, and therefore all mass lost from one star is accreted by the other star. During conservative mass transfer, the orbital elements of the binary can change. Consider a system with total mass M = M1 + M2 and semi-major axis a. The total orbital angular momentum

$$J = {\left[ {{{GM_1^2M_2^2a} \over M}} \right]^{1/2}}$$
(24)

is a constant, and we can write a ∝ (M1M2)−2. Using Kepler’s third law and denoting the initial values by a subscript i, we find:

$${P \over {{P_i}}} = {\left[ {{{{M_{1i}}{M_{2i}}} \over {{M_1}{M_2}}}} \right]^3}.$$
(25)

Differentiating Eq. (25) and noting that conservative mass transfer requires \({{\dot M}_1} = - {{\dot M}_2}\) gives:

$${{\dot P} \over P} = {{3{{\dot M}_1}({M_1} - {M_2})} \over {{M_1}{M_2}}}.$$
(26)

Note that if the more massive star loses mass, then the orbital period decreases and the orbit shrinks. If the less massive star is the donor, then the orbit expands. Usually, the initial phase of RLOF takes place as the more massive star evolves. As a consequence, the orbit of the binary will shrink, driving the binary to a more compact orbit.

In non-conservative mass transfer, both mass and angular momentum can be removed from the system. There are two basic non-conservative processes that are important for the formation of relativistic binaries — the common-envelope process and the supernova explosion of one component of the binary. The result of the first process is often a short-period, circularized binary containing a white dwarf. Although the most common outcome of the second process is the disruption of the binary, occasionally this process will result in an eccentric binary containing a neutron star or a black hole.

Common envelope scenarios result when one component of the binary expands so rapidly that the mass transfer is unstable and the companion becomes engulfed by the donor star. This can happen if the mass transfer rate is so great that it exceeds the Eddington mass accretion rate of the accretor, or when the donor expands past the outer Lagrange point [361, 481]. The companion can then mechanically eject the envelope of the donor star. There are two proposed approaches to determining the outcome of the process of ejection of the envelope.

The most commonly used approach is the “α-prescription” of Webbink [482], in which the energy required to eject the envelope comes from the orbital energy of the binary and thus the orbit shrinks. The efficiency of this process determines the final orbital period after the common envelope phase. This is described by the efficiency parameter

$${\alpha _{{\rm{CE}}}} = {{\Delta {E_{{\rm{bind}}}}} \over {\Delta {E_{{\rm{orb}}}}}},$$
(27)

where ΔEbind bind is the binding energy of the mass stripped from the envelope and ΔEorb is the change in the orbital energy of the binary. The result of the process is the exposed degenerate core of the donor star in a tight, circular orbit with the companion. This process can result in a double degenerate binary if the process is repeated twice or if the companion has already evolved to a white dwarf through some other process such as standard stellar evolution. If enough orbital energy is lost it can also lead to a merger of the binary components. A brief description of the process is outlined by Webbink [482], and a discussion of the factors involved in determining αCE is presented in Sandquist et al. [417].

The other approach that has been suggested is the “γ-prescription” of Nelemans et al. [351]. In this approach, other energy sources such as tidal heating or luminosity of the donor may assist in the unbinding of the envelope. The material lost by the envelope is assumed to carry away angular momentum and reduce the total angular momentum of the system, such that:

$${J_f} = {J_i}\left({1 - \gamma {{\Delta M} \over {{M_{{\rm{tot}}}}}}} \right),$$
(28)

where Ji and Jf are the initial and final angular momenta of the binary, Mtot is the total mass of the binary prior to mass loss, and ΔM is the mass lost in the ejection of the envelope. for most reasonable models of post common envelope binaries, the value of γ was found to lie in the range of 1.5 ≤ γ ≤ 1.7 and this was interpreted to mean that γ was a more constrained parameter to determine outcome of a common envelope phase. Recently, however, Webbink showed that for reasonable initial and final masses, γ will lie within this range if the angular momentum loss is required to lie between the minimum given by the Jeans mode mass loss and total loss of angular momentum of the system [483]. This has been borne out by observation [494].

The effect on a binary of mass loss due to a supernova can be quite drastic. Following Padmanabhan [363], this process is outlined using the example of a binary in a circular orbit with a semi-major axis a. Let v be the velocity of one component of the binary relative to the other component. The initial energy of the binary is given by

$$E = {1 \over 2}\left({{{{M_1}{M_2}} \over {{M_1} + {M_2}}}} \right){v^2} - {{G{M_1}{M_2}} \over a} = - {{G{M_1}{M_2}} \over {2a}}.$$
(29)

Following the supernova explosion of M1, the expanding mass shell will quickly cross the orbit of M2, decreasing the gravitational force acting on the secondary. The new energy of the binary is then

$$E^{\prime} = {1 \over 2}\left({{{{M_{{\rm{NS}}}}{M_2}} \over {{M_{{\rm{NS}}}} + {M_2}}}} \right){v^2} - {{G{M_{{\rm{NS}}}}{M_2}} \over a},$$
(30)

where MNS is the mass of the remnant neutron star. We have assumed here that the passage of the mass shell by the secondary has negligible effect on its velocity (a safe assumption, see Pfahl et al. [372] for a discussion), and that the primary has received no kick from the supernova (not necessarily a safe assumption, but see Davies and Hansen [96] or Pfahl et al. [373] for an application to globular cluster binaries). Since we have assumed that the instantaneous velocities of both components have not been affected, we can replace them by v2 = G (M1 + M2)/a, and so

$$E^{\prime} = {{G{M_{{\rm{NS}}}}{M_2}} \over {2a}}\left({{{{M_{\rm{1}}} + {M_2}} \over {{M_{{\rm{NS}}}} + {M_2}}} - 2} \right).$$
(31)

Note that the final energy will be positive and the binary will be disrupted if MNS < (1/2)(M1 + M2). This condition occurs when the mass ejected from the system is greater than half of the initial total mass,

$$\Delta M > {1 \over 2}({M_1} + {M_2}),$$
(32)

where ΔM = M1MNS. If the binary is not disrupted, the new orbit becomes eccentric and expands to a new semi-major axis given by

$$a^{\prime} = a\left({{{{M_1} + {M_2} - \Delta M} \over {{M_1} + {M_2} - 2\Delta M}}} \right),$$
(33)

and orbital period

$$P^{\prime} = P{\left({{{a^{\prime}} \over a}} \right)^{3/2}}{\left({{{2a^{\prime} - a} \over {a^{\prime}}}} \right)^{1/2}}.$$
(34)

Note that we have not included any mention of the expected velocity kick that newly born neutron star or black hole will receive due to asymmetries in the supernova explosion. These kicks can be quite substantial, up to several hundred kilometers per second and, at least for observed pulsars, seem to be drawn from a Maxwellian distribution with a peak at 265 km s−1 [222]. In most cases, the kick will further increase the likelihood that the binary will become unbound, but occasionally the kick velocity will be favorably oriented and the binary will remain intact. If the kick is higher than the escape velocity of the cluster (typically less than 50 km s−1), it will also remove the remnant from the system. This mechanism may be very important in depleting the numbers of neutron stars and black holes in globular clusters.

We have seen that conservative mass transfer can result in a tighter binary if the more massive star is the donor. Non-conservative mass transfer can also drive the components of a binary together during a common envelope phase when mass and angular momentum are lost from the system. Direct mass loss through a supernova explosion can also alter the properties of a binary, but this process generally drives the system toward larger orbital separation and can disrupt the binary entirely. With this exception, the important result of all of these processes is the generation of tight binaries with at least one degenerate object.

The processes discussed so far apply to the generation of relativistic binaries anywhere. They occur whenever the orbital separation of a progenitor binary is sufficiently small to allow for mass transfer or common envelope evolution. Population distributions for relativistic binaries are derived from an initial mass function, a distribution in mass ratios, and a distribution in binary separations. These initial distributions are then fed into models for binary evolution, such as StarTrack [46] or SeBa [391, 352] in order to determine rates of production of relativistic binaries. The evolution of the binary is often determined by the application of some simple operational formulae, such as those described by Tout et al. [462] or Hurley et al. [230]. For example, Hils, Bender, and Webbink [220] estimated a population of close white dwarf binaries in the disk of the galaxy using a Salpeter mass function, a mass ratio distribution strongly peaked at 1, and a separation distribution that was flat in ln(a). Other estimates of relativistic binaries differ mostly by using different distributions [45, 242, 352, 350]. Because of the uncertainties in the form of the initial distributions, results from these simulations can differ by an order of magnitude or more.

4.3 Globular cluster processes

In the galactic field stellar densities are low enough that stars and binaries rarely encounter each other. In this environment binaries evolve in isolation with their properties and fate determined solely by their initial conditions and the processes described in Sections 4.1 and 4.2. In star clusters, stellar densities are much higher and close encounters between stars and binaries are common. Such encounters can affect a binary’s parameters and dramatically alter the evolution that would otherwise occur if the same binary were isolated. Some outcomes that are of particular interest for relativistic binaries include:

  • Reduction or increase of the period during distant encounters.

  • Exchange interactions where binary membership is changed in close encounters.

  • Binary formation during strong few-body interactions.

  • Encounter-induced binary mergers. Different types of interactions can produce these outcomes and we briefly describe the most important in the terms normally used by stellar dynamicists — the number of bodies and type of objects involved in the interaction. In all cases these interactions must be distinguished from the distant, weak interactions that drive two-body relaxation of clusters as described in Section 2.3. The interactions affecting binary parameters result from close encounters, the outcomes of which cannot be described statistically using the language of relaxation theory. Note also that the term binding energy is normally spoken of as though it were a positive quantity. Thus binaries with the largest (most negative) binding energy are the most tightly bound.

4.3.1 Single-single interactions

As the name suggests, single-single interactions are encounters between two individual stars. If the periastron of the encounter is sufficiently small, the stars may excite tidal oscillations in each other at the cost of some of their relative kinetic energies. If sufficient kinetic energy is dissipated, the stars can become bound and form a new binary.

The exact nature of the oscillations excited in the stars is not important, only that they dissipate sufficient kinetic energy. Furthermore, to obtain a bound orbit it is only necessary to dissipate sufficient kinetic energy that the total energy at apocenter is reduced to less than zero. The basic condition for tidal capture is:

$$\Delta {E_{{\rm{T1}}}} + \Delta {E_{{\rm{T2}}}} \geq {1 \over 2}\mu v_a^2$$
(35)

where ΔET1,T2 is the energy associated with the tidal oscillation in each star, μ is the reduced mass, and va is the relative velocity at apocenter (or infinity) [121]. Only a fraction ∼ (va/vp)2, where vp is the velocity at pericenter, of the energy must be converted to tidal oscillations in order for a capture to occur. For example, an encounter with a pericenter velocity of ∼ 100 km s−1 in a cluster with a velocity dispersion of ∼ 10 km s−1 only needs to dissipate ∼ 1% of the apocenter kinetic energy in order for a capture to occur. This mechanism favours the creation of binaries in encounters with large ratios between apocenter and pericenter velocity. Thus it tends to produce highly eccentric binaries with very small pericenter separations. These binaries will circularize over time due to continued excitation of tidal oscillations during the pericenter passages.

This process was once thought to be the dominant path for dynamically creating binaries in star clusters [121, 57] since two-body interactions are more likely than higher-order encounters [443]. Giant stars are the most likely to tidally capture a companion because their large radii give them a large cross section for interactions and their low densities make exciting tidal oscillations relatively easy. In particular binaries of neutron stars with degenerate companions were though to originate from the capture of the neutron star by a red giant [121]. It has been realized, however, that tidalcapture binaries may be rare because even giant stars must approach very closely, to within a few stellar radii, in order for tidal oscillations to dissipate enough energy for a capture to occur [121, 290, 393]. In these situations it is more likely that, depending upon the exact stellar equation of state, tidal capture events more commonly lead to a merger than a binary [407, 324, 24, 241, 405]. Because of the difficulty of exciting tidal oscillations in degenerate matter, forming a compact binary directly by tidal capture between two compact objects is highly unlikely.

4.3.2 Three-body interactions

During a close three-body encounter it is possible for one of the stars to gain kinetic energy and, if the relative kinetic energy of the remaining two stars is sufficiently reduced, they can form a bound pair. In general the least massive body will gain the highest velocity, since it is the easiest to accelerate, and the two more massive bodies will remain (this can also be seen as a consequence of energy equipartition). Heggie [196] showed that if Eb is the binding energy of the new binary, the rate of three-body binary formation is \( \propto E_b^{- 7/2}\). Thus three-body interactions tend to favour the creation of binaries with low binding energies. The binding energy can be increased by later encounters (see Section 4.3.3). For the case of an interaction where all stars have the same mass there is about a 10% probability that a binary formed by three-body interactions will gain sufficient binding energy to survive permanently [170, 443].

Three-body binary formation interactions are less common than tidal interactions because the probability of a close encounter between three stars is smaller than the probability of an encounter between two. Therefore stellar densities need to be higher for three-body binary formation to be efficient than they do for tidal interactions [170, 235, 443]. For the case of equal masses Spitzer [443] estimates a three-body binary formation rate of

$${\left({{{d{n_b}} \over {dt}}} \right)_{3{\rm{bdy}}}} = 1.91 \times {10^{- 13}}{\left({{{{\rho _n}} \over {{{10}^4}{\rm{p}}{{\rm{c}}^{- 3}}}}} \right)^3}{\left({{m \over {{M_ \odot}}}} \right)^5}{\left({{{10{\rm{km}}\,{{\rm{s}}^{- 1}}} \over {{v_m}}}} \right)^9}{\rm{p}}{{\rm{c}}^{- 3}}{\rm{y}}{{\rm{r}}^{- 1}},$$
(36)

where ρn is the number density of stars, m is their mass and vm is the three dimensional RMS random velocity. By contrast the rate of formation due to tidal capture in the same system (assuming a polytropic stellar equation of state with = 3) is

$${\left({{{d{n_b}} \over {dt}}} \right)_{{\rm{tidal}}}} = 1.52 \times {10^{- 8}}{\left({{{{\rho _n}} \over {{{10}^4}{\rm{p}}{{\rm{c}}^{- 3}}}}} \right)^2}{\left({{m \over {{M_ \odot}}}} \right)^{1.09}}{\left({{{{r_s}} \over {{R_ \odot}}}} \right)^{0.91}}{\left({{{10{\rm{km}}\,{{\rm{s}}^{- 1}}} \over {{v_m}}}} \right)^{1.18}}{\rm{p}}{{\rm{c}}^{- 3}}{\rm{y}}{{\rm{r}}^{- 1}},$$
(37)

rs is the radius of the stars. For m = 1 and vm = 10 km s−1 it is clear that the rate of binary formation by tidal capture is much higher at low ρn and that three-body binary formation only begins to dominate at ρn > 7.96 × 108 pc−3, a density higher than that found in most globular clusters. However, thanks to the strong scaling with mass, three-body binary formation becomes much more important for more massive stars. For the same vm but with m =10 M and assuming rsm0.8 then

$${\left({{{d{n_b}} \over {dt}}} \right)_{{\rm{3bdy}}}} = 1.91 \times {10^{- 8}}{\left({{{{\rho _n}} \over {{{10}^4}{\rm{p}}{{\rm{c}}^{- 3}}}}} \right)^3}$$
(38)

while

$${\left({{{d{n_b}} \over {dt}}} \right)_{{\rm{tidal}}}} = 1.0 \times {10^{- 6}}{\left({{{{\rho _n}} \over {{{10}^4}{\rm{p}}{{\rm{c}}^{- 3}}}}} \right)^2}$$
(39)

and three-body interactions begin to dominate at pn > 5.23 × 105 pc−3, much closer to the densities found in the cores of Galactic globular clusters. Because it is not necessary to excite tides during three-body interactions, pericenter passages need not be as close in order to form binaries so mergers are a less common outcome than they are for tidal interactions. Furthermore, because it is not necessary to excite tidal oscillations, it is possible for three-body interactions to form compact binaries directly from degenerate objects. It is worth noting that, in equal-mass systems, neither tidal nor three-body interactions are likely to produce many binaries over the lifetime of a star cluster [235, 170, 57, 443, 244]. For a system with a mass function the situation may be slightly different. Ivanova et al. showed that while the formation rate of binaries due to three-body interactions is negligable for solar-mass stars, for more massive stars the rate could be much higher [244, 245]. In particular, Ivanova et al. (2010) [245] showed that for very large mass ratios three-body binary formation could be important for forming certain types of compact binaries. Nevertheless, primordial binaries and the interactions involving them are critical for producing large numbers of relativistic binaries in globular clusters.

4.3.3 Binary-single interactions

Binary-single interactions, although formally still three-body interactions, differ because two of the stars are already bound. Several outcomes are possible in such an encounter depending on the relative kinetic energy of the star and binary and on the binding energy of the binary itself.

If the kinetic energy of the star is less than the binding energy of the binary, energy equipartition requires that the star accelerate and the binding energy of the binary increase (the period will become shorter). If, however, the kinetic energy of the star is greater than the binding energy of the binary, the star will donate kinetic energy to the binary and the binding energy will decrease (the period will become longer). For equal-mass systems this introduces a simple yet important classification for binaries. If |Eb| < kBT (Section 2.3) then the binary is said to be “soft” and, on average, will tend to have its period increased by encounters (softening). If |Eb| > kBT than the binary is said to be “hard” and will tend to have its period decreased by encounters (hardening). Assuming a Maxwellian velocity distribution, a hard binary gains on average 〈ΔEb〉 ∼ 0.2–0.4Eb per interaction [196, 239] (although this is only fully valid in the case of very hard binaries [235]). Since the encounter rate is proportional to the semi-major axis of the binary (∝ 1/Eb), the energy of a hard binary increases by ΔEb,rlx ∼ 0.6kbT/trlx per relaxation time [57]. This leads to the “Heggie-Hills Law” [196, 216] that states “hard binaries get harder while soft binaries get softer”. The end result of binary hardening can be a merger while the end result of binary softening can be the disruption of the binary.

In a multi-mass system the division between hard and soft binaries is not so clear since the relative energies depend upon the specific masses of all the bodies involved. For an individual encounter between a binary with member masses m1 and m2 and a single star of mass m3 travelling at velocity it is possible to define a critical velocity [237]:

$$v_c^2 = {{2({m_1} + {m_2} + {m_3})} \over {{m_3}({m_1} + {m_2})}}{E_b},$$
(40)

such that for vc < v the binary cannot capture the single star and the star can disrupt the binary (softening) and for vc > v then the binary can capture the single star and cannot be disrupted in the process (hardening). Thus the distinction between hard and soft binaries still exists in multi-mass systems. It has been shown by Hills that it is also possible to use the ratio of orbital speeds to define whether the binary will gain or lose energy [218].

The reference to capture alludes to another important process that can occur in binary-single interactions: exchange. In an exchange interaction one of the original binary members is replaced by the single star so the binary membership changes. As with other encounters, equipartition of energy favours the ejection of the lowest-mass object and exchange encounters are a way of introducing massive objects into binaries. It has been shown that in the limit of m3m1 or m2, the probability for a massive object to be exchanged into the binary is ∼ 1 [219]. Thus it is possible to create a compact binary from a binary where neither member was originally massive enough to become a compact object. This may be particularly effective for BHs since they are the most massive objects in evolved star clusters [436].

4.3.4 Binary-binary interactions

There are many possible outcomes for binary-binary interactions, especially in multi-mass systems, and they depend rather sensitively upon the initial conditions of the encounter. Therefore a quantitative theory of these interactions is lacking. For sufficiently distant encounters both binaries will “see” each other as single centers of mass and a hardening or softening interaction will occur in both binaries. There are also some general results for close encounters where all stars have equal mass. If one of the binaries is much harder than the other, it can be exchanged into the other binary as a single star. This produces a hierarchical three-body system, the future of which will be decided by further interactions [443]. One or both binaries can also be disrupted during the encounter. Numerical experiments show that at least one of the binaries is disrupted in ∼ 88% of cases [329, 330, 331, 144]. It is also possible for one or, more rarely, both binaries to exchange members with each other. Thus binary-binary interactions can provide all of the effects of binary-single interactions but can result in extensive binary destruction as well.

Higher-order interactions (those involving more than four bodies) are also possible but will be quite rare and are even less amenable to general quantitative analysis than are binary-binary encounters. Because binaries with small semi-major axes have a small interaction cross section, the majority of these encounters affect the binary while it is still quite wide by relativistic standards. Therefore it is unlikely to have a three-body (or higher-order) interaction while a binary is in the gravitational wave emission regime so globular clusters dynamics do not make it necessary to calculate merger waveforms perturbed by a third body [14].

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:

$${\ddot \vec r_i} = - \sum\limits_{i \neq j} {G{{({m_j}{{\vec r}_i} - {{\vec r}_j})} \over {\vert {{\vec r}_j} - {{\vec r}_i}{\vert ^3}}}}.$$
(41)

Here, mj is the mass of particle j and \({{\vec r}_{i,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. <monospace>Kira</monospace> integrates the N-body equations of motion using a 4th-order Hermite predictor-corrector scheme [317]. <monospace>Kira</monospace> 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. <monospace>Kira</monospace> 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 <monospace>Kira</monospace> 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 N3–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 1014 while a cluster may need to be advanced on timescales 1020 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,\,{f_m}(\vec x,\vec \upsilon, t)\), is 6N-dimensional in position-velocity phase space. The number of stars in the phase space volume d3 x d3 v is given by dN = fmd3 x d3 v. This description requires that either the phase space volume d3 x d3 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 d3 x d3 v 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. \(\phi (\vec x)\) is found by solving the Poisson equation, which can be written in terms of the distribution function as:

$${\nabla ^2}\phi = 4\pi G\sum\limits_i {\left[ {{m_i}\int {{f_{{m_i}}}(\vec x,\vec v,t){d^3}\vec v}} \right]}.$$
(42)

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 is calculated using a modified form of the Boltzmann Equation:

$${{\partial {f_m}} \over {\partial t}} + \vec v \cdot \nabla {f_m} - \nabla \phi \cdot {{\partial {f_m}} \over {\partial \vec v}} = \Gamma [f].$$
(43)

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 \((\vec x,\vec \upsilon \in \vec w)\) and takes on the form [57]:

$$\Gamma [f] = - \sum\limits_{i = 1}^6 {{\partial \over {\partial {w_i}}}[f(\vec w)\langle \Delta {w_i}\rangle ] + {1 \over 2}} \sum\limits_{i,j = 1}^{6} {{{{\partial ^2}} \over {\partial {w_i}\partial {w_j}}}[f(\vec w)\langle \Delta {w_i}\Delta {w_j}\rangle ]},$$
(44)

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 tcrtrlx 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, nonspherical 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 Stódoł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 N3–4 dependence of a direct N-body code. Thus a simulation of N ∼ 106 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 reproduce 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 highermass 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, rc the core radius, rt the tidal radius and trh the current estimated half-mass relaxation time.

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
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
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.

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

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.

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
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 ∼ 2 × 107 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
figure 15

Contours of constant number of X-ray binaries as a function of the “Verbunt parameters”. Here \(\Gamma = {\rho ^2}r_c^3/{\upsilon _c},\gamma = \rho/{\upsilon _c},\rho\), γ = ρ/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 2 × 10−6 < 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 TH. 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]

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–105 M. Ultimately IMBH is a catch-all term for BHs in the mass range between the stellar mass BHs (< 102 M) that are known to exist from stellar evolution and and the super-massive BHs (> 105 M) 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 ∼ 105 M 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.

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 ∼ 240 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.

Figure 16
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 103 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.

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].

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].

Figure 17
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.

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.

6 Prospects of Gravitational Radiation

In the coming decades, several gravitational wave observatories will begin detecting signals from relativistic binaries. It is likely that several of these binaries will have their origins within globular clusters. In this section, we review the detectors and their expected sources from globular clusters.

Ground-based interferometers such as the Laser Interferometer Gravitational-wave Observatory (LIGO) [299] in the United States, Virgo [476] in Europe, or the Large-scale Cryogenic Gravitational wave Telescope (LCGT) [287] in Japan, will be sensitive to gravitational radiation in the frequency range from a few Hz to a few kHz. In this frequency range, these detectors will be sensitive to the inspiral and merger of binaries containing neutron stars or black holes. These systems are somewhat rare in the Galactic globular cluster system, and the probability of a coalescence occurring within our Galaxy while the detectors are operating is vanishingly small. The enhanced design of LIGO was operational until the end of 2010, but did not detect any known gravitational wave signal. The enhanced design of Virgo continued operating into early 2012. Both Virgo and LIGO will be upgrading to advanced designs which will provide significant increases in sensitivity. Like advanced LIGO (aLIGO) and advanced Virgo (AdV), LCGT is a next-generation detector. It is unique in that it will be in an underground facility with and it will operate cryogenically. All of these detectors should be sensitive to double neutron star inspirals at distances of a few hundred Mpc, and double black hole mergers at distances of about a Gpc. At these distances, the number of extragalactic globular clusters is large enough that coalescences of neutron star or black hole binaries may be commonplace. Unfortunately, the angular resolution of the network of next-generation detectors will be much too coarse to identify any detected signal with a globular cluster. However, the gravitational waveform of the coalescence of these objects depends upon the orientation of the spins of the components relative to the orbital angular momentum. Although it is not entirely clear if field binaries will tend to have their spins aligned, there is no reason for this to be the case with binaries that have been formed dynamically through exchange of compact objects. Estimates of the event rates of compact object binary coalescence during the advanced detector era tend to weight the binary neutron star coalescences toward field binaries and the binary black hole coalescences toward cluster binaries. Successful parameter estimation for the properties of coalescing binaries in advanced detectors will require accurate modeling of the wave-form, including spin effects for both black holes and neutron stars, and tidal effects for neutron stars.

Binary systems containing white dwarfs will be brought into contact at orbital periods of a few hundred seconds. At contact, the system will either coalesce on a dynamical timescale if the mass transfer is unstable, or it will begin to spiral out if the mass transfer is stable. In any case, these systems will never produce gravitational radiation in the frequency band of ground-based interferometers. Space-based interferometers are capable of achieving sufficient sensitivity in the millihertz band. Thus, the sources for space-based interferometers will include white dwarf, neutron star, and black hole binaries. At these orbital periods, the systems are not very relativistic and the gravitational waveforms can safely be described by the quadrapole formula of Peters & Mathews [371, 370]. An angle averaged estimate of the typical strain amplitude is [55]

$${h_0} = 1.5 \times {10^{- 21}}{\left({{{{f_{{\rm{gw}}}}} \over {{{10}^{- 3}}{\rm{Hz}}}}} \right)^{2/3}}\left({{{1\,{\rm{kpc}}} \over r}} \right){\left({{\mathcal{M} \over {{M_ \odot}}}} \right)^{5/3}}.$$
(45)

At a typical globular cluster distance of ∼ 10 kpc and typical chirp mass of \({\mathcal M} \sim 0.5{M_ \odot}\), a relativistic WD-WD or WD-NS binary with Porb = 400 s will have a gravitational wave amplitude of 10−22. This value is within the range of proposed space-based gravitational wave observatories of the LISA class [55].

Many globular clusters lie off the plane of the galaxy and are relatively isolated systems with known positions. The angular resolution of LISA improves with signal strength. By focusing the search for gravitational radiation using known positions of suspected sources, it is possible to increase the signal-to-noise ratio for the detected signal. Thus, the angular resolution of LISA for globular cluster sources can be on the order of the angular size of the globular cluster itself at fgw > 1 mHz. Consequently, the orbital period distribution of a globular cluster’s population of relativistic binaries can be determined through observations in gravitational radiation. We will discuss the prospects for observing each class of relativistic binaries covered in this review.

WD-WD binaries that are formed from a common envelope phase will be briefly visible while the recently revealed hot core of the secondary cools. These objects are most likely the “nonflickerers” of Cool et al. [82] and Edmonds et al. [116]. WD-WD binaries formed through exchange interactions may very well harbor white dwarfs which are too cool to be observed. In either case, hardening through dynamical interactions will become less likely as the orbit shrinks and the effective cross section of the binary becomes too small. These objects will then be effectively invisible in electromagnetic radiation until they are brought into contact and RLOF can begin. During this invisible phase, the orbital period is ground down through the emission of gravitational radiation until the orbital period is a few hundred seconds [51]. With a frequency of 1 to 10 mHz, gravitational radiation from such a binary will be in the band of LISA [55].

WD-NS binaries that are expected to be progenitors of the millisecond pulsars must pass through a phase of gravitational radiation after the degenerate core of the donor star emerges from the common envelope phase and before the spin-up phase begins with the onset of mass transfer from the white dwarf to the neutron star. The orbital period at the onset of RLOF will be on the order of 1 to 2 minutes and the gravitational wave signal will be received at LISA with a signal-to-noise of 50–100 at a frequency of around 20 mHz for a globular cluster binary.

Binaries with significant eccentricity will have a spectrum of harmonics of the orbital frequency, with the relative strength of the nth harmonic for eccentricity e given by [371]

$$\begin{array}{*{20}c} {g(n,e) = {{{n^4}} \over {32}}{{\left\{{\left[ {{J_{n - 2}}(ne) - {J_{n - 1}}(ne) + {2 \over n}{J_n}(ne) + {J_{n + 1}}(ne) - {J_{n + 2}}(ne)} \right]} \right.}^2}} \\ {\quad \quad \quad \quad \quad \quad\,\, + \left. {(1 - {e^2}){{[{J_{n - 2}}(ne) - 2{J_n}(ne) + {J_{n + 2}}(ne)]}^2} + {4 \over {3{n^2}}}{{[{J_n}(ne)]}^2}} \right\}}, \\ \end{array}$$
(46)

where Jn is the Bessel function. The higher harmonics of sufficiently eccentric binaries (e > 0.7) can be detected by LISA even though the fundamental orbital frequency is well below its sensitivity band of 1–100 mHz [53]. Dynamical interactions may produce an eccentric population of 30–140 white dwarf binaries that would be present in the LISA data after a 5 year observation [487].

The Galactic globular cluster population of NS-NS binaries is expected to be quite small (≲ 10), they may have high eccentricities. As such, their higher harmonics may have sufficient power to be detectable in the LISA band even if their orbital period is too low. The binary pulsar B2127+11C is an example of a NS-NS binary in a globular cluster. In terms of the unknown angle of inclination i, the companion mass to the pulsar is M2 sin i ∼ 1 M and its eccentricity is e = 0.68 [296]. These binaries may also be detectable by LISA. If the globular cluster systems of other galaxies follow similar evolution as the Milky Way population, these binaries may be potential sources for LIGO as gravitational radiation grinds them down to coalescence. With an expected reach of over 300 Mpc, aLIGO can be expected to observe 40 field DNS systems per year [7]. It is uncertain how many DNS coalescences would be added to this rate by globular clusters and it would be difficult to distinguish between the two populations.

With their high eccentricities and large chirp mass, black hole binaries will also be good potential LISA sources for gravitational radiation from the galactic globular cluster system [52, 53]. Although the orbits are expected to have circularized by the time the enter the advanced LIGO frequency band, aLIGO is expected to observe DBH systems out to a typical distance of 1 Gpc, with optimally oriented systems detectable out to over 2 Gpc. Thus, it will be able to probe an enormous population of globular cluster systems. Estimates of the expected rate of DBH coalescences of field binaries range from a pessimistic 0.4 per year to an optimistic 1000 per year, with a realistic estimate of 20 per year [7]. Given the wide range of predicted enhancements in the DBH population in globular clusters, it is possible that this rate may be significantly enhanced by globular cluster binaries. As noted previously, if there is a preference for spin alignments in field binaries, the spin distribution of observed binaries may be used to distinguish globular cluster binaries from field binaries. Finally, we note that IMBHs may be potential sources of gravitational radiation for both space-based and ground-based detectors. The low mass end of the IMBH class (∼ 100 M) will be observable from ground-based detectors out to very large distances, while the extreme mass ratio inspiral of stellar mass compact objects around the high mass end of the IMBH class will be observable from space-based detectors.

The relatively close proximity of the galactic globular cluster system and the separations between individual globular clusters allows for the identification of gravitational radiation sources with their individual host clusters. Although the expected angular resolution of LISA is not small enough to allow for the identification of individual sources, knowledge of the positions of the clusters will allow for focused searches of the relativistic binary populations of the majority of the galactic globular clusters. Armed with a knowledge of the orbital periods of any detected binaries, concentrated searches in electromagnetic radiation can be successful in identifying relativistic binaries that may have otherwise been missed.

There may be a few double black hole systems within the Galactic globular cluster system that are either within the LISA band, or are sufficiently eccentric that their peak gravitational wave power is in a high harmonic that is within the LISA band. In either case, such systems would be easily detectable by a space-based gravitational wave detector due to their large chirp masses. However, a greater number of double black hole systems are expected to have been ejected from their birth globular clusters and these would now be inhabitants of the Galactic halo. If these systems have been formed dynamically, they are likely to be composed of two relatively massive black holes, and they may have chirp masses that are high enough for them to be detected at extragalactic distances.

7 Summary

Relativistic binaries are tracers of the rich dynamical evolution of globular clusters. The properties of this binary population are the result of an interplay between the gravitational dynamics of large N-body systems, the dynamics of mass transfer, the details of stellar evolution, and the effect of the gravitational field of the galaxy. The gravitational dynamics of globular clusters can enhance the population of short period binaries of main-sequence stars as well as inject compact objects such as white dwarfs and neutron stars into stellar binary systems. Once they are in such systems, the details of stellar evolution and mass transfer in close binary systems govern the likely end products of the dynamical interaction between the two stars. Furthermore, most models of the evolution of the core of a globular cluster rely on the gradual hardening and ejection of binary systems to delay the onset of core collapse. The hardening of binaries in the core of globular clusters will produce relativistic binaries, but it will also eventually eject these systems as they gain larger and larger recoil velocities in each subsequent encounter. The threshold for ejection from a globular cluster depends both upon the gravitational potential of the cluster itself and the gravitational potential of its environment generated by the Milky Way. As the globular cluster orbits the Milky Way, its local environment changes. Consequently, if other dynamical processes (such as gravothermal oscillations) do not dominate, the globular cluster’s population of relativistic binaries may also reflect the past orbital history of the globular cluster.

Over the last decade, observational techniques and technology have improved to the extent that significant discoveries are being made regularly. At this point, the bottleneck in observations of binary millisecond pulsars, low-mass X-ray binaries, and cataclysmic variables is time, not technology. As these observational techniques are brought to bear on more clusters, more discoveries are bound to be made. In the next decade, the possibility of using gravitational wave astronomy to detect relativistic binaries brings the exciting possibility of identifying the populations of electromagnetically invisible objects such as detached white dwarf and neutron star binaries and black hole binaries in globular clusters. These observations can only help to improve the understanding of the complex and interesting evolution of these objects and their host globular clusters.