Dynamically formed black hole+millisecond pulsar binaries in globular clusters

Abstract

1 INTRODUCTION

Radio pulsars in binary systems provide constraints on the processes that drive binary stellar evolution and unparalleled tests of general relativity in the strong-field regime. In most cases, the pulsars in these binaries are ‘recycled’. That is, the neutron star (NS) has been spun-up by accreting mass and angular momentum from its companion (Alpar et al. 1982). Compared to ‘normal’ pulsars, recycled pulsars exhibit greater stability and have much shorter spin periods (P_S ≲ 100 ms), both of which facilitate high precision measurements of the pulse arrival times (Lorimer 2008). Recycled pulsars with P_S ≲ 20 ms are referred to as millisecond pulsars (MSPs). The outcomes of binary evolution can be probed by using the recycled pulsar in such a system as a stable clock to precisely determine the binary's Keplerian orbital parameters and the properties of its component stars. If the recycled pulsar's companion is another neutron star, it is possible to measure post-Keplerian orbital parameters in a model-independent fashion and then compare these measurements with the predictions of various theories of gravity (Stairs 2004). The post-Keplerian parameters measured in the double pulsar binary PSR J0737−3039 offer the best test of gravity in the strong-field limit, to date, and are in excellent agreement with the predictions of general relativity (Kramer et al. 2006). New and better tests of general relativity may be possible by applying these techniques to a binary comprising a black hole (BH) and a recycled pulsar; however, such a system is yet to be discovered.

It is possible to produce a BH+recycled pulsar binary through standard evolutionary processes in an isolated, high-mass binary. The scenario requires that the primary (the initially more massive member of the binary, which evolves faster than its companion) produces an NS at the end of its lifetime and that the secondary produces a BH. This can occur if the primary transfers enough material to its companion to drive the companion's mass above the threshold for BH production. The NS created by the primary could then be recycled by accreting material from the companion before it evolves to become a BH. Under the assumption that these recycled pulsars would have lifetimes longer than 10¹⁰ yr, Narayan, Piran & Shemi (1991) placed an empirical upper limit on the formation rate of such BH+recycled pulsar binaries of 10⁻⁶ yr⁻¹ within the Galaxy. Population synthesis models by Lipunov et al. (1994) found a comparable formation rate, while similar studies by Sipior & Sigurdsson (2002), Voss & Tauris (2003), Sipior, Portegies Zwart & Nelemans (2004), and Pfahl, Podsiadlowski & Rappaport (2005) favoured lower BH+recycled pulsar binary formation rates of ∼10⁻⁷ yr⁻¹. Additionally, Pfahl et al. (2005) argued that the NSs in these systems would only be mildly recycled. Due to the rapid evolution of its massive companion, the NS, accreting at the Eddington limit, would only have time to accrete 10⁻³–10⁻⁴ M_⊙ of material before the secondary collapsed into a BH. These mildly recycled pulsars would only have lifetimes of 10⁸ yr. Even if the pulsars in these systems were completely recycled into long-lived MSPs, the population synthesis calculations suggest that most of these systems will undergo a gravitational wave-driven merger within ∼10⁸ yr. With formation rates of 10⁻⁷ yr⁻¹ and lifetimes of ≲ 10⁸ yr the number of BH+recycled pulsar binaries expected to exist in the Milky Way is only ∼10.

In a globular cluster (GC), a BH+recycled pulsar binary need not form directly from a primordial binary. The high stellar density in GCs leads to dynamical encounters between cluster members, which opens a wide array of evolutionary pathways that are inaccessible to isolated binaries. For example, a single NS in a GC can gain a companion by exchanging into a primordial binary during a three-body interaction. Subsequent evolution of these newly created binaries can result in the NS being spun-up into a recycled, MSP (Hills 1976; Sigurdsson & Phinney 1995; Ivanova et al. 2008). These encounters are evidenced by the enhanced formation rates of low mass X-ray binaries (LMXBs) and their progeny, MSPs, observed in GCs (e.g. Katz 1975; Verbunt & Hut 1987; Pooley et al. 2003; Camilo & Rasio 2005). Any of the BHs present in the cluster could acquire an MSP companion through similar interactions (Sigurdsson 2003). However, uncertainties in the size and nature of the BH population present in GCs complicate investigations of this formation channel.

Kulkarni, Hut & McMillan (1993) and Sigurdsson & Hernquist (1993) argued that the stellar mass BHs formed in a GC would rapidly sink to the centre of the cluster and eject one another in a phase of intense self-interaction. The frenzy of ejections results in a substantial depletion of the cluster's stellar mass BH population during the first Gyr of evolution. The fact that a firm BH candidate had not been identified in a GC during decades of observational study was in line with this theoretical picture. Given the meagre BH populations implied by these investigations, the dynamical formation of BH+MSP binaries in GCs has received little attention. After all, this channel closes if there is no population of BHs present in GCs. Nevertheless, the production of BH+MSP binaries through multibody interactions has been considered in dense stellar environments, analogous to GCs, that are likely to harbour BHs. Faucher-Giguère & Loeb (2011) showed that a few dynamically formed BH+MSP binaries should be present in the Galactic Center, where a cluster of ∼10⁴ stellar mass BH is expected to exist. This result indicates that BH+MSP binaries might be produced in GCs if the clusters retained some of their stellar mass BHs.

Recent observational efforts have shown that there are BHs present in some GCs, prompting a renewed interest in the nature of GC BH populations. A number of promising BH candidates have been discovered in X-ray observations of extragalactic GCs (Maccarone et al. 2007, 2011; Brassington et al. 2012, 2010; Irwin et al. 2010; Shih et al. 2010; Barnard, Garcia & Murray 2012; Roberts et al. 2012). Furthermore, three BH candidates have been identified in deep radio observations of Milky Way GCs; two candidates reside in M22 and one candidate is in M62 (Strader et al. 2012; Chomiuk et al. 2013). There is also a growing body of theoretical work suggesting that it may be possible for GCs to retain a substantial fraction of their stellar mass BH populations, under certain circumstances (Mackey et al. 2008; Morscher et al. 2013; Breen & Heggie 2013a,b; Sippel & Hurley 2013; Heggie & Giersz 2014). Motivated by these new results, we set out to explore how efficiently three-body exchanges produce BH+MSP binaries in GCs.

It has also been suggested that GCs may harbour intermediate-mass BHs (IMBHs; M ∼ 10²–10⁴ M_⊙). Previous studies have considered the consequences of interactions between MSPs and these IMBHs. The encounters could result in an MSP being significantly displaced from the GC core (Colpi, Mapelli & Possenti 2003), produce an IMBH+MSP binary (Devecchi et al. 2007), or populate the Milky Way halo with several high-velocity MSPs (Sesana et al. 2012). We will not include IMBHs in the models presented here, and instead focus on stellar mass BHs.

This paper is organized as follows. We describe the features of our models and motivate the range of initial conditions that our simulations explore in Section 2. In Section 3, we discuss the orbital parameters of the BH+MSP binaries produced in our models. We discuss the size of the BH+MSP binary population and the possibility of detecting such a binary in Section 4. Finally, in Section 5, we summarize and discuss our findings.

2 METHOD

Several characteristics of a GC can influence the nature of the BH+MSP binaries produced within it. The observed diversity in GC structure combined with uncertainties in the clusters’ binary and BH populations produce a vast range of plausible model inputs. Exploring this parameter space efficiently necessitates the use of fast, approximate methods, which themselves require a substantial amount of computation time. The study presented here will motivate and inform observations necessary for constraining the parameter space, and identify interesting regions in this space for follow-up with detailed N-body or Monte Carlo calculations.

We simulated the dynamical formation of BH+MSP binaries by evolving a variety of BH-binaries in fixed background GC models using the method described in Sigurdsson & Phinney (1995) and Clausen, Sigurdsson & Chernoff (2013, hereafter CSC13). Here, we briefly list the key physical effects included in our models:

• Fokker–Planck advection and diffusion of single BHs or binaries containing at least one BH in static, multimass background clusters;

• direct integration of binary–single encounters to track changes in the orbital parameters of the binary, exchanges of members of the binary, destruction of the binary, and physical collisions between stars and/or BHs;

• changes in the semimajor axis and eccentricity of the binary on account of gravitational wave emission;

• changes in the masses of binary components following mergers and recoil of BHs after gravitational wave-driven mergers.

We will describe our implementation of these effects in more detail below, but for a complete description we refer the reader to Sigurdsson & Phinney (1995) and CSC13.

We focus on multimass King models whose total stellar content derives from a simple, evolved mass function. We start once the static, equilibrium approximation implicit in the King model treatment is reasonably accurate. The calculation begins only after stars of several solar masses have evolved off the main sequence because prior evolutionary phases involved significant mass-loss. Collisional truncation of the wide binaries has already occurred. The cluster's BHs and NSs of interest in this paper were born much earlier. Many important and complicated processes have run to completion, e.g. dynamical friction segregated heavy objects and promoted formation of a self-gravitating core. Interactions involving binaries have ejected all but a few of the heaviest remnants. We do not model any of these processes. We assume a priori (and without the guidance of direct observations) (1) the single versus binary content and (2) the number and mass of resident BHs. We focus on evolutionary stages when cluster changes occur on the characteristic half-mass relaxation time-scale. We specify structural properties that match observed GCs, apparently in the late time, stable phase seen in fully dynamical models of GC evolution (e.g. Hurley, Aarseth & Shara 2007; Chatterjee et al. 2010). These assumptions capture the outcomes of all the early, complicated phases of evolution and provide the initial conditions for our calculations. We generated four background cluster families to study the impact of GC structure on the formation rate and properties of the BH+MSP binaries. The names, central number densities n_c, mean central velocity dispersions |$\bar{v}_{\rm m}$|⁠, concentrations c, and total masses M_GC of the these cluster families are listed in Table 1.

Each cluster's stellar population was assumed to have an initial mass function of the form ξ(m) ∝ m^−α. We used a broken power law similar to the one given in Kroupa (2001) with α values of 1.3 and 2.35 for stars with mass m < 0.35 M_⊙ and m > 0.35 M_⊙, respectively. We chose a main-sequence turn-off mass of 0.85 M_⊙ and assumed that stars above the turn-off mass had evolved completely to white dwarfs (WDs), NSs, or BHs. The evolved stellar population was then binned into 10 mass groups. Main-sequence stars and WDs with masses in the range 0.08–1.2 M_⊙ were grouped into eight bins with widths of ∼0.15 M_⊙. A ninth mass bin for stars with mass 1.2 M_⊙ < m < 1.5 M_⊙ contained the NSs. The final, high-mass bin held the BHs. WD masses were computed using the semi-empirical initial–final mass relation derived by Catalán et al. (2008). We used an NS mass of 1.4 M_⊙ and considered NS retention fractions, f_ret, of 5, 10, and 20 per cent (Sigurdsson & Phinney 1995; Pfahl, Rappaport & Podsiadlowski 2002; Ivanova et al. 2008).

The BH populations in GCs are poorly constrained, so we considered a range of BH masses, M_BH, and population sizes. By analogy to the BHs found in the Galaxy, we used M_BH = 7 M_⊙ (Özel et al. 2010; Farr et al. 2011). A second value of M_BH = 15 M_⊙ was motivated by Strader et al. (2012), who used mass segregation arguments to estimate the masses of the BH candidates in M22 to be ∼15 M_⊙. We also used M_BH = 35 M_⊙ in some of our simulations. Observations of the stellar mass BH candidates in extragalactic GC indicate that these BHs may have masses ≳ 30 M_⊙ (Maccarone et al. 2007; Irwin et al. 2010; Barnard et al. 2012; Clausen et al. 2012). We controlled the number of BHs present in each simulation by truncating the high-mass end of the initial mass function. In a subset of our simulations, we assumed that nearly every BH formed in the cluster was ejected, leaving a lone BH.¹ In other simulations, we allowed the clusters to retain several BHs, with the number of BHs (N_BH) in the range 7-191. In such simulations, the additional BHs were included as part of the background cluster and comprised the high-mass group of our multimass King model. As such, all of these BHs were assumed to be single and held in mass-segregated equilibrium with the rest of the cluster. Including these BH populations in our static background clusters had a substantial impact on the mass and structure of these clusters. To facilitate comparisons between models with different values of N_BH, we adjusted the free parameters in the cluster models, slightly, to ensure that n_c, |$\bar{v}_{\rm m}$|⁠, c, and M_GC changed by ≲ 20 per cent when N_BH was varied. The ranges of these parameters are given in Table 1.

We generated 2000 BH-binaries and followed their evolution, one-by-one, in each of the background clusters. We refer to this group of calculations as one simulation. The initial configuration of the BH-binary used in each run was determined using the following distributions. The BH's companion was randomly selected from the cluster's evolved mass distribution. We selected the initial eccentricity, e, from a thermal distribution, f(e) ∝ 2e. The distribution of initial semimajor axes was assumed to be flat in log a between 10⁻³ au and a_max. The initial BH-binaries in our simulations are the products of complicated stellar and dynamical evolution, the latter of which occurs at different rates in different clusters. We set a_max to 100, 33, 15, and 10 au for simulations in cluster A, B, C, and D, respectively, in an attempt to capture the broad range of the possible outcomes of this evolution (CSC13). Finally, the initial position and velocity of each BH-binary was selected from the radial density and velocity distributions of the third most massive mass group (CSC13).

We sought to determine how often the BHs in these binaries acquired MSP companions. Each binary was evolved in the cluster potential with dynamical friction and Fokker–Planck diffusion calculated explicitly. Each run continued until either the BH was ejected from the cluster or the run had covered a maximum time of t_max = 10¹⁰ yr. The probability that the binary would undergo a close encounter with a background star was calculated continuously (see Sigurdsson & Phinney 1995 for a detailed description). In simulations of clusters that contained more than one BH, this included the possibility that the binary would interact with a background BH. The probability of such an encounter was effectively 0 in the clusters that only contained one BH. If an encounter was deemed to have occurred, we integrated the three-body interaction between the binary and the background star. These encounters can alter the orbital parameters and the components of the binary, due to a collision/merger or an exchange. The interaction may also result in the disruption of the binary. If the binary survived the encounter, we continued to follow its evolution within the cluster.

Although we were primarily concerned with interactions between BH-binaries and single background objects (stars and remnants), some situations required us to model encounters between single BHs and background binaries. If the BH that we were following lost its companion as the result of a merger or a disruptive encounter, we used the method described in CSC13 to determine whether or not the BH would be able to exchange back into a binary. This physical process depends upon some assumptions about the nature and size of the binary populations present in GCs. We considered three binary populations with different compositions. One population, labelled OBS, was motivated by observational studies of photometric binaries in GCs, and was composed primarily of binaries that contained two main-sequence stars (e.g. Milone et al. 2012). The second population, labelled FIR, was motivated by the theoretical study presented in Fregeau, Ivanova & Rasio (2009). In addition to showing that the hard binary fraction in the core of a cluster will increase over time, this work found that most binaries in the core will harbour at least one non-luminous member, making it difficult for observational studies to identify these systems. The final binary population was optimized to interact with the single BHs and is labelled the OPT population. None of the background binary populations included BHs. We varied the size of the binary population by adjusting the binary fraction, f_b, within observational constraints.

A single BH could also be incorporated into a binary that is produced as the result of an encounter between the BH and two single stars. The cross-section for such an encounter is typically very small, but it scales as roughly |$M_{\rm BH}^5$|⁠. We used the time-scale for three-body binary formation given in Ivanova et al. (2005) to assess whether this processes is likely to occur before a single BH would exchange into an existing binary, as described above. In most cases, the exchange time-scale of a few Gyr is orders of magnitude shorter than the three-body binary formation time-scale (CSC13). However, in simulations of cluster D that used a 35 M_⊙ BH, or several lower mass BHs, the two time-scales are within a factor of ∼2 of each other. Although three-body binary formation becomes competitive with exchanges in these extreme scenarios, we did not include this process in our models because it does not significantly reduce the amount of time it takes a BH to regain a binary companion.

3 BH+MSP BINARY ORBITAL PARAMETERS

The results of our simulations are summarized in Table 2. The first column lists an identification number for each simulation. The next six columns describe the initial conditions used for the simulations, noting the background cluster, binary population, binary fraction (f_b), NS retention fraction (f_ret), the BH mass (M_BH), and the number of BHs (N_BH). The final six columns list the number of BH+NS binaries produced in that simulation (N_BH+NS), the number of BH+BH binaries produced in the simulation (N_BH+BH), the median orbital period (P_B) of the BH+NS binaries, the standard deviation of the BH+NS binary orbital period distribution (s), the fraction of time that a BH+NS binary exists (τ_BH+NS, see Section 4), and the number of BH+NS binaries expected to be present in the cluster at any given time (p_BH+NS, see Section 4). We note that our code did not track whether or not an NS had been recycled into an MSP. Thus, in the discussion that follows we will examine the nature of all of the BH+NS binaries produced in our simulations. We will consider the fraction of NSs that are MSPs in Section 4.

Encounters between BH-binaries that do not contain an NS and background main-sequence stars, giants, or WDs can influence the nature of the BH+NS binary population. These interactions occur more frequently because these objects are far more common than NSs. Such encounters can limit the size of the BH+NS binary population by removing BHs from the cluster as a result of a superelastic binary–single collision. However, encounters of this strength were rare in our simulations. In most cases, more that 75 per cent of the 2000 BHs that we followed during a simulation were retained by the cluster for the entire 10¹⁰ yr run. In the rare cases where most of the BHs were ejected, the majority of BHs were ejected after merging with an NS and receiving a recoil that exceeded the cluster's escape velocity (Shibata et al. 2009). While encounters with non-NSs did not efficiently remove BHs from the cluster, these interactions did impact the size of the BH+NS star population by destroying BH-binaries. An interaction can either disrupt the binary into three single stars or induce a physical collision and merger between the BH and its companion. Either way, at the end of the encounter the BH has lost its companion and must begin the slow process of exchanging back into a binary, delaying the formation of a BH+NS binary. We defer the discussion of these encounters and BH-binaries, which are potential X-ray sources, to a future publication. Below we will focus on the properties of the BH+NS binary population and encounters involving these binaries.

In our simulations, there were two processes that drove the evolution of a BH+NS binary's orbital parameters. Encounters between the binary and background stars changed the semimajor axis (a) and eccentricity (e) of the binary impulsively. In most cases, an encounter resulted in the binary becoming more tightly bound or ‘hardened’. The emission of gravitational radiation also modified the BH+NS binaries’ orbital parameters, driving a and e towards zero. Since the orbital parameters of the binaries were constantly changing, we resampled the simulation output in even time intervals to ensure that each orbital configuration that the BH+NS binaries evolved through was properly weighted. We chose a resampling time step of 10⁷ yr. We checked that our choice of 10⁷ yr intervals did not bias the resampled orbital parameter distributions by repeating the analysis of a subset of our simulations with finer (10⁵ yr) time resolution. A coarser resampling, with steps of 10⁸ yr, failed to capture the wings of the semimajor axis distribution, where rapid evolution occurs. Examples of the resampled data from three representative simulations are shown in Fig. 1. The figure shows the joint distribution for the eccentricity and semimajor axis of the BH+NS binaries produced in simulations using background clusters B, C, and D.

Figure 1.

Joint e-a distribution for BH+NS binaries produced in simulations 42 (left-hand panel), 26 (centre panel), and 10 (right-hand panel). The eccentricity distribution is similar in all of these simulations. The semimajor axis distribution shifts to larger a as the cluster central density decreases from n_c = 10⁶ in the left-hand panel to n_c = 10⁵ in the right-hand panel.

Open in new tabDownload slide

3.1 Eccentricity distribution

In nearly all of our simulations, the eccentricity distribution of the BH+NS binaries was roughly thermal for e ≲ 0.9. At higher eccentricity the distribution flattened out and turned over. That is to say, there were fewer binaries with e ≳ 0.9 than expected given the f(e) = 2e distribution. This is because the emission of gravitational radiation acted to quickly circularize such high-eccentricity systems. Given these results, we expect that the mean e of dynamically formed BH+NS binaries in GCs will be in the range 0.6–0.7. In the skewed eccentricity distributions seen in our simulations, the median value of e exceeds the mean so most binaries will have eccentricities larger than these average values. The only simulations that did not result in a thermal eccentricity distribution were those run in background cluster A. As can be seen in Table 2, very few BH+NS binaries were produced in this low-density cluster. We were, therefore, unable to study the eccentricity distributions in these poorly sampled cases.

3.2 Orbital separation distribution

Of all of the parameters varied in our study, the background cluster's structural properties had the strongest impact on the semimajor axes of the BH+NS binaries. For this reason, we will describe the distributions of semimajor axes and orbital periods cluster-by-cluster. Once we have described the results for each cluster, we will investigate the origin of the observed trends.

GC A was the least massive cluster considered in our study. Furthermore, because of its low density, this cluster also had the lowest encounter rate. NSs only exchanged into 45 of the 1.8 × 10⁴ BH-binaries that we evolved in this cluster. The corresponding specific BH+NS binary formation rate of ∼3 × 10⁻¹⁰ yr⁻¹ M_⊙⁻¹ is smaller than the specific BH+recycled pulsar formation rate predicted for the field. Although the formation rate of BH+MSP binaries is not enhanced by dynamical processes in this low-density cluster, given the short lifetimes of BH+recycled pulsar binaries formed through standard evolutionary channels (see Introduction), it is likely that any BH+MSP binaries present in such a cluster today were formed through dynamical interactions. Thus, the properties of these clusters’ present-day BH+NS binary populations should be described by the systems formed in our simulations. The semimajor axes of these rare BH+NS binaries spanned a range of nearly two orders of magnitude, from 1.2 to 113 au. The BH+NS binaries produced in this cluster were extremely wide because neither of the hardening mechanisms described above were effective. The binaries underwent very few encounters because of the cluster's low density, and the effects of gravitational radiation were negligible at such large orbital separations. Accordingly, we predict that any BH+MSP binaries present in low-density clusters are likely to have orbital periods of several decades.

The dynamically formed BH+NS binaries in globular cluster B had significantly smaller orbital separations. In this higher density cluster, frequent encounters hardened some of the BH+NS binaries to small enough a that gravitational radiation effects became important. Fig. 2 shows the cumulative distributions of a and P_B amongst present-day BH+NS binaries for all of our simulations in cluster B (i.e. simulations 10–25). We constructed a present-day population by only selecting binaries that exist at t > 8 × 10⁹ yr. The median values of a in these distributions fell in the range 0.42–1.6 au. It is evident from Fig. 2 that the orbital separation is influenced by the value of M_BH. Simulations with higher mass BHs produced tighter BH+NS binaries. This trend is amplified in the orbital period distributions because the semimajor axis and orbital period are linked by an additional factor of |${\sim } M_{\rm BH}^{-1/2}$|⁠, specifically |$P_{\rm B} = \scriptstyle\sqrt{a^{3}/(M_{\rm BH} + M_{\rm NS})}$|⁠. At a given value of M_BH, the orbital period distributions are fairly similar to one another, despite significantly different assumptions about the number of BHs in the clusters and the clusters’ binary populations. This suggests that these properties do not strongly impact the orbital periods of the BH+NS binaries produced in GCs with structure similar to that of cluster B. The median BH+NS binary orbital periods in all of the simulations in cluster B fell between 16 and 260 d. The standard deviation of the P_B distribution within a particular simulation was much wider (see Table 2).

Figure 2.

Cumulative distributions of the semimajor axis (left) and orbital period (right) for a subset of the BH+NS binaries produced in each of the simulations in cluster B. Each curve shows the distribution of a or P_B in a particular simulation. The distributions shown only include systems that exist at t > 8 × 10⁹ yr, so they correspond to the present day. The colour of each curve denotes the mass of the BH(s) used in a particular simulation, with 35, 15, and 7 M_⊙ BHs denoted by red, green, and blue, respectively. Note the clear relationship between M_BH and orbital separation: simulations with more massive BHs produced BH+NS binaries with smaller semimajor axes. This is the opposite of the trend seen in higher density clusters.

Open in new tabDownload slide

For simulations in cluster C, the BH+NS binary orbital separations were smaller still. The cumulative distributions of the semimajor axis and orbital periods for the BH+NS binaries formed in this cluster are shown in Fig. 3. The median values of a in these simulations ranged between 0.15 and 0.43 au. Simulation 32 resulted in the largest median a, and from Fig. 3 it is clear that the cumulative a and P_B distributions for simulation 32 deviate from distributions seen in other simulations. The critical difference between this simulation and the others run in cluster C was the large number of BHs (N_BH = 19) retained by the cluster. As we will discuss in Section 4, the evolution of the BH+NS binaries can be significantly altered when there are several relatively massive BHs present in the cluster. If we exclude the simulations that used more than one 15 M_⊙ BH (i.e. simulations 31, 32, and 38) from our analysis, the range of median semimajor axes reduces to 0.15–0.20 au. This range is much narrower than that observed in our simulations in cluster B. The semimajor axes of the BH+NS binaries produced in cluster C only depended weakly on many of the input parameters, including M_BH and f_b. In addition to being a much smaller effect, the relationship between BH mass and semi-axis seen in cluster C is the reverse of what was seen in cluster B. Here, the simulations with the lowest mass BHs produced the BH+NS binaries with the smallest semimajor axes. Most of the BH+MSP binaries in GCs with structures similar to cluster C will have orbital periods shorter than 10 d.

Figure 3.

Cumulative distributions of the semimajor axis (left) and orbital period (right) for a subset of the BH+NS binaries produced in each of the simulations in cluster C. Each curve shows the distribution of a or P_B in a particular simulation. The distributions shown only include systems that exist at t > 8 × 10⁹ yr, so they correspond to the present day. The colour of each curve denotes the mass of the BH(s) used in a particular simulation, with 35, 15, and 7 M_⊙ BHs denoted by red, green, and blue, respectively. Note the clear relationship between M_BH and orbital separation: simulations with more massive BHs produce BH+NS binaries with larger semimajor axes. This is contrary to the trend seen in the lower density cluster B. Note, however, that simulations with M_BH = 15 M_⊙andN_BH > 1 do not follow this trend.

Open in new tabDownload slide

Finally, the shortest period BH+NS binaries observed in our simulations were produced in cluster D. The encounter rate was highest in cluster D, so BH-binaries were rapidly hardened to small orbital separations in the simulations that used this background cluster. Fig. 4 shows the present-day cumulative distributions of a and P_B for each of the simulations performed in cluster D (i.e. simulations 42–51). As was seen in the simulations that used cluster C, the simulation with M_BH = 7 M_⊙ produced BH+NS binaries with the smallest semimajor axes. However, the median value of a only varied slightly from simulation to simulation in cluster D. Simulation 42 had the smallest median a = 0.11 au. The median value of a was largest in simulations 47 and 48, both of which had median a = 0.17 au. We expect that many BH+MSP binaries present in densest GCs will have P_B ≲ 5 d.

Figure 4.

Cumulative distribution of the semimajor axis (left) and orbital period (right) for a subset of the BH+NS binaries produced in each of the simulations in cluster family D. Each curve shows the distribution of a or P_B in a particular simulation. The distributions shown only include systems that exist at t > 8 × 10⁹ yr so they correspond to the present day. Note that simulations with 7 M_⊙ BHs produce binaries with smaller semimajor axes than the simulations that used 35 M_⊙ BHs. This is contrary to the trend seen in the lower density cluster B.

Open in new tabDownload slide

To review, we have shown that the BH+NS binaries produced in high-density clusters have smaller semimajor axes than those produced in low-density clusters. Furthermore, in the high-density GCs, the binaries with lower mass BHs had smaller orbital separations. However, the size of a BH+NS binary's semimajor axis only depended weakly on the mass of the BH. When we changed M_BH by a factor of 5, the median semimajor axis only changed by a factor of ∼1.5. The opposite trend was seen in lower density clusters (e.g. cluster B). In such clusters, binaries with higher mass BHs tended to have smaller semimajor axes than those with low-mass BHs. Additionally, the mass of the BHs had a stronger impact on the orbital separations of the binaries produced in these simulations (compare Figs 2 and 4). It is clear that the BH+NS binaries present in high-density clusters (similar to clusters C and D) are in a different evolutionary phase than those in lower density clusters (similar to cluster B).

Now let us compare these rates. For wide binaries, |〈da/dt〉_enc| ≫ |〈da/dt〉_GW| and encounters with single stars will be the dominant hardening mechanism. On the other hand, the evolution of binaries with small orbital separations will occur more rapidly through the emission of gravitational radiation than through encounters with background stars. The initial semimajor axes of the dynamically formed BH+NS binaries considered here will be relatively large, as is typical of systems produced in exchange interactions. Thus, the semimajor axis of a newly formed BH+NS binary will shrink fairly rapidly as a result of encounters with single stars. As encounters reduce the binary's orbital separation, |〈da/dt〉_enc| also declines. Eventually, a (and |〈da/dt〉_enc|) will be such that encounters do not efficiently harden the binary, and the orbital evolution will effectively stall. The binary will spend a large fraction of its lifetime at roughly constant semimajor axis, awaiting a final encounter that will reduce its semimajor axis to the point at which gravitational radiation takes over. At this point, the binary's orbital evolution will occur more rapidly as a continues to shrink. Therefore, the semimajor axis distributions will be dominated by BH+NS binaries with orbital configurations that fall in the transition between the encounter and gravitational radiation dominated regimes.

We explored the transition between the encounter dominated phase and the gravitational radiation-dominated phase by examining the ratio 〈da/dt〉_GW/〈da/dt〉_enc. Using equations (2) and (3), we calculated the ratio for each point in the BH+NS binaries’ resampled evolutionary tracks. The cumulative distributions of 〈da/dt〉_GW/〈da/dt〉_enc from a subset of our simulations are shown in Fig. 5. Here, we have only plotted simulations in clusters C and D. As discussed above, the semimajor axis distributions in these clusters exhibit similar trends. The left-hand panel shows how 〈da/dt〉_GW/〈da/dt〉_enc is distributed at t < 1 × 10⁹ yr. The centre panel shows how these distributions have changed after ∼5 Gyr of evolution. Note that the distribution of 〈da/dt〉_GW/〈da/dt〉_enc for simulations of cluster D that used 35 M_⊙ BHs had hardly changed. The distributions in all of the other simulations, with lower encounter rates, evolved towards these stationary curves. The right-hand panel shows how 〈da/dt〉_GW/〈da/dt〉_enc is distributed at the present day, after an additional ∼3 Gyr of evolution. Most of the 〈da/dt〉_GW/〈da/dt〉_enc distributions from simulations in cluster D (solid lines) lie on top of one another. Simulations in cluster C (dashed lines) that used M_BH = 35 M_⊙ also closely follow this trend, while some of the simulations with lower mass BHs exhibit smaller values of 〈da/dt〉_GW/〈da/dt〉_enc.

Figure 5.

Cumulative distributions of the ratio of the rate of semimajor axis evolution due to the emission of gravitational radiation to rate of semimajor axis evolution due to encounters with background stars (given in equations 3 and 2, respectively). Each curve shows the distribution of 〈da/dt〉_GW/〈da/dt〉_enc in a particular simulation in cluster C (dashed lines) or cluster D (solid lines). The three panels illustrate how the distribution evolves over time. The distributions shown are for t < 1 × 10⁹ yr (left-hand panel), 4 × 10⁹ < t < 6 × 10⁹ yr (centre panel), and t > 8 × 10⁹ yr (right-hand panel). The colour of each curve denotes the mass of the BH(s) used in a particular simulation, with 35, 15, and 7 M_⊙ BHs denoted by red, green, and blue, respectively. The binaries evolve towards a steady distribution in 〈da/dt〉_GW/〈da/dt〉_enc. Simulations with the largest encounter rates (those in cluster D and those in cluster C with M_BH = 35 M_⊙) have reached the steady configuration by the current epoch. The simulations with lower encounter rates are still approaching this configuration, from the left, but have not reached this state after 10¹⁰ yr of evolution.

Open in new tabDownload slide

We interpret this evolution as follows. Given sufficient time, the distribution of 〈da/dt〉_GW/〈da/dt〉_enc in any cluster will approach the steady configuration seen in the right-hand panel of Fig. 5. At this stage, the orbital evolution of most of the BH+NS binaries has slowed as they make the transition from the encounter dominated phase to the gravitational radiation-dominated phase. The median value of 〈da/dt〉_GW/〈da/dt〉_enc in this steady configuration is ∼2 × 10⁻³. Even though |〈da/dt〉_GW| is still a few orders of magnitude smaller than |〈da/dt〉_enc| at this point, this ratio corresponds to the beginning of transition between the encounter and gravitational wave-dominated phases. Given the steep dependence of 〈da/dt〉_GW on a, such a binary is only about one encounter away from an orbit in which 〈da/dt〉_GW ∼ 〈da/dt〉_enc. We note that despite the similarities illustrated in Fig. 5, the distributions from each simulation are formally distinct from one another, in most cases. We performed the two-sample Kolmogorov–Smirnov test between each pair of distributions and found that only those from simulations 44, 45, and 46 were consistent with coming from the same parent distribution. However, our interpretation does not require that the distributions be identical, it only requires that each distribution is dominated by systems that are in the transition phase.

Notably, many of the simulations that used 15 M_⊙ BHs do not seem to follow the evolutionary path sketched above. Fig. 5 shows that the BH+NS binaries produced in these simulations do not settle into the steady 〈da/dt〉_GW/〈da/dt〉_enc distribution seen in simulations that used both higher and lower mass BHs. Furthermore, the semimajor axes of the BH+NS binaries produced in most of the simulations that used 15 M_⊙ BHs do not appear to follow the monotonic relationship with M_BH implied by equation (4). If the BH+NS binaries did follow this relationship, one would expect the curves corresponding to simulations with 15 M_⊙ BHs (shown in green) to fall between the curves for simulations that used 7 M_⊙ (shown in blue) and 35 M_⊙ (shown in red) BHs in the left-hand panels of Figs 3 and 4. The only simulation that obeys this relation is simulation 30, which considered a cluster that only contained one 15 M_⊙ BH. This suggests that the larger BH populations used in all of the other simulations with 15 M_⊙ BHs are likely to blame for the modified orbital evolution of the BH+NS in these models. In fact, as we will discuss in Section 4, many of the BH+NS binaries produced in background clusters that contained a large population of fairly massive BHs were not hardened, but rather quickly destroyed when a second BH exchanged into the binary and ejected the NS. This accounts for the lack of BH+NS binaries with small a in the simulations that used 15 M_⊙ BHs.

In clusters with low encounter rates, many of the BH+NS binaries will not be hardened fast enough to reach the transition from the encounter-dominated regime to the gravitational radiation-dominated regime within a Hubble time. Therefore, the orbital parameters of most BH+NS binaries in these clusters are determined by binary–single encounters alone. Accordingly, we expect that the semimajor axes of the binaries in these clusters will follow a different scaling relationship than those in the high encounter rate clusters explored above. The rate at which encounters harden binaries is proportional to the mass of the BH (see equation 2), so we would expect for the simulations that used higher mass BHs to produce BH+NS binaries with smaller semimajor axes in the encounter dominated regime. This is what we observed in simulations in cluster B; the binaries with higher mass BHs were hardened to smaller a than the binaries with low-mass BHs.

4 BH+MSP BINARY POPULATION SIZE

The 51 simulations presented here investigated how several parameters impact the dynamical formation of BH+NS binaries in GCs. Below, we will describe how each of these traits affects the likelihood that a BH+NS binary exists within a cluster. We will use two metrics to characterize the likelihood that a cluster harbours a BH+NS binary. The first metric involves the average time fraction that a BH in our simulations had an NS companion, τ_BH+NS = ∑_it_BH+NS,i/(N_runst_max). Here, t_BH+NS,i is the amount of time that a BH had an NS companion during a particular run, the sum is over all of the runs in a single simulation, t_max = 10¹⁰ yr is the duration of each run, and N_runs is the number of runs in the simulation. In practice, N_runs is slightly less than 2000 because we rejected a small number of runs in which the initial BH-binary was randomly selected to contain an NS. The second metric that we use is the number of BH+NS binaries a cluster is expected to contain at any given time p_BH+NS = N_BHτ_BH+NS. In simulations with multiple BHs, τ_BH+NS corresponds to the probability that any one of these BHs has an NS companion at a given time, and must therefore be scaled by N_BH to find the number of BH+NS binaries expected to exist in the cluster. The values of these metrics for each of our simulations are listed in Table 2.

As was the case with the orbital parameter distributions, the structure of the background GC had a large impact on the values of τ_BH+NS. In cluster A (n_c = 10⁴ pc⁻³), the BH+NS binaries were long lived, with a mean lifetime of 7 × 10⁹ yr. However, because of the low encounter rate, these binaries were produced so rarely that τ_BH+NS ≲ 2 × 10^{− 3} in every simulation in this cluster. Simulations in cluster D (n_c = 10⁶ pc⁻³) had the opposite problem. BH+NS binaries were formed efficiently in this cluster, but they were also rapidly driven to coalescence because of the high encounter rate. The intermediate encounter rates in clusters B (n_c = 10⁵ pc⁻³) and C (n_c = 5 × 10⁵ pc⁻³) struck a balance between the production and destruction of BH+NS binaries. All other parameters held equal, simulations in clusters B or C had the largest values of τ_BH+NS. The BHs in these clusters typically spent a few percent of their lifetimes with an NS companion.

The amount of time that the BHs in our simulations had NS companions decreased with increasing M_BH. Both of the evolutionary processes that drive the BH+NS binaries to merge speed up when the mass of the BH is increased. A more massive BH leads to a larger gravitational focusing cross-section, which increases the encounter rate. Furthermore, orbital decay is more rapid through gravitational wave emission when M_BH increases.

Our models captured how the size and composition of a cluster's binary population will influence the production of BH+NS binaries. The size of the background binary population was changed by varying the binary fraction. Three different binary populations were considered, an observationally motivated population (OBS), a theoretically motivated population (FIR), and a population optimized to interact with the BHs (OPT; see Section 2 and CSC13 for a more complete description of these binary populations). The assumed background binary population did not influence τ_BH+NS in clusters A and B. In these clusters, mergers and disruptive encounters rarely occurred. Accordingly, there were few single BHs in these clusters that needed to interact with the background binary population to acquire a new companion. In clusters C and D, on the other hand, the nature of the background binary population had a substantial influence on τ_BH+NS. In both clusters, it appears that there is a threshold around f_b ∼ 0.2. Simulations with f_b < 0.2 all exhibit similar values of τ_BH+NS. At larger f_b, τ_BH+NS increased because single BHs were able to quickly exchange back into a binary. The composition of the binary population did not seem to impact the results much. As we explained in CSC13, the main advantage of the FIR population over the OBS population is that it allows for larger binary fractions within observational constraints. At the same value of f_b, simulations using the FIR population and the OBS population produced BH+NS binaries with equal efficiency. Use of the OPT population did, however, lead to larger τ_BH+NS. in clusters B, C, and D.

The simulations presented here also explored how changing the size of the NS and BH populations impacted the formation of BH+NS binaries. The value of τ_BH+NS responded linearly to changes in the NS retention fraction (f_ret). In every simulation with N_BH > 1, the increased size of the BH population lead to a reduction in the value of τ_BH+NS. When there was more than one BH in the simulation, it was common for the BH-binaries that we were evolving to interact with the other BHs in the cluster. Nearly half of all three-body encounters between a BH+star binary and a second BH will result in the formation of a BH+BH binary (Sigurdsson & Phinney 1993). The number of BH+BH binaries formed in each simulation is listed in Table 2. In a future paper, we will compute both the BH+BH and the BH+NS merger rates using our simulations with more than one BH. Once a BH+BH binary formed, it was nearly impossible for an NS to exchange into the system. These BHs were essentially locked up for the rest of the run, leading to the reduction in τ_BH+NS. In addition to preventing the formation of BH+NS binaries, the presence of several BHs can also result in the destruction of BH+NS binaries. In our simulations with multiple BHs, 10–50 per cent of BH+NS binaries were destroyed when another BH exchanged into the binary.

For relatively small BH populations, the value of p_BH+NS increased compared to simulations with a single BH. In these simulations, the small decline in τ_BH+NS was outpaced by the increase in N_BH. However, the growth in p_BH+NS quickly flattened out and/or turned over as N_BH was increased. Comparing simulations 14, 15, and 16 shows that p_BH+NS increased as we added more BHs to the cluster. Extrapolating the trend seen in τ_BH+NS to larger values of N_BH, we found that increasing N_BH to 100 would only boost the expected number of BH+NS binaries in this cluster from 0.051 to 0.06. Of course, one should use caution in drawing conclusions from such an extrapolation, but it seems unlikely that the presence of a substantial BH population would increase the size of a cluster's BH+NS binary population. In fact, such a large number of BHs could reduce the size of the BH+NS binary population. This behaviour is seen in some of our simulations. In simulations 12 and 13, and 31 and 32, an increase in N_BH resulted in a reduction of p_BH+NS. When there were many BHs in a cluster, the production of BH+BH binaries was favoured over the production of BH+NS binaries.

We construct our best estimate for the total number of BH+MSP binaries by generating several Monte Carlo realizations of the Milky Way GC population and counting the number of BH+MSP binaries in each realization. In these models, we used N_GC = 150. The mass of each cluster was drawn from the GC mass function presented in McLaughlin & Pudritz (1996). We then assigned each cluster a p_BH+NS by randomly selecting one of our simulations that was done in a GC of similar mass. Next, we chose f_MSP for each cluster from a normal distribution with a mean of 0.13 and a standard deviation of 0.07. We chose to cast our estimates in terms of the fraction of GCs that retain at least one BH (f_BH) because this quantity is so poorly constrained. Based on 10⁴ realizations, we found N_BH+MSP = (0.6 ± 0.2)f_BH. In computing this number, we assumed that the size of each GC's BH population was random. It is possible that the sizes of the BH populations are correlated, i.e. the GCs that retain BHs either all retain several or all retain ∼1. If we recompute our Monte Carlo realizations and require that N_BH > 1, then we find N_BH+MSP = (1.3 ± 0.3)f_BH. If we only consider simulations that used N_BH = 1, the size of the BH+MSP binary population is reduced to N_BH+MSP = (0.2 ± 0.1)f_BH. We have not accounted for the fact that many components of these estimates (e.g. N_BH, f_BH, and f_MSP) are likely to be functions of the GCs’ structural parameters. However, it is unlikely that including these dependences will significantly alter our conclusion that N_{BH + MSP} ≲ 1. This estimate suggests that dynamically formed BH+MSP binaries in GCs may be even rarer than those that are likely to be produced through standard binary evolution in the disc of the Galaxy.

5 DISCUSSION AND CONCLUSIONS

We have presented a study of the dynamically formed BH+MSP binaries in GCs. We found that in the highest density clusters (n_c ≳ 5 × 10⁵ pc^{− 3}), the semimajor axis distribution of the BH+MSP binaries is nearly independent of all of the parameters that we varied in our study. This property of the BH+NS binary populations is beneficial for observers who hope to identify such systems. Regardless of the nature of many uncertain characteristics, including the GC BH and binary populations, the vast majority of BH+MSP binaries produced in dense GCs will have 2 < P_B < 10 d. In lower density clusters, M_BH does influence the expected orbital periods of the BH+MSP binaries. In clusters with n_c ∼ 10⁵ pc⁻³ BH+MSP binaries with massive BHs (M_BH = 35 M_⊙) will typically have orbital periods around 20 d. For BH+MSP binaries with 7 M_⊙ BHs, the expected orbital periods are much longer, with typical periods in the 150–250 d range.

Importantly, we have also found that dynamically formed BH+MSP binaries are quite rare. We estimated that the maximum number of detectible BH+MSP binaries produced through this channel in the Milky Way GC system is approximately 3–5. Comparing the size of the BH+MSP binary population predicted by our models to population synthesis models of such binaries in the field, we find that the dynamical encounters result in a factor of ∼100 enhancement in BH+MSP binary production per unit mass in GCs. The birthrates of other exotic objects (e.g. LMXBs, MSPs) in GCs receive a similar boost over the field due to the additional, dynamical formation channels available to the members of a dense stellar system.

The small size of the population is not a consequence of our assumption that most stellar mass BHs are ejected from the cluster early in its evolution. The presence of many BHs will also reduce the probability that a cluster harbours a BH+MSP binary. BH+MSP binary formation can be stifled by as few as 19 BHs. If there are several BHs in the cluster the BHs will preferentially interact with each other and not the NSs. Furthermore, any BH+NS binaries that are formed may be destroyed when another BH exchanges into the binary. This behaviour has also been seen in models that considered the evolution of the BH population as a whole. Sadowski et al. (2008) and Downing et al. (2010) found that very few BH+NS binaries were produced in their simulations, which included several hundred to over one thousand BHs. We expect dynamically formed BH+MSP binaries to be rare regardless of the size of the retained BH population.

Another factor that played a surprisingly small role in limiting the size of the BH+MSP binary population were the large post-merger recoil that BH+NS binaries are expected to receive. Several recent studies have used numerical relativity to simulate BH+NS binary mergers, and these models show that the remnant BH will receive a kick of more than 50 km s⁻¹ when 3 ≲ M_BH/M_NS ≲ 10 (e.g. Etienne et al. 2009; Shibata et al. 2009; Foucart et al. 2011). These kicks exceed the escape velocities of all but the most massive GCs, so BHs of ∼7 M_⊙ will be ejected from the cluster once they merge with an NS. These post-merger ejections reduce τ_BH+NS because they act to remove single BHs from the cluster. However, there are three ways to avoid the large recoils. First, at smaller mass ratios, the NS is tidally disrupted before the merger, which halts the anisotropic emission of gravitational radiation and suppresses the kick. Unfortunately, it is unlikely that a BH+NS binary with such a small mass ratio exists in nature. Secondly, the linear momentum flux responsible for the kick declines for larger values of M_BH/M_NS, again reducing the magnitude of the kick (see, e.g. Fitchett 1983). Finally, a large BH spin could also decrease the magnitude of the post-merger kick. Foucart et al. (2013) showed that for BHs with dimensionless spin parameters of 0.9, the recoil would be smaller than typical GC escape velocities for BH+NS binaries with mass ratios as small as 7. We tested how these latter two, plausible scenarios would affect our results.

Many of the models discussed above already illustrate the case in which the post merger kick is suppressed because M_BH ≫ M_NS. In the simulations that used 35 M_⊙ BHs, the kick is small enough that most of the BHs are retained by the cluster after merging with an NS. However, as we previously discussed in Section 4, the BH+NS binaries produced in these models do not live as long as the binaries produced in simulations with lower mass BHs. The presence of a more massive BH accelerates both of the hardening mechanisms that drive the BH+NS binaries to merge. We also tested a more extreme mass ratio by running a simulation with a 100 M_⊙ BH. We found τ_BH+NS = 0.036 in this model, which is similar to the value of τ_BH+NS for a simulation in the same background cluster with a 35 M_⊙ BH. Furthermore, the IMBH+MSP formation rate implied by this value τ_BH+NS is consistent with previous work on the formation of such binaries by Devecchi et al. (2007). Increasing the mass of the BH to suppress the kick does not significantly increase τ_BH+NS, and accordingly the size of the BH+MSP binary population.

To test the high spin scenario, we reran simulations using one 7 M_⊙ BH in the high-density clusters (n_c ≥ 5 × 10⁵ pc⁻³) with the post-merger kicks switched off. In both cases, τ_BH+NS increased by nearly a factor of 2, to 0.13 and 0.08 in clusters C and D, respectively. Despite the significant increase in τ_BH+NS, we still find N_{BH + MSP} ≲ 1 in the Milky Way GC system. It should also be noted that these new simulations actually overestimate the number of rapidly spinning BHs retained by the GCs. The kick is only reduced significantly if the misalignment between the BH's spin and the angular momentum of the BH+NS binary is ≲ 60°. For dynamically formed binaries, it is likely that the orbital angular momentum and the spins of their components will have random orientations. Furthermore, we also counted the number of times the BHs in these simulations merged with non-NSs to estimate the degree to which the BHs would be spun up by thin disc accretion. We will discuss these results more broadly in the context of X-ray binary production in a future publication. For now, we will only examine whether the BHs are able to accrete a substantial amount of angular momentum. In both of the background clusters considered, the BHs merged with an average of 0.97 non-NSs during their lifetimes. It was only in rare cases that the BHs merged with 7–9 stars. Thus, it is unlikely that the BH's spin will increase substantially during its evolution in the cluster. The BHs must be born with large spins for this post-merger recoil mechanism to be effective. Even if this is the case, we expect that BH+MSP binaries will be extremely rare in GCs.

Some limitations of our method will impact the results of our simulations. Because our simulations do not include binary–binary encounters, they do not capture several processes that affect the formation of BH+MSP binaries. As discussed above, binary–binary interactions open up additional BH+NS binary formation channels. Furthermore, in clusters with multiple BHs, collisions between pairs of BH+BH binaries could eject or disrupt many BH+BH binaries (e.g. O'Leary et al. 2006; Banerjee, Baumgardt & Kroupa 2010; Downing et al. 2010). Reducing the number of BHs in the cluster and freeing BHs from otherwise impenetrable BH+BH binaries would increase the likelihood that BH+MSP binary is produced. However, binary–binary interactions will also disrupt and eject BH+MSP binaries. Models that include binary–binary interactions are needed to see which processes dominate.

We also neglected long-range interactions between BH-binaries and background stars. These interactions do no perturb the binary's orbital parameters as strongly as close encounters, but they do occur more frequently. Of particular concern is whether the change in eccentricity resulting from these encounters will accelerate the rate of orbital contraction through the emission of gravitational waves (see equation 3). For eccentric binaries, the change in eccentricity induced by an encounter declines as |$r_p^{-3/2}$|⁠, where r_p is separation between the binary and single star at pericentre (Heggie 1975). The change in eccentricity declines even faster with increasing r_p for circular binaries (Hut & Paczynski 1984; Rappaport, Putney & Verbunt 1989; Phinney 1992; Rasio & Heggie 1995). These encounters drive a random walk in eccentricity because they are just as likely to increase a binary's eccentricity as they are to decrease it. We used the cross-sections for eccentricity change derived by Heggie & Rasio (1996) to estimate the root mean square rate of change in eccentricity induced by the distant encounters that our models did not include. Given these rates, we found that the eccentricity of the binaries in our simulations would change by <0.05 over their entire lifetimes. In simulation 10 (n_c = 10⁵ pc⁻³, one 7 M_⊙ BH), for example, the median eccentricity change amongst the BH+NS binaries was 0.008. Allowing for this modest change in eccentricity changed the median gravitational wave merger time for the BH+NS binaries by ±7 per cent. The median merger time for BH+WD binaries changed by ±9 per cent. Distant encounters will only have a small effect on the BH+MSP binary population.

Finally, as a consequence of the assumption that the background cluster was static, our simulations were unable to capture some dynamical processes. In the static cluster models, we forced the BHs to remain in equilibrium with the rest of the cluster. If we had allowed for the dynamical evolution of the BHs, they might have decoupled from the cluster and produced a dense, inner subcluster. This would have further reduced the number of encounters between BHs and NSs. Alternatively, heating of the cluster by the BHs could result in expansion of the core, which would result in longer lifetimes for any BH+MSP binaries that managed to form (Mackey et al. 2008; Heggie & Giersz 2014). Primordial binaries are an additional source of heating that we were unable to include in our models. We showed that efficient BH+MSP binary formation required substantial binary fractions, but we did not account for the impact heating by these large binary populations could have on the structure of the cluster. Clearly, more detailed models are needed to determine how these additional processes impact the BH+MSP populations in GCs.

Although the number of BH+MSP binaries in the Milky Way GC system is expected to be small, searching for these binaries is still warranted. We know that these binaries are in GCs, and our models make specific predictions about the types of GCs that are likely to harbour BH+MSP binaries. Given the rarity of BH+MSP binaries predicted by our models, the discovery of a BH+MSP pulsar binary in the Milky Way GC system would imply that the fraction of clusters that retain at least one stellar mass BH is large. With the potential scientific payoff, continued deep radio observations of the cores of the ∼20 Milky Way GCs with appropriate structural properties may be justified. The clusters most likely to host a BH+MSP binary include 47 Tuc, Terzan 5, NGC 1851, NGC 6266, NGC 6388, and NGC 6441. Intriguingly, previous theoretical and observational studies have suggested that NGC 6388 and NGC 6441 may harbour BHs (Lanzoni et al. 2007; Moody & Sigurdsson 2009). Finally, even though there might not be any BH+MSP binaries in the Milky Way GC system, such binaries could be detected in extragalactic GCs with the Square Kilometre Array (SKA). SKA should be able to detect most pulsars within 10 Mpc (Cordes 2007), and our models predict there could be ∼100 dynamically formed BH+MSPs binaries within this volume.

We thank E. S. Phinney, Francois Foucart, and Lawrence E. Kidder for valuable discussion. We also thank the referee for an insightful review that helped us improve the paper. SS and DFC thank the Kavli Institute for Theoretical Physics for their hospitality. SS thanks the Aspen Center for Physics.

REFERENCES