Low-mass X-ray binaries from black-hole retaining globular clusters
MMNRAS , 1–27 (2017) Preprint 8 October 2018 Compiled using MNRAS L A TEX style file v3.0
Low-mass X-ray binaries from black-hole retaining globular clusters
Matthew Giesler, (cid:63) Drew Clausen, Christian D Ott TAPIR, Walter Burke Institute for Theoretical Physics, California Institute of Technology, Pasadena, CA 91125, USA
ABSTRACT
Recent studies suggest that globular clusters (GCs) may retain a substantial population ofstellar-mass black holes (BHs), in contrast to the long-held belief of a few to zero BHs.We model the population of BH low-mass X-ray binaries (BH-LMXBs), an ideal observableproxy for elusive single BHs, produced from a representative group of Milky Way GCswith variable BH populations. We simulate the formation of BH-binaries in GCs throughexchange interactions between binary and single stars in the company of tens to hundredsof BHs. Additionally, we consider the impact of the BH population on the rate of compactbinaries undergoing gravitational wave driven mergers. The characteristics of the BH-LMXBpopulation and binary properties are sensitive to the GCs structural parameters as well asits unobservable BH population. We find that GCs retaining ∼ ∼
150 ejected BH-LMXBs whereas GCs retaining only ∼
20 BHs produce zeroejected BH-LMXBs. Moreover, we explore the possibility that some of the presently knownBH-LMXBs might have originated in GCs and identify five candidate systems.
The fate of the population of stellar-mass black holes (BH) in glob-ular clusters (GCs) is still widely uncertain. It is expected that tensto hundreds and possibly thousands of BHs are formed in GCs, ofwhich some fraction might be ejected early due to a kick at forma-tion (Belczynski et al. 2006). In the standard GC evolution picture,the remainder of the BHs should rapidly sink to the core due tomass segregation. There they are subject to a high rate of dynamicalinteractions that are likely to eject the BHs as singles or in bina-ries. It was long accepted that this process would lead to repeatedejections from the GC leaving a few to zero BHs (e.g., Kulkarniet al. 1993; Sigurdsson & Hernquist 1993). Historically, this wassupported by the lack of observational evidence for a BH in a GC;however, BHs are difficult to observe unless they are actively ac-creting from a stellar companion.In order to explore the population of BHs within and outside ofGCs, black-hole low-mass X-ray binaries (BH-LMXBs) can serveas an ideal proxy. In an evolved cluster, a main-sequence star (MS)will necessarily be less than the MS turnoff mass, yielding an abun-dance of potential low-mass companions. This, coupled with a highrate of encounters due to the high-density environment of GCs,makes GCs ideal BH-LMXB factories. However, this assumes thata significant number of BHs are retained by GCs and that the BHsavoid segregating completely from the lower-mass stars.The discovery of two BH-LMXB systems in the Milky WayGC M22 (Strader et al. 2012) has led to a renewed interest in GCBH retention. This observation coupled with an estimate for thefraction of the BH population expected to be in accreting bina- (cid:63)
Email: [email protected] ries (Ivanova et al. 2010) suggests that M22 may contain between5 −
100 BHs (Strader et al. 2012). Recent theoretical studies, includ-ing some detailed N -body simulations (e.g., Aarseth 2012; Wanget al. 2016), support the idea that GCs are capable of retaining froma few to hundreds of BHs (e.g., Breen & Heggie 2013; Morscheret al. 2013; Sippel & Hurley 2013; Rodriguez et al. 2016b).There is an increasing number of BH-LMXB candidates iden-tified in the Milky Way galaxy. BlackCAT (Corral-Santana et al.2016), a catalog of BH-LMXBs, has to date identified 59 candidateMilky Way BH-LMXBs. An LMXB is identified as a candidateBH-LMXB if the X-ray spectrum rules out a neutron star (NS) asthe compact accretor (McClintock & Remillard 2006). Of the 59candidate BH-LMXBs in BlackCAT, 22 are currently considered tobe ‘confirmed’ BH-LMXBs. A BH-LMXB labeled as ‘confirmed’has a dynamical measurement of the primary mass or mass-function f ( M BH ) (see, e.g., Casares & Jonker 2014).Roughly one-fifth of the observed BH-LMXBs reside at anabsolute distance | z | perpendicular to the galactic plane greaterthan 1 kpc (e.g., Jonker & Nelemans 2004; Corral-Santana et al.2016). The distribution of the candidate and confirmed BH-LMXBswithin the Milky Way gives rise to the idea that BHs might besubject to high-velocity kicks at formation (e.g., Gualandris et al.2005; Fragos et al. 2009; Repetto et al. 2012; Repetto & Nelemans2015). In some cases, the velocity needed for the binary to reachlarge | z | exceeds the contribution from a Blaauw kick (Blaauw1961). This is the velocity imparted to a binary in the case ofsudden mass loss, i.e. in the BH progenitor’s supernova explosion.The exceptional high-velocity BH-LMXB cases have led to the ideaof high-velocity formation kicks, also known as ‘natal’ kicks, wherethe binary receives a large kick through an asymmetric explosionlaunched prior to BH formation (Janka 2013; Janka 2017). Due to © 2017 The Authors a r X i v : . [ a s t r o - ph . H E ] A ug Giesler, Clausen, & Ott the long-held assumption that GCs maintain a near-zero populationof BHs, the possibility that some of these systems originated in GCshas been largely ignored. BH-LMXBs sourced by BH-retaining GCsmight help to explain some of the peculiar properties of the observedMilky Way BH-LMXB population. Although GCs are not likely todescribe the entire population of BH-LMXBs, the halo-orbits ofGCs in the Milky Way make GCs ideal candidate sources for thehigh- | z | systems. In light of the recent studies that suggest GCsmight harbor a large number of BHs, we revisit in this paper thepossibility of GCs as a potential origination point for a subset of theobserved BH-LMXB systems.In addition to using BH-LMXBs as probes of BH retentionin GCs, the BH-BH merger rates might also serve to place someconstraints on GC BH retention. The recent success in observingmerging BH-BH binaries by advanced LIGO (aLIGO) makes thisa realistic possibility (Abbott et al. 2016a; Abbott et al. 2016b; Ab-bott et al. 2016c). Furthermore, binary BH mergers occurring inGCs may be characteristically eccentric due to dynamical forma-tion channels. Although these eccentric systems are likely to havecircularized by the time they are visible in the aLIGO frequencyband, the eccentricity is potentially detectable at lower frequencies.The addition of a space-based gravitational wave observatory (e.g.,LISA) in the future, designed for sensitivity at lower frequencies,further improves the prospect of using BH-BH mergers to probe GCdynamics.In this study, we explicitly evolve ‘test’ binaries in a fixedcluster background subject to dynamical friction and single-binaryinteractions. Additionally, we include an updated prescription forallowing single BHs to exchange into existing binaries. The GCs arechosen to represent a realistic subset of Milky Way GCs with varyingBH populations in order to investigate the effects of BH retentionin clusters. Each GC background is described by an isotropic multi-mass King model. We produce a large number of realizations foreach set of initial parameters to obtain statistical distributions of thenumber of ejected binaries and their relevant properties. Using thestatistics from the GC simulations, we then perform Monte Carlosimulations to obtain a population of BH-LMXBs produced by GCs.The GCs and the ejected binaries are evolved in time through theMilky Way potential while simultaneously accounting for the stellarevolution of the ejected binaries. The resulting mass-transferringsystems make up a previously unexplored galactic population of BH-LMXBs from GCs. We investigate the distribution and propertiesof the resulting population and its dependence on BH retention inGCs. Specifically, we find that in the case of minimal BH retention( N BH =
20) no observable BH-LMXBs are produced, while the N BH =
200 and N BH = + − and 156 + − BH-LMXBs. Furthermore, we usethe resulting population to determine the most likely candidates fora GC origin in the population of observed Milky Way BH-LMXBs:the five systems that are compatible with our simulated population ofBH-LMXBs from GCs are SWIFT J1357.2-0933, SWIFT J1753.5-0127, XTE J1118+480, and GRO J0422+32. Future measurementswill be necessary to increase support for a GC origin theory, butif we can confidently attribute a BH-LMXB to a GC, this wouldprovide strong evidence for significant BH retention in GCs.The remainder of this paper is structured as follows. In sec-tion 2, we describe our model for the GCs and the evolution of atest-binary in a static cluster background. In section 3, we lay outhow we generate the present-day BH-LMXB population from oursimulations of Milky Way GCs. In section 4, we review the prop-erties of the ejected BH binaries along with the distribution andproperties of the present-day BH-LMXBs from GCs. Additionally, we explore the effects of BH retention on the BH-BH merger ratein GCs. We conclude the section by comparing our results withobservations and previous work. Finally, in section 5, we provideconcluding remarks.
GCs typically contain ∼ − stars, which makes them accessi-ble to modern N -body simulations (e.g., Zonoozi et al. 2011; Wanget al. 2016) that can track GC evolution. However, full N -body clus-ter evolution simulations are still very computationally expensive,making this method poorly suited for studying many realizations ofdifferent GCs necessary for building statistics on the evolution of BHbinaries inside clusters. Fokker-Planck methods are more approx-imate and describe GCs with a phase-space distribution functionfor its constituent stars that evolves via the Fokker-Planck equation,a Boltzmann equation with a small local collision term that mod-ifies only velocities (see, e.g., Spitzer 1987). The Fokker-Planckequation can be numerically integrated directly (e.g., Cohn 1979;Chernoff & Weinberg 1990) or, more commonly, integrated withMonte Carlo methods (see, e.g., Hénon 1971; Spitzer & Hart 1971and Rodriguez et al. 2016b for a comparison between N -body andthe Monte Carlo approaches). However, here we are concerned withthe evolution of BH binaries in GCs and not with the GC evolu-tion itself. Hence, we adopt the approach of Sigurdsson & Phinney(1995) and model the evolution of binaries in a fixed cluster back-ground. We approximate the collision term in the Fokker-Planckequation analytically to model the effects of distant encounters asthe binary evolves through the GC. Near encounters are accountedfor by explicitly integrating the three-body equations of motion.We build up statistics by carrying out simulations of many randomrealizations of binaries for a given GC background model. In thefollowing sections, we describe our method in detail. Our model, based on Sigurdsson & Phinney (1995), incorporates anumber of assumptions that simplify the simulations and allow usto perform ∼ realizations for a given cluster model with rela-tively minimal computational needs. The three key assumptions are:(i) GCs are well described by a ‘lowered Maxwellian’ distributionfunction, (ii) the gravitational potential and distribution functionsare stationary, and (iii) the effect of distant interactions is well de-scribed by the leading order terms in the Fokker-Planck equation.The ‘lowered Maxwellian’ distribution function, which eliminatesthe tail of the Maxwellian velocity distribution, introduces a max-imum energy for stars within the cluster to remain bound. Thismaximum energy φ ( r t ) implies a finite mass and a maximum radius r t , commonly referred to as the ‘tidal’ radius, as stars beyond thisdistance are pulled from the cluster by the galactic tidal field. Mod-els based on a ‘lowered Maxwellian’, commonly referred to as Kingmodels, readily describe many observed clusters (Peterson & King1975; Bahcall & Hausman 1977; Spitzer 1987).We evolve a single ‘test binary’, initialized according to sec-tion 2.2.5, in a static cluster background described by an isotropicmulti-mass King model (King 1966) defined by single particle dis-tribution functions f α ( r , v , m α ) for a discrete set of mass groups.Here, r and v are the radius and velocity in the cluster center-of-mass frame and m α is the representative mass of group α . The MNRAS , 1–27 (2017)
H-LMXBs from BH retaining GCs distribution function for a given mass group is given by f α ( ε ) = (cid:40) n α ( πσ α ) / ( e − ε / σ α − ) ε < ε (cid:62) . (1)Here, ε is the energy per unit mass, ε = v / − Ψ ( r ) , and Ψ ( r ) ≡ φ ( r t ) − φ ( r ) is the gravitational potential relative to thatat the tidal radius r t . Additionally, σ α is the group’s velocity dis-persion at the core of the cluster and n α is a normalization fac-tor. For an isotropic cluster, the velocity dispersion reduces to theone-dimensional mean-square velocity, such that 3 σ α = ¯ v α . Thenormalization factor in its full form is n α = η α n o e Ψ ( )/ σ α erf (cid:18)(cid:114) Ψ ( ) σ α (cid:19) − (cid:114) Ψ ( ) πσ α (cid:18) + Ψ ( ) σ α (cid:19) , (2)where η α = N α / N is the number fraction for mass group α and n o = n ( ) is the central density.The free structural parameters necessary to specify a modelcluster, with specified mass groups, are the mean core veloc-ity dispersion ¯ σ , the core number density n o , and the potentialdepth, which is specified by the dimensionless King parameter W o = Ψ ( )/ ¯ σ . The remaining structural parameters, which arefully determined by the free parameters, are: total mass M c , core ra-dius r c , tidal radius r t , and concentration c = log ( r t / r c ) . The coreradius r c is defined as the radius at which the surface brightness hasdropped to half the value at the core.For a given set of masses with corresponding distribution func-tions, the cluster satisfies Poisson’s equation for the relative poten-tial ∇ Ψ ( r ) = − π G (cid:205) α ρ α . Here, ρ α = m α n α , where n α is thenumber density of mass group α given by n α = ∫ v ( r t ) f α ( r , v , m α ) π v d v . (3)The upper limit of the integral is the maximum allowed velocity v ( r t ) = (cid:112) Ψ ( r t ) , i.e. the escape velocity. The object masses m α andnumber fraction η α are determined by the evolved mass function,discussed in section 2.2.1. We generate a model cluster that satisfiesPoisson’s equation for the specified masses and number fractions inan iterative fashion. We begin by integrating Poisson’s equation outto a radius r t , implicitly determined by Ψ ( r t ) =
0, with boundaryconditions Ψ ( ) = W o and ∇ Ψ ( ) =
0, and take η α = η α asour initial guess. The actual number fraction of each mass group, η α = N α / N , is then calculated using N α = ∫ r t n α ( r ) π r dr , (4)along with N = (cid:205) α N α . We then update our guess to η α = ( η α new + η α old )/
2, where η α new → η α old × ( η α / η α ) . We repeatthe above steps until ( η α − η α )/ η α < δ is satisfied for all massgroups, where we have made the somewhat arbitrary choice of δ = . × − for our convergence threshold. This iterative pro-cedure determines the normalization constant n α and r t . Once r t is found, the concentration c = log ( r t / r c ) is determined and thetotal mass of the cluster M c is obtained from −∇ Ψ ( r t ) = GMr . (5)The evolution of our ‘test binary’ in the cluster background isaffected by long-range and short-range interactions, which modifythe magnitude and direction of the binary’s velocity. The short-range encounters are accounted for by fully resolving the three-body interactions, detailed in section 2.3.4. We account for the velocity fluctuations due to long-range interactions with ‘field stars’,distant cluster stars, through the diffusion coefficients D ( ∆ v i ) and D ( ∆ v i ∆ v j ) in the Fokker-Planck equation, D fDt = (cid:32) ∂ f ∂ t (cid:33) enc = (cid:213) i , j (cid:26) − ∂∂ v i ( D ( ∆ v i ) f ) + ∂ ∂ v i ∂ v j ( D ( ∆ v i ∆ v j ) f ) (cid:27) . (6)In this context, a diffusion coefficient D ( X ) for a variable X , cor-responds to the average change in X per unit time. Here, we focuson velocity changes per unit time as experienced by the binary dueto interactions with the ‘field stars’. The form of the coefficientscan be derived, for a simple case, by first considering the change invelocity of a mass m , initially at rest, due to an encounter with asecond mass m at a relative velocity v with impact parameter p , ( ∆ v ) = m ( m + m ) v ( + ( pp o ) ) , (7)where p o ≡ G ( m + m )/ v is a reference impact parameterwhich causes a deflection of π /
2, consistent with close encoun-ters (e.g., Spitzer 1987). The average rate of change of the quantityin Equation 7, per unit time, due to encounters is then obtained byintegrating over the possible impact parameters for a given densityof field stars n , D ( ∆ v ) = π ∫ p max ∆ v pn v dp , (8)up to a maximum allowable impact parameter p max . The maxi-mum impact parameter is required to suppress the divergence ofthe integral and essentially determines the maximum distance oflong-range encounters that contribute to the velocity perturbations.This maximum value, p max , is not explicitly specified, but finds itsway into the coefficient calculations through the so-called Coulomblogarithm, ln Λ ≡ ln ( p max / p o ) , which appears as a result of theintegration.We work out the details for the case of an isotropic velocitydispersion with a density of field stars given by Equation 3 andrestate the relevant coefficients we use in our model (cf. Binney &Tremaine 2008). These coefficients, which describe the average rateof change in the velocity of the binary due to long-range encounters,are used to update the velocity of the binary at each time step.The implementation is described further in section 2.3. A detailedderivation and a more general form of the coefficients can be foundin Spitzer (1987).By choosing a coordinate system in which one axis is alignedwith the velocity of the binary, we can decompose D ( ∆ v i ) into acoefficient parallel to the binary’s velocity D ( ∆ v (cid:107) ) and two mutuallyorthogonal coefficients perpendicular to the velocity, D ( ∆ v ⊥ ) and D ( ∆ v ⊥ ) . In an isotropic cluster, there is no preferred direction withregard to the two perpendicular components, so the contributionsfrom D ( ∆ v ⊥ ) and D ( ∆ v ⊥ ) tend to cancel each other out. Theirsquares, D ( ∆ v ⊥ ) and D ( ∆ v ⊥ ) , on the other hand, do not and arenon-vanishing. Additionally, we include a quadratic term for theparallel component D ( ∆ v (cid:107) ) and in consideration of the symmetrywe retain only the sum of the perpendicular components D ( ∆ v ⊥ ) = D ( ∆ v ⊥ ) + D ( ∆ v ⊥ ) .The diffusion coefficient D ( ∆ v (cid:107) ) parallel to the binary’s motionis by analogy often referred to as the coefficient of dynamical frictionas it opposes the binary’s direction of motion, D ( ∆ v (cid:107) ) = − (cid:213) α γ α (cid:18) + m b m α (cid:19) ∫ v (cid:18) v α v (cid:19) f α ( v α ) d v α . (9) MNRAS000
2, consistent with close encoun-ters (e.g., Spitzer 1987). The average rate of change of the quantityin Equation 7, per unit time, due to encounters is then obtained byintegrating over the possible impact parameters for a given densityof field stars n , D ( ∆ v ) = π ∫ p max ∆ v pn v dp , (8)up to a maximum allowable impact parameter p max . The maxi-mum impact parameter is required to suppress the divergence ofthe integral and essentially determines the maximum distance oflong-range encounters that contribute to the velocity perturbations.This maximum value, p max , is not explicitly specified, but finds itsway into the coefficient calculations through the so-called Coulomblogarithm, ln Λ ≡ ln ( p max / p o ) , which appears as a result of theintegration.We work out the details for the case of an isotropic velocitydispersion with a density of field stars given by Equation 3 andrestate the relevant coefficients we use in our model (cf. Binney &Tremaine 2008). These coefficients, which describe the average rateof change in the velocity of the binary due to long-range encounters,are used to update the velocity of the binary at each time step.The implementation is described further in section 2.3. A detailedderivation and a more general form of the coefficients can be foundin Spitzer (1987).By choosing a coordinate system in which one axis is alignedwith the velocity of the binary, we can decompose D ( ∆ v i ) into acoefficient parallel to the binary’s velocity D ( ∆ v (cid:107) ) and two mutuallyorthogonal coefficients perpendicular to the velocity, D ( ∆ v ⊥ ) and D ( ∆ v ⊥ ) . In an isotropic cluster, there is no preferred direction withregard to the two perpendicular components, so the contributionsfrom D ( ∆ v ⊥ ) and D ( ∆ v ⊥ ) tend to cancel each other out. Theirsquares, D ( ∆ v ⊥ ) and D ( ∆ v ⊥ ) , on the other hand, do not and arenon-vanishing. Additionally, we include a quadratic term for theparallel component D ( ∆ v (cid:107) ) and in consideration of the symmetrywe retain only the sum of the perpendicular components D ( ∆ v ⊥ ) = D ( ∆ v ⊥ ) + D ( ∆ v ⊥ ) .The diffusion coefficient D ( ∆ v (cid:107) ) parallel to the binary’s motionis by analogy often referred to as the coefficient of dynamical frictionas it opposes the binary’s direction of motion, D ( ∆ v (cid:107) ) = − (cid:213) α γ α (cid:18) + m b m α (cid:19) ∫ v (cid:18) v α v (cid:19) f α ( v α ) d v α . (9) MNRAS000 , 1–27 (2017)
Giesler, Clausen, & Ott
Here, m b is the mass of the binary and γ α ≡ ( π Gm α ) ln Λ , wherewe have chosen to set ln Λ =
10, a value typical for GCs (Spitzer1987). The two remaining coefficients, D ( ∆ v (cid:107) ) = (cid:213) α v γ α (cid:40) ∫ v (cid:18) v α v (cid:19) + ∫ ∞ v (cid:18) v α v (cid:19)(cid:41) f α ( v α ) d v α (10)and D ( ∆ v ⊥ ) = (cid:213) α v γ α × (cid:40) ∫ v (cid:20) (cid:18) v α v (cid:19) − (cid:18) v α v (cid:19) (cid:21) + ∫ ∞ v (cid:18) v α v (cid:19)(cid:41) f α ( v α ) d v α , (11)are strictly positive. These coefficients are responsible for thestochastic perturbations to the parallel and perpendicular compo-nents of the velocity, which take the binary on a random walkthrough velocity space and compete with the slowing due to dy-namical friction. We implement these ‘random kicks’ as discretechanges to the binary’s velocity by sampling from a normalizeddistribution of the velocity perturbations, described in section 2.3. We obtain an initial distribution of masses in the range 0 . M (cid:12) < m < M (cid:12) from the broken-power-law initial mass function(IMF) ξ ( m ) ∝ (cid:40) m − . m . − x ∗ x m < m x m − . − x ∗ m (cid:62) m x , (12)with x ∗ = .
35 and m x = . M (cid:12) chosen to incorporate a SalpeterIMF (Salpeter 1955) for masses above m x and a Kroupa ‘correc-tion’ (Kroupa 2001) to masses below m x along with a normalizationfactor for continuity. Stars with masses below the main-sequenceturn-off, which we set to m to = . M (cid:12) (Meylan & Heggie 1997),are assumed not to evolve significantly on the timescale of the sim-ulations, while masses above m to are assumed to be completelyevolved according to a specified evolved mass function (EMF). Theevolved mass m e is determined by the EMF m e = m MS = m . M (cid:12) < m (cid:54) m to m WD = . + . ( m − ) m to < m < M (cid:12) m NS = . M (cid:12) (cid:54) m < M (cid:12) m BH = m BH ( m , f s BH ) M (cid:12) < m < M (cid:12) , (13)where the mass subscripts label the object type and refer to mainsequence (MS), white dwarf (WD), neutron star (NS), and black hole(BH). We occasionally refer to the set of MS and WD objects as thenon-compact (NC) population. The MS stars below the turnoff massare set to their zero-age main-sequence (ZAMS) mass, the WD starsare a linear function of their ZAMS mass (Catalán et al. 2008), NSare simply set to 1 . M (cid:12) . Following the work of Sana et al. (2012),the BHs are assumed to have formed from two possible channels:stars with companions that significantly affect the evolution of thestar and those stars that are ‘effectively single.’ Effectively singleis used to describe stars that evolve in isolation as well as those M ZAMS [ M (cid:12) ] M H e [ M (cid:12) ] Figure 1.
The He core mass (marked by circles) as a function of zero-agemain-sequence mass from the
MESA (Paxton et al. 2011) runs, along with thefit (blue, dashed line) given by Equation 14. For the ∼
70 per cent of BHsformed in binaries, we approximate the remnant BH mass with the He coremass of the progenitor. The remnant mass for the remaining ∼
30 per cent ofBHs is approximated by 0 . M ZAMS , which accounts for the hydrogen masslost to stellar winds at low metallicity. stars that evolve in wide binaries with minimal interaction. Sanaet al. (2012) estimate that ∼
70 per cent of massive stars will havetheir final state impacted by a companion, which motivates setting f s BH = . ∼
10 per cent of the initial massand set m e = . m . The remaining 70 per cent of BHs formed willhave evolved with a companion and likely passed though a commonenvelope phase, stripping the stars down to their helium (He) cores(Sana et al. 2012; de Mink et al. 2014). Using MESA (Paxton et al.2011) to evolve masses in the range 20 M (cid:12) < m < M (cid:12) , weobtained the He core mass as a function of the ZAMS mass in orderto determine the remnant mass for the remaining ( − f s BH ) fractionof BHs m e = m He = . (cid:0) m ZAMS (cid:1) . M (cid:12) . (14)The stellar evolution performed using MESA version 6794, followsthe procedure laid out in Morozova et al. (2015). Figure 1 displaysthe resulting He core mass as a function of the ZAMS mass from the
MESA runs with metallicity Z = × − , along with the power-lawfit of Equation 14. This metallicity corresponds to the higher peakin the bimodal, GC metallicity distribution (Harris 1996). We useEquation 14 for all clusters, since the He core masses from modelswith Z = × − , corresponding to the secondary peak, differ by (cid:46)
10 per cent.In addition to specifying the evolved masses, it is also necessaryto specify the number of NS and BH objects retained by the clusterin its static state. We specify the retained population of compactobjects, comprised of NSs and BHs, through the retention fractions f r NS and f r BH , respectively. This is necessary since we are modelingthe cluster in its evolved state, a time at which many of the NSand BHs formed within the cluster have already been ejected due MNRAS , 1–27 (2017)
H-LMXBs from BH retaining GCs Mass group m min [ M (cid:12) ] m max [ M (cid:12) ] ¯ m [ M (cid:12) ] f m f n f L Table 1.
Evolved mass groups for NGC 6121 ( N BH = m min , the upper boundary bin mass m max , the average mass of the group ¯ m , the fraction of the total mass in the cluster f m , the number fraction with respect to the total number of objects in thecluster f n , and the fraction of luminous objects in the group f L . For reference, the BH masses occupy the top three mass groups with mean masses of 8 . M (cid:12) ,20 . M (cid:12) , and 57 . M (cid:12) . to formation kicks. Studies of the proper motion of pulsars suggestthat NSs receive kicks in the range of 200 −
450 km s − (Lyne& Lorimer 1994), easily exceeding the typical escape velocity ofclusters, which is on the order of tens of km s − . However, theobservations of pulsars in GCs implies a ‘retention problem,’ sincethe observed fraction retained is inconsistent with the average natalkick velocities being significantly greater than GC escape velocities.This issue is somewhat reconciled by assuming some NSs form inbinaries, which dampen the kick and allow the GC to maintain ahold on the NS and companion (Pfahl et al. 2002). In considerationof these observations, for the case of NSs, we retain a constantfraction, f r NS = .
1, of those produced by the IMF (Sigurdsson &Phinney 1995; Pfahl et al. 2002; Ivanova et al. 2008). In the BHcase, the distribution of natal kicks is highly uncertain. Rather thantake the retention fraction f r BH to be a constant across clusters, asin the NS case, we utilize this fraction as a free parameter in ourmodels to control the number of retained BHs in each modeled GC.Once we have determined the evolved masses from the IMF, themasses are binned into 12 groups. The small number of bins allowsfor a proper representation of the true distribution while keeping thecomputational costs to a minimum. Poisson’s equation is then inte-grated to determine the final structural parameters as discussed insection 2.1. For illustrative purposes, the evolved mass distributionfor NGC 6121 with 200 retained BHs is given in Table 1. The binsfor each mass group, the mean mass in each bin, and the fraction ofluminous objects are constant across simulations, however the massfraction and number fraction depend on the structure of the clusterand the number of BHs. As discussed in section 2.1, one of the free parameters in our modelwhen specifying a cluster’s structure is the core number density n o .However, because this parameter is not easily observable, a GC’sdensity is often reported in terms of a central luminosity density ρ L .For each mass group we determine a central luminous number den-sity n L α = f L α ¯ n α , where f L α and ¯ n α are the fraction of luminousobjects and the core density, respectively, of mass group α . The cen-tral luminosity density is then given by ρ L = (cid:205) α L α n L α . In order toaccount for the variability in the mass-luminosity relation with stel-lar mass, we use a parameterized luminosity for each group of theform L α = a ( m α ) b , with luminosity coefficients a = . , b = . m α < . M (cid:12) and a = . , b = . ρ L for eachintegrated cluster and adjust n o accordingly to match the observedquantity. In order to account for the uncertainty in the size of the binary pop-ulation within a cluster, we allow for a specifiable binary fraction.The fraction of objects that are binaries is f b = N b N s + N b , (15)where N s and N b are the number of single objects and binary objects,respectively, and the total number of objects in our model clusters isthen N = N s + N b . Observations of the binary fraction are limitedto the luminous objects within the cluster. Due to this restriction,we take the observed fraction to be determined solely by the MSstar binary fraction f obs = N MS b /( N MS s + N MS b ) , where, as above,we respectively refer to N MS s and N MS b as the number of single andbinary MS stars. Using the above definitions along with the fractionof all binaries that are MS-MS binaries, f MS b = N MS b / N b , and thefraction of objects that are MS stars, f MS = N MS / N , we convert theobserved binary fraction into a uniform total binary fraction for usein our models through the relation f b = (cid:32) f MS b f MS ( f obs + ) f obs − (cid:33) − . (16)The number of MS stars N MS is determined solely by the IMFand for the simulations in this study we use f MS b = .
23 (Fregeauet al. 2009). We perform our simulation with f obs covering a rangeof values, consistent with observational constraints, between 5 − f obs taking values from the set { . , . , . } . However, we find that this parameter has a negli-gible effect on the quantities of interest, so for conciseness, it is notspecified in the simulation parameters. Recent studies of BH retention in GCs have shown clusters initiallyretain between 65 −
90 per cent of the BHs formed in cluster, withthe remainder being lost due to formation kicks (Morscher et al.
MNRAS000
MNRAS000 , 1–27 (2017)
Giesler, Clausen, & Ott N -bodysimulation (Rodriguez et al. 2016b).In the standard King model, it is common to assume that themass groups satisfy an equipartition of energy. Specifically, m α σ α = ¯ m ¯ σ , (17)where m α and σ α are the mass and velocity dispersion of massgroup α , ¯ m is the mean mass of all objects in the cluster, and ¯ σ isthe mean velocity dispersion. However, with this equipartition ofkinetic energy amongst all mass groups, the heavier objects thennecessarily have lower random velocities compared to the lighterobjects and become trapped deep in the gravitational potential at thecore of the cluster. With an equipartition of kinetic energy in place,the much more massive BHs densely populate the central regionof the cluster, driving the core radius to a small fraction of thetidal radius. This disparity between the core radius and tidal radiusleads to concentrations that deviate from observations, limiting themodeled clusters to supporting only a small number of BHs. Inorder to generate clusters with a significant BH population that arestill representative of observed GCs, motivated by Morscher et al.(2015), we implement a velocity dispersion for the BHs away fromenergy equipartition, We maintain an equipartition of energy amongthe lower-mass objects and use a modified energy partitioning forthe BHs of the form m β σ β = (cid:205) m β (cid:205) m α f s ¯ m ¯ σ , (18)where the indices β and α label the mass groups corresponding toBHs and non-BHs, respectively. Here, f s is a specifiable scale factorof order unity. The f s parameter is enough to rescale the velocitydispersion for the BHs, however, the factor involving the mass ratiocontributes substantially and f s remains of order unity and does notvary wildly across the GCs we consider.With this modified BH velocity dispersion in place, we findthat we can match the observed structural parameters of a specificcluster for zero BHs up to ∼
20 BHs, in the case of more massiveclusters up to ∼
100 BHs, and in the most massive clusters up to ∼ f s in Equation 18 and the fraction retained, f r BH , introduced in section 2.2.1. To illustrate the spreading of theBHs, we present in Figure 2 the radial density profiles for the BHsand the non-BH objects for different populations of retained BHs inthe cluster model representing NGC 6656. In the case of minimalBH retention, the BH number density falls off quickly outside of thecore, which for our model of NGC 6656 is located at r c = .
73 pcand is marked by a vertical line in Figure 2 for reference. However,in the case of many BHs, the modified velocity dispersion extendsthe number density profile radially, spreading the BHs throughoutthe cluster, without affecting the central density. The distribution ofnon-BH objects is largely unaffected by the change in BH numbers.
We choose the initial masses for our ‘test binary’ by randomly sam-pling from the evolved mass distribution and reject those that do notcontain at least one BH. If one of the component masses falls withina mass bin with a non-zero luminous population, we then sample r [ pc ] − − − l o g ( n [ N p c − ] ) non-BH ( N BH = N BH = N BH = N BH = N BH = N BH = Figure 2.
Radial number density profiles for the BH subgroup (solid lines)and the non-BH objects (dashed lines) in NGC 6656 for the three consideredvalues of N BH . The vertical line (red, dashed), at r c = .
73 pc, marks thecore radius for this cluster. The non-BH objects are largely unaffected by thedifferent numbers of BHs added to the cluster and the necessary modificationto the velocity dispersion. For N BH =
20, the BHs are concentrated in thecore region, whereas to accommodate N BH (cid:62) from the luminous mass fraction to determine whether the low-massobject is an MS star or WD. Additionally, if the selected mass isin the turnoff group, 0 . M (cid:12) (cid:54) m (cid:54) . M (cid:12) , then the object ischosen to be a giant with probability P = . f L , where f L is theluminous fraction for the turnoff-mass group. The probability forgiants is adopted from Sigurdsson & Phinney (1995) and representsthe approximate fraction of the cluster age that giants in this massrange survive. Once the masses and object types are established,the BH radii are set to the Schwarzschild radius R BH = GM / c ,while the stellar radii are determined as described in Sigurdsson& Phinney (1995). The eccentricity of the binary e is specified bysampling from the probability density function f ( e ) = e (Jeans1919), commonly referred to as a ‘thermal’ eccentricity distribu-tion. The semi-major axis a is obtained from a distribution uniformin log a in the range − (cid:54) log ( a au − ) (cid:54)
1. To avoid animmediate merger of the objects in our initial binary, we enforce a > f tid ( R + R )/( − e ) , where R i are the radii of each componentof the binary and f tid = .
1, by letting a → a until this condition issatisfied. The factor f tid is chosen based on the separation at whichtidal effects would induce a merger (Lee & Ostriker 1986). Oncethe binary parameters are set, we sample the primary-mass num-ber density profile n α ( r ) to determine the binary placement withinthe cluster and obtain a velocity from the primary-mass velocitydistribution function at r . Once we have an appropriate model, which satisfies the structuralparameters for a specific cluster and an initial binary, we then evolvethis single binary within the cluster background. In addition to thestatic potential, we include the interaction terms discussed in sec-tion 2.1. To account for dynamical friction, the diffusion coefficient D ( ∆ v (cid:107) ) is added to the potential gradient to create a smooth ef- MNRAS , 1–27 (2017)
H-LMXBs from BH retaining GCs fective acceleration a eff = ∇ Ψ ( r ) + D ( ∆ v (cid:107) ) . This smooth force isintegrated using a 4th order Runge-Kutta integrator, which is dis-cussed in detail in section 2.3.2. The quadratic scattering terms, orrandom ‘kicks’, are implemented by discretely updating the corre-sponding velocity components at each time step ∆ t . As discussedin section 2.1, the diffusion coefficient for ∆ v , of the form D ( ∆ v ) ,represents the change in this quantity per unit time, i.e. ∆ v / ∆ t . Weupdate the velocity at each time step by sampling from the normaldistribution of kicks through ∆ v (cid:107) = X (cid:113) D ( ∆ v (cid:107) ) ∆ t , ∆ v ⊥ = Y (cid:114) D ( ∆ v ⊥ ) ∆ t , ∆ v ⊥ = Y (cid:114) D ( ∆ v ⊥ ) ∆ t , (19)where X and Y are random numbers with mean values of zero andstandard deviations of one.At each time step, we also consider the evolution of the binary’ssemi-major axis a and its eccentricity e due to gravitational wave(GW) emission. If the BH is in a binary with another compact object— which includes BHs, NSs, and WDs — then we implement theevolution of a and e according to the gravitational radiation formal-ism of Peters (1964). In these cases, we also calculate the time untilcoalescence t d due to the decay of a , and if this will occur withinthe current time step, t d < ∆ t , we consider this a GW merger. Ifthe merger is of a BH-BH or BH-NS binary, we add a recoil veloc-ity, or ‘kick’, based on the fits to numerical relativity simulationsgiven by Campanelli et al. (2007) with initial spin magnitudes andorientations assigned as in Clausen et al. (2013). As the binary moves throughout the cluster, at each time step, wecheck for the possibility of a short-range encounter with a singlestar. Since the effects of long-range interactions are accounted for bythe diffusion coefficients (section 2.1), here we focus on capturingthe effects due to strong three-body interactions with much smallerimpact parameters. We limit the range of encounters to include onlythose three-body interactions that result in a resonance, exchange,ionization, or the occasional flyby. We accomplish this by choosingthe maximum impact parameter to be p = a [ B + C ( + e )] , (20)where we have set B = C = . m α with velocity v α is σ ( v , v α ) = π p + π G ( m b + m α ) p | v − v α | , (21)(see, e.g., Spitzer 1987). We then calculate the expected encounterrate between the binary and each mass group Γ ( r , v , α ) = ∫ σ ( v , v α )| v − v α | f α ( v α ) d v α , (22) and from this assign the probability of interacting with mass group α to be P α = Γ ( r , v , α ) ∆ t . (23)An encounter is deemed to have occurred, based on a random gen-erated number Z from a uniform distribution between 0 and 1, if Z is less than the total probability P = (cid:205) α P α . The total probability isimplicitly constrained to be less than unity by controlling the time-step size ∆ t , which is discussed in more detail in the subsequentsection. In the case that Z < P , we select the third star m based onthe relative probabilities P α and initiate our three-body integrationscheme explained in section 2.3.3. We use a 4th order Runge-Kutta integrator to evolve the effectiveacceleration introduced in section 2.3 as well as the three-bodyinteractions described in section 2.3.3. During integrations, we uti-lize a time step reduction scheme requiring that the accuracy of thesolution does not vary by more than a tolerance of (cid:15) rk = − when the time step is halved. The initial integration time step ∆ t = λ ( + r )/( + v ) is dynamically determined to account forthe position and velocity of the binary in the cluster, with λ = . r c / ¯ σ for a binary at rest in the core. This time stepperaccounts for the higher density in the core and the enlarged crosssection at small velocities. Although this choice of time step is usu-ally sufficient, some extra care needs to be taken when using ∆ t inEquation 23 to determine the encounter probability, so that the totalprobability does not exceed unity. To ensure that we correctly sam-ple the encounter probabilities, by satisfying the constraint P (cid:28) P max = . P < P max by reducing the timestep ∆ t when necessary. For the case P > P max , we decrease thesucceeding time step by letting λ → . ( P max / P ) λ . During thesubsequent step, if λ < λ o , where λ o = . P < P max , we allow the time step to increase slowly by setting λ → . λ . Once λ > λ o and the probability is satisfactorily small,which often occurs once the binary migrates out of the problematicdense region, we reset the time step factor to λ = λ o . In the case of an encounter, the relative probabilities described insection 2.3.1 determine the mass and velocity of the third object.We take this sampled velocity v to be the velocity of the thirdbody at infinity and calculate the relative velocity at infinity forthe encounter from v ∞ = | v − v | = (cid:113) v + v − vv cos χ . Given v and the sampled m and v , the relative velocity at infinity isdetermined up to the cos χ term, which for an isotropic King modeldistribution function can be sampled from an analytic expression for χ ∈ [ , π ] as in Sigurdsson & Phinney (1995). With the mass of thethird body and the relative velocity known, the maximum impactparameter is obtained from the cross section for the encounter π p = σ ( v , v ) = π p (cid:18) + G ( m b + m ) p v ∞ (cid:19) , (24)with p defined in Equation 20. The actual impact parameter forthe encounter is sampled from a uniform distribution in the areaspanned by the maximum impact parameter π p . The anglesthat comprise the remaining free variables necessary to specify theinitial conditions are the projected true anomaly f of the binary at MNRAS000
H-LMXBs from BH retaining GCs fective acceleration a eff = ∇ Ψ ( r ) + D ( ∆ v (cid:107) ) . This smooth force isintegrated using a 4th order Runge-Kutta integrator, which is dis-cussed in detail in section 2.3.2. The quadratic scattering terms, orrandom ‘kicks’, are implemented by discretely updating the corre-sponding velocity components at each time step ∆ t . As discussedin section 2.1, the diffusion coefficient for ∆ v , of the form D ( ∆ v ) ,represents the change in this quantity per unit time, i.e. ∆ v / ∆ t . Weupdate the velocity at each time step by sampling from the normaldistribution of kicks through ∆ v (cid:107) = X (cid:113) D ( ∆ v (cid:107) ) ∆ t , ∆ v ⊥ = Y (cid:114) D ( ∆ v ⊥ ) ∆ t , ∆ v ⊥ = Y (cid:114) D ( ∆ v ⊥ ) ∆ t , (19)where X and Y are random numbers with mean values of zero andstandard deviations of one.At each time step, we also consider the evolution of the binary’ssemi-major axis a and its eccentricity e due to gravitational wave(GW) emission. If the BH is in a binary with another compact object— which includes BHs, NSs, and WDs — then we implement theevolution of a and e according to the gravitational radiation formal-ism of Peters (1964). In these cases, we also calculate the time untilcoalescence t d due to the decay of a , and if this will occur withinthe current time step, t d < ∆ t , we consider this a GW merger. Ifthe merger is of a BH-BH or BH-NS binary, we add a recoil veloc-ity, or ‘kick’, based on the fits to numerical relativity simulationsgiven by Campanelli et al. (2007) with initial spin magnitudes andorientations assigned as in Clausen et al. (2013). As the binary moves throughout the cluster, at each time step, wecheck for the possibility of a short-range encounter with a singlestar. Since the effects of long-range interactions are accounted for bythe diffusion coefficients (section 2.1), here we focus on capturingthe effects due to strong three-body interactions with much smallerimpact parameters. We limit the range of encounters to include onlythose three-body interactions that result in a resonance, exchange,ionization, or the occasional flyby. We accomplish this by choosingthe maximum impact parameter to be p = a [ B + C ( + e )] , (20)where we have set B = C = . m α with velocity v α is σ ( v , v α ) = π p + π G ( m b + m α ) p | v − v α | , (21)(see, e.g., Spitzer 1987). We then calculate the expected encounterrate between the binary and each mass group Γ ( r , v , α ) = ∫ σ ( v , v α )| v − v α | f α ( v α ) d v α , (22) and from this assign the probability of interacting with mass group α to be P α = Γ ( r , v , α ) ∆ t . (23)An encounter is deemed to have occurred, based on a random gen-erated number Z from a uniform distribution between 0 and 1, if Z is less than the total probability P = (cid:205) α P α . The total probability isimplicitly constrained to be less than unity by controlling the time-step size ∆ t , which is discussed in more detail in the subsequentsection. In the case that Z < P , we select the third star m based onthe relative probabilities P α and initiate our three-body integrationscheme explained in section 2.3.3. We use a 4th order Runge-Kutta integrator to evolve the effectiveacceleration introduced in section 2.3 as well as the three-bodyinteractions described in section 2.3.3. During integrations, we uti-lize a time step reduction scheme requiring that the accuracy of thesolution does not vary by more than a tolerance of (cid:15) rk = − when the time step is halved. The initial integration time step ∆ t = λ ( + r )/( + v ) is dynamically determined to account forthe position and velocity of the binary in the cluster, with λ = . r c / ¯ σ for a binary at rest in the core. This time stepperaccounts for the higher density in the core and the enlarged crosssection at small velocities. Although this choice of time step is usu-ally sufficient, some extra care needs to be taken when using ∆ t inEquation 23 to determine the encounter probability, so that the totalprobability does not exceed unity. To ensure that we correctly sam-ple the encounter probabilities, by satisfying the constraint P (cid:28) P max = . P < P max by reducing the timestep ∆ t when necessary. For the case P > P max , we decrease thesucceeding time step by letting λ → . ( P max / P ) λ . During thesubsequent step, if λ < λ o , where λ o = . P < P max , we allow the time step to increase slowly by setting λ → . λ . Once λ > λ o and the probability is satisfactorily small,which often occurs once the binary migrates out of the problematicdense region, we reset the time step factor to λ = λ o . In the case of an encounter, the relative probabilities described insection 2.3.1 determine the mass and velocity of the third object.We take this sampled velocity v to be the velocity of the thirdbody at infinity and calculate the relative velocity at infinity forthe encounter from v ∞ = | v − v | = (cid:113) v + v − vv cos χ . Given v and the sampled m and v , the relative velocity at infinity isdetermined up to the cos χ term, which for an isotropic King modeldistribution function can be sampled from an analytic expression for χ ∈ [ , π ] as in Sigurdsson & Phinney (1995). With the mass of thethird body and the relative velocity known, the maximum impactparameter is obtained from the cross section for the encounter π p = σ ( v , v ) = π p (cid:18) + G ( m b + m ) p v ∞ (cid:19) , (24)with p defined in Equation 20. The actual impact parameter forthe encounter is sampled from a uniform distribution in the areaspanned by the maximum impact parameter π p . The anglesthat comprise the remaining free variables necessary to specify theinitial conditions are the projected true anomaly f of the binary at MNRAS000 , 1–27 (2017)
Giesler, Clausen, & Ott the time that the incoming third body reaches pericenter, two angles θ and φ specifying the initial location of the third body with respectto the binary center-of-mass, and the impact orientation ψ , whichspecifies the angle of the impact parameter in a plane transverseto the incoming velocity of the third body. Theses four angles aresampled in a manner consistent with Hut & Bahcall (1983). Withthe initial conditions specified, the explicit integration is performedwith a modified scheme based on Sigurdsson & Phinney (1993).We modify the original method of a fixed initial distance ofthe third star, at R in = a , to one of variable distance to improveefficiency and to prevent the case of long three-body interactionsthat can exceed the cluster time step. The addition of massive BHsintroduces the possibility for wide binaries with orbital separationsmuch greater than those for which the previous method was suited tohandle. With a fixed choice for the distance of the third star from thebinary, interactions such as distant flybys, which are the quickest toresolve computationally and have little impact on the binary, oftentake a time that exceeds the cluster evolution time step and leads tothe possibility of missing other probable encounters.To represent the three-body system as an isolated one, and toreduce excessive time spent integrating long approaches, we requirethat R in (cid:54) R max ( n ) , where R max ( n ) = ( π n / ) − / is the ‘interpar-ticle’ distance and is a function of the local density n ( r ) . Once R in is specified, we determine the relative velocity v in at R in based onthe relative velocity at infinity. With these two quantities specified,we approximate the time for a flyby as δ t = R in / v in . For the casein which δ t > ∆ t , we let R in → ( ∆ t / δ t ) R in , calculate v in at thenew initial distance and recompute the new estimated time. We re-peat this procedure until the estimated time is roughly the sameas the cluster time step, 0 . < δ t / ∆ t < .
1. One important caveatis that this could lead to placing the third object too close to thebinary, spoiling the assumption of an object at infinity approach-ing a well defined binary. To address this issue, we maintain oneextra condition on the initial distance specification, a considerationfor which we are willing to forgo our time step restrictions: that ( a / R in ) (cid:54) . δ T max to be an arbitrarily smallfraction (cid:15) = . × − of the binary period T b , i.e. δ T max = (cid:15) T b .At the end of each integration step, we update the time step to δ T = (cid:15) ( r min / v max ) , where r min is the minimum separation betweenany pair of the three objects and v max is the largest velocity of thethree bodies. This sets the time step to the maximum allowable valuein consideration of the need to resolve the dynamics of the threeobjects or any potentially bound pair. In some instances, a resonancecan form a temporarily bound triple system, causing the integrator toreach the maximum number of steps N max = × or to exceed thearbitrarily specified maximum allowable time of 5 ∆ t . Under theserare circumstances, we reinitialize the system with newly sampledinitial angles and restart the integration. In addition to the occasionallong-lasting semi-stable triples that form, there are also instanceswhen a binary makes its way to the core where the average timescalenecessary to resolve the three body encounters begins to approachthe timescale for the evolution of the binary in the cluster. Sincewe calculate three-body encounters decoupled from the binary’sevolution in the cluster, we are forced to terminate the run in suchcases. As the cluster timescale is inversely proportional to the clusterdensity, this situation is most likely to occur in the densest clusters.As a result of this timescale termination criterion, although a similarnumber of realizations are performed for each cluster, the highestdensity clusters have noticeably fewer runs than the lower density clusters, as is observable in the rightmost column of Table 2. Forstandard encounters, which are often much shorter than the clustertime step, we periodically check whether the interaction has resolved— according to the criteria discussed in the following section —and in the case that a new binary has formed, even temporarily, weupdate δ T max with the period of this new binary. We first identify a potential binary among the triple system com-posed of the original binary, m and m , and the third mass m , byselecting the pair with the largest gravitational binding energy. Werefer to the masses in the potential binary as ¯ m and ¯ m , which mayno longer correspond to the original binary composed of m and m . The remaining object, which is not part of the potential binary,is labeled ¯ m which is distinct from m . All unbarred variables rep-resent the initial configuration where the third object is incoming,while barred variables refer to the system where a binary has beenidentified and the encounter is nearly resolved. The encounter canbe resolved in three ways: (I) there is a well defined bound binarysystem with the third object unbound and moving off to infinity, (II)a merger has occurred or (III) the system is completely ionized.For case (I), we terminate the integration once the followingcriteria are all satisfied: (i) the third body has achieved the minimumrequired separation from the binary, | ¯ r −( ¯ m ¯ r + ¯ m ¯ r )/( ¯ m + ¯ m )| > max { R max ( n ) , . R in } , (ii) the eccentricity ¯ e of ¯ m and ¯ m is lessthan unity, (iii) ¯ m and ¯ m are bound, specifically ¯ E b <
0, and (iv)¯ m is unbound, i.e. ¯ E >
0. Here, ¯ E b is the total energy of the finalbinary and ¯ E is the total energy of the third body. In addition tothe above requirements, to determine the final state of the ‘isolated’binary, we continue the integration until the total potential energybetween ¯ m and each mass in the binary is a fraction of the totalenergy of the system E , specifically G ¯ m ¯ m | ¯ r − ¯ r | + G ¯ m ¯ m | ¯ r − ¯ r | > . E . (25)In case (II), two of the bodies merge and the third body iseither unbound or forms a new binary with the merger product. Thecriteria for mergers is based on the distance of nearest approach d between two bodies during the three body encounter. In the case ofa potential merger between two BHs, the merger criterion is d (cid:54) R + R . For the remaining merger situations, the criterion remains d = f tid ( R + R ) , as adopted from Sigurdsson & Phinney (1995),using the same value for f tid as introduced in section 2.2.5. Whena merger occurs between the BH companion and the third body,if the merger product remains bound to the BH, this dynamicallyformed binary becomes our new ‘test binary’, which we continueto follow and evolve within the cluster. Similarly, if the BH mergeswith a third body and we still have a bound binary system, weagain continue to follow this binary. However, if the BH becomesunbound by merging with another body or becomes unbound froma merger product, we handle the newly single BH as described in thesubsequent section. In each of these cases, the position of the newbinary, or single BH, is updated by continuing along the originalbinary trajectory and the velocity is updated by converting from thethree-body center-of-mass frame, where the three-body integrationis performed, back to the cluster frame.The result of the encounter can also end in complete ionization,case (III). Ionization occurs in the case of ill-defined binaries thatwill inevitably be unbound if, given that all previous criteria aresatisfied, either (a) the eccentricity of ¯ m and ¯ m satisfies 1 − ¯ e < × − or (b) | ¯ r − ¯ r | > R max ( n ) is satisfied. Additionally, ionization MNRAS , 1–27 (2017)
H-LMXBs from BH retaining GCs Name M [ M (cid:12) ] σ [ cm s − ] ρ L [ L (cid:12) pc − ] c N BH N runs Pal 13 5.12 × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × Table 2.
Summary of simulations. Listed are the 15 GCs modeled for evolution along with the total cluster mass M c , squared velocity dispersion σ , theluminous core density ρ L and concentration c . The clusters are ordered by total mass. There are 39 independent models after taking into account the numberof BHs retained by the cluster. Medium to high mass clusters can accommodate large BH populations without disrupting the listed structural parameters. Thesize of the BH population in lower-mass clusters is either (1) limited in number by the IMF or (2) by the ability of the cluster to maintain the model structuralparameters in their presence; in these cases, the cluster is not used for evolutions and is omitted from the table. In the final column we list the total number ofevolutions performed for each case. occurs if ¯ m i ¯ v i > (cid:0) ¯ m i ¯ m j /| ¯ r i − ¯ r j | + ¯ m i ¯ m k /| ¯ r i − ¯ r k | (cid:1) is true for allmasses at any time, with i (cid:44) j (cid:44) k taking on values { , , } . Thislast criterion is a straightforward definition for a totally unboundtriple. In addition to these choices for ionization during three-bodyencounters, there is one other instance in which the binary can bedissociated. For very wide binaries, the encounters are dominatedby repeated grazing encounters with low mass stars, which tend tofurther widen the orbital separation. As a result, strong interactionsbecome less likely and the binary will inevitably be dissociatedby the increasing occurrence of these slowly ionizing encounters.For this reason, we use the encounter rate to define a maximumsemi-major axis of dynamically formed binaries as a max ( Γ ) = (cid:18) Gm b ( π Γ ) (cid:19) / , (26)which is equivalent to requiring a minimum of three orbits be-tween encounters. Here, the total encounter rate Γ = (cid:205) α Γ ( r , v , a ) is a sum over the rate associated with each mass group de-fined by Equation 22). The final criterion for ionization is then a > min { a max ( Γ ) , R max ( n )} . As described in the previous section, a BH can become single dueto three-body dynamics such as exchange, merger, or through thedismantling of a binary that exceeds our large a or large e criteria.In the case of a single BH, we allow for the solitary BH to form anew binary by interacting with existing binaries within the cluster.In order to accomplish this, we need to know the probabilityfor the following encounter, ( m , m ) + m BH → ( m BH , m ) + m , (27)in which the BH exchanges with m into a binary originally com-posed of masses m and m . We also consider the possibility that MNRAS000
Summary of simulations. Listed are the 15 GCs modeled for evolution along with the total cluster mass M c , squared velocity dispersion σ , theluminous core density ρ L and concentration c . The clusters are ordered by total mass. There are 39 independent models after taking into account the numberof BHs retained by the cluster. Medium to high mass clusters can accommodate large BH populations without disrupting the listed structural parameters. Thesize of the BH population in lower-mass clusters is either (1) limited in number by the IMF or (2) by the ability of the cluster to maintain the model structuralparameters in their presence; in these cases, the cluster is not used for evolutions and is omitted from the table. In the final column we list the total number ofevolutions performed for each case. occurs if ¯ m i ¯ v i > (cid:0) ¯ m i ¯ m j /| ¯ r i − ¯ r j | + ¯ m i ¯ m k /| ¯ r i − ¯ r k | (cid:1) is true for allmasses at any time, with i (cid:44) j (cid:44) k taking on values { , , } . Thislast criterion is a straightforward definition for a totally unboundtriple. In addition to these choices for ionization during three-bodyencounters, there is one other instance in which the binary can bedissociated. For very wide binaries, the encounters are dominatedby repeated grazing encounters with low mass stars, which tend tofurther widen the orbital separation. As a result, strong interactionsbecome less likely and the binary will inevitably be dissociatedby the increasing occurrence of these slowly ionizing encounters.For this reason, we use the encounter rate to define a maximumsemi-major axis of dynamically formed binaries as a max ( Γ ) = (cid:18) Gm b ( π Γ ) (cid:19) / , (26)which is equivalent to requiring a minimum of three orbits be-tween encounters. Here, the total encounter rate Γ = (cid:205) α Γ ( r , v , a ) is a sum over the rate associated with each mass group de-fined by Equation 22). The final criterion for ionization is then a > min { a max ( Γ ) , R max ( n )} . As described in the previous section, a BH can become single dueto three-body dynamics such as exchange, merger, or through thedismantling of a binary that exceeds our large a or large e criteria.In the case of a single BH, we allow for the solitary BH to form anew binary by interacting with existing binaries within the cluster.In order to accomplish this, we need to know the probabilityfor the following encounter, ( m , m ) + m BH → ( m BH , m ) + m , (27)in which the BH exchanges with m into a binary originally com-posed of masses m and m . We also consider the possibility that MNRAS000 , 1–27 (2017) Giesler, Clausen, & Ott m BH and m undergo an exchange, which contributes to the totalprobability that the BH will exchange into the binary. However, forconciseness in deriving the probability of exchange, we will focusspecifically on the encounter described by Equation 27, later addingthe contribution from the reaction where the subscripts are inter-changed. Unfortunately, we can no longer compute the probabilityfor encounter as in section 2.3.1, since we do not possess a distri-bution function for binaries. However, by considering the reversereaction of Equation 27, given by ( m BH , m ) + m → ( m , m ) + m BH , (28)and relating this to the one of interest, we can obtain the encounterprobability for the BH to exchange into an existing binary in thesame way that we compute encounters for a binary composed of aBH and a companion.We use the seminumerical fit of Heggie et al. (1996),¯ σ , = (cid:18) M M (cid:19) / (cid:18) m M (cid:19) / (cid:18) M M (cid:19) / (cid:18) M M (cid:19) g ( , , ) , (29)as the dimensionless cross section for a generically labeled singlemass m to exchange into a binary of masses m and m to forma new binary composed of m and m , with m being ejected. Inthis notation, uppercase masses represent the sum of the mass sub-scripts, i.e. M ij = m i + m j . The coefficient g ( , , ) is a numericalfitting factor designed to improve the analytically derived fit. Thisdimensionless cross section ¯ σ , is related to the dimensionful crosssection for exchange Σ , through¯ σ , = | v , − v | π GM a , Σ , . (30)The existing binaries that the BH is likely to encounter, which haveremained intact in the cluster over long timescales, can be considered‘hard’. These ‘hard’ binaries are characterized by having a bindingenergy U bin that exceeds the average energy of the other stars inthe cluster | U bin | > ¯ m ¯ σ and this is what allows them to stayintact over such long timescales. In this case, we approximate thetotal encounter cross section by the dominant gravitational focusingterm in Equation 21, explicitly: σ , (cid:39) π GM a , | v , − v | . (31)Finally, relating Equation 30 and Equation 31 allows us to expressthe cross section for exchange in terms of the total encounter crosssection σ , through Σ , = (cid:0) ¯ σ , / (cid:1) σ , . (32)Evidently, the dimensionless cross section for exchange is related tothe fractional probability that the total encounter ends in the specificexchange we previously described. Considering Equation 31 andassuming the relative velocities are similar for the forward andreverse reactions, we can relate the forward and backward totalcross sections through σ , = ( a , a , ) σ , . Since the energy givento the binary is comparable to the energy required to destroy it, m m / a , ∼ m m / a , , we can recast the relation in terms of themasses alone: σ , = (cid:18) m m (cid:19) σ , . (33)The cross section for the specific exchange of m for m in terms ofthe total encounter cross section of the original binary is found bysubstituting Equation 33 into Equation 32, yielding Σ , = (cid:18) ¯ σ , m m (cid:19) σ , . (34) By writing the exchange probability in terms of the post-exchangebinary, we can now utilize the same procedure described in sec-tion 2.3.1. In this formalism, m represents the BH and we returnto referring to this body as m BH , while m goes to m α , a variablecompanion used for computing the relative probabilities for eachmass group α . First we select a companion object m for the BHon the left-hand side of Equation 28. We obtain m by samplingfrom the local number density and determine a and e for the binaryas in section 2.2.5. The probability of the encounter described byEquation 27, where the BH exchanges places with m α in a binarycomposed of m and m α is then, P α, = ∆ t ∫ (cid:18) ¯ σ α, m α m BH (cid:19) σ , BH ( v , v α )| v − v α | f α ( v α ) d v α . (35)The usefulness of the manipulations in this section is most clearlyseen by writing this in terms of Equation 23: P α, = (cid:18) ¯ σ α, m α m BH (cid:19) P α , (36)which in practice makes computing the exchange probabilities aseasy as rescaling our standard encounter computations by the par-enthetical factor. Since we also allow for the BH to exchange with m , we also consider the probability P ,α = (cid:18) ¯ σ ,α m m BH (cid:19) P .We apply one final rescaling to account for the density ofbinaries that are of type m and m α . We assume that the fractionof objects that are binaries f b is constant throughout the clusterwith the value specified by Equation 16. The density of binariesis then n b ( r ) = ( f b + f b ) n ( r ) , which is derived from Equation 15.Additionally, we also assume that the fraction of binaries of a giventype is constant at all cluster radii, n ij ( r ) = f i / j n b ( r ) . Here, f i / j represents the fraction of binaries that have a star of type i and astar of type j , e.g. f NS / MS is the fraction of all binaries that arecomposed of an NS and an MS star. For binaries composed of onlyMS or WD we use values of f MS / MS = . f MS / WD = .
44, and f WD / WD = .
32 (Fregeau et al. 2009). The remaining one percentof binaries contain at least one BH or NS, for which we computethe binary fraction through f i / j = . ( N i N )( N j N BH + NS ) , where i canbe any object type, j is limited to BH or NS, N is the total numberof objects in the cluster, and N BH + NS is the total number of BHsand NSs.The final total probability for the BH to exchange into a binary,given the sampled mass m , is then P exch ( r ) = (cid:213) α n α ( r ) (cid:18) P α, n α ( r ) + P ,α n ( r ) (cid:19) . (37)Here we divide out the respective local density picked up in theintegration of the distribution function in order to enforce our as-sumption of a uniform binary fraction. If an exchange is determinedto occur based on this total probability, we select a specific binaryfor the encounter based on the relative probabilities of exchange foreach mass group m α . With a binary in hand, we initiate our three-body system, which is run until we get the proper outcome dictatedby the encounter cross section — i.e. that m BH exchanges with theappropriate mass in the binary. We present 698,486 realizations from 15 GC models with totalmasses in the range of 5 . × – 2 . × M (cid:12) , velocity disper-sions covering 9 × – 1 . × cm s − , core densities of 1 . MNRAS , 1–27 (2017)
H-LMXBs from BH retaining GCs – 1 . × pc − , and concentrations spanning 0 .
66 – 2 .
07. Thesimulations are summarized in Table 2, which includes the cata-log name for the modeled cluster, total mass, velocity dispersionsquared, central luminosity density, concentration, the number ofretained BHs in the model, and the total number of completed runs.The simulations are run for t = years or until the single/binaryis ejected from the cluster, when r > r t . In our framework, a GC’s structure is determined by four parame-ters: the total cluster mass M c , the core velocity dispersion σ , thecore luminosity density ρ L , and the concentration c . McLaughlin(2000) finds that GCs described by single-mass isotropic King mod-els are fully defined by four independent physical parameters: themass-to-light ratio Υ v , , total binding energy E b , central concen-tration c , and total luminosity L . Furthermore, McLaughlin (2000)shows that Milky Way GCs lie in a ‘fundamental plane’ and thuscan be fully described by just two independent parameters, c and L . A face-on view of the fundamental plane is defined by the axes (cid:15) = .
05 log E ∗ b + log L and (cid:15) = c . The apparent dependenceon the third quantity log E ∗ b is due to a rotation in the larger threedimensional space in order to remove projection effects. However,this is reconciled by showing that this third parameter, E ∗ b , is fullydescribed by the luminosity, such that E ∗ b ( L ) (McLaughlin 2000).With the space of physical clusters reduced to the fundamental plane,we determine a representative group of 15 Milky Way clusters bysampling from the two-dimensional distribution. A face-on viewof the fundamental plane is given in Figure 3, which includes allGCs from the Harris catalog (Harris 1996, 2010 edition) for whichobserved concentrations are available. We omit clusters identifiedin the catalog as core-collapsed, since these are not generally welldescribed by King models. This includes those with c = .
5, anarbitrary value assigned to clusters in the catalog with central den-sity cusps indicative of core collapse. There are 125 Milky WayGCs remaining after core-collapse pruning; of these, 15 GCs arechosen as representative models, in an attempt to properly coverthe fundamental parameter space. The 15 Milky Way GC modelsrepresentative of the 125 Milky Was GCs are described in Table 2and represented by stars in Figure 3 to visualize our coverage of thefundamental parameter space.As stated in section 2.1, our input parameters for specifying thestructure of a cluster are the core velocity dispersion ¯ σ , the centraldensity n o , and the King parameter W o . The mean core velocitydispersion ¯ σ is chosen to be the observed value listed in the Harriscatalog. The core number density n o is adjusted until the centralluminosity density ρ L is consistent with observation. Finally, theKing parameter W o , which sets the depth of the potential, is varieduntil the cluster has the desired total mass M c and concentration c .Once we have a model for a given GC, we add BHs by increasing thefraction of retained BHs f r BH , where a value of unity correspondsto retention of all BHs produced according to the IMF. For a givennumber of BHs in the cluster, we use the parameter f s in Equation 18to adjust the BH velocity dispersion such that the overall structureof the cluster is unaffected by the presence of a significant numberof BHs. However, we find that there is a limit to the number of BHseach cluster can harbor. For the lowest-mass clusters, such as Pal 13,setting the retention factor to unity, f r BH =
1, in order to maximizethe number of BHs retained by the cluster produces a peak numberof ∼
20 BHs. In this case, the number of BHs retained by the clusteris inherently limited by its structure. More generally, for lower-massclusters that allow for more BHs, the large number of BHs can
95 100 105 110 e e modeled GCunmodeled GC Figure 3.
The distribution of non core-collapsed Milky Way GCs in a face-on view of the fundamental plane. The color of each unmodeled GC (markedby circles) indicates the corresponding modeled GC (marked by stars) thatserves as its proxy for determining the properties of the ejected binaries.The plane is defined by (cid:15) = .
05 log E ∗ b + log L and (cid:15) = c , withthe dashed line corresponding to the fit (cid:15) = − . + . (cid:15) . Here c is theconcentration, L is the total luminosity, and E ∗ b is an additional parameterrelated to L (see section 3.2 for additional details). become problematic as they become a more significant part of thetotal mass of the cluster. As the fraction of the total mass in BHsincreases, the BHs begin to affect the structural parameters such thatno set of initial parameters exists that satisfy the observed structureof the GC. We find that for many of the lower-mass clusters we areonly able to simulate populations of 20 or 200 BHs (cf. Table 2). The GC evolution models, described in detail in section 2, computethe properties of the BH binaries at the moment they are ejected froma GC. Determining the present day properties of potentially observ-able, ejected BH binaries requires further modeling that tracks boththe evolution of ejected binaries in the Milky Way potential and theinternal evolution of each binary. In this section, we describe MonteCarlo models for the subsequent evolution of the ejected binariesthat are seeded with results from our GC models.
We first build a sample of GCs to include in our galactic evolutionsimulations. The orbit of a cluster is specified by its location onthe sky (right ascension and declination), distance from the Sun D (cid:12) , radial velocity v r , and proper motion µ α and µ δ . Of the 125non core-collapsed GCs in the Harris catalog (Harris 1996, 2010edition), we are able to find literature values for the orbital pa-rameters of 106 of these clusters in the catalogs of Moreno et al.(2014) and Kharchenko et al. (2013). For clusters appearing in bothcatalogs we use the values given in Moreno et al. (2014).To begin each realization in our Monte Carlo ensemble, weinitialize the GC orbits by sampling the uncertainty in their currentpositions and velocities. We assume normally distributed errors and MNRAS000
We first build a sample of GCs to include in our galactic evolutionsimulations. The orbit of a cluster is specified by its location onthe sky (right ascension and declination), distance from the Sun D (cid:12) , radial velocity v r , and proper motion µ α and µ δ . Of the 125non core-collapsed GCs in the Harris catalog (Harris 1996, 2010edition), we are able to find literature values for the orbital pa-rameters of 106 of these clusters in the catalogs of Moreno et al.(2014) and Kharchenko et al. (2013). For clusters appearing in bothcatalogs we use the values given in Moreno et al. (2014).To begin each realization in our Monte Carlo ensemble, weinitialize the GC orbits by sampling the uncertainty in their currentpositions and velocities. We assume normally distributed errors and MNRAS000 , 1–27 (2017) Giesler, Clausen, & Ott use the quoted uncertainties in v r , µ α , and µ δ . Following Krauss &Chaboyer (2003), we assume a 6 per cent error in D (cid:12) . After the orbitis specified, we integrate it 10 Gyr backward in time, correspondingto the duration of our GC dynamical simulations.The orbits of the GCs, and the ejected binaries, are integratedusing the python galactic dynamics library galpy (Bovy 2015).We model the Milky Way gravitational potential using the builtin MWPotential2014 . The potential includes contributions fromthe galactic bulge, disk, and halo, which have been fit to observa-tional data to provide a realistic model of the Milky Way potential.The physical scale of the potential is set using the distance fromthe center of the Galaxy to the Sun and the circular velocity of theSun, which we set to 8 kpc and 220 km s − , respectively. For allcalculations, we use the dopr54_c integrator, a fast implementationof a high order Dormand-Prince method included with galpy .Now that we have calculated the positions and velocities of theMilky Way GCs during the past 10 Gyr, the next step is to determinethe properties of any potential BH-LMXBs ejected by these clus-ters. Since our dynamical simulations only include a subset of thegalactic GCs, we use the results from the 15 GCs simulated in Ta-ble 2 as proxies for the ejected binary populations produced by theremaining 110 clusters in our galactic evolution models. For eachof the unmodeled clusters, a proxy cluster is selected by finding thenearest simulated cluster in the fundamental plane (see section 3.1).Specifically, we find min (cid:104) ( (cid:15) (cid:48) , i − (cid:15) (cid:48) , j ) + ( (cid:15) (cid:48) , i − (cid:15) (cid:48) , j ) (cid:105) , where the i index runs over all 106 clusters in the galactic evolution models,the j index runs over the 15 clusters included in our GC dynamicsmodels, and the primes denote the normalized versions of (cid:15) and (cid:15) restricted to the range [ , ] . Figure 3 shows the proxy cluster chosenfor each GC, by assigning the same color marker to each GC as thecolor of the proxy cluster used, which are marked by colored stars.To ensure the robustness of this method for choosing a proxy clus-ter, we assign a proxy by two additional methods. One secondarymethod is to assign the proxy cluster based on the minimum dis-tance in the fundamental plane using the unnormalized axes (cid:15) and (cid:15) . The second alternative is by identifying the most similar clusterusing the structural parameters M c , σ , and ρ L weighted accordingto the strengths of the correlations between these parameters and theejected binary populations, which are explored in 4.1. Selecting theproxy cluster by any of these three methods gives similar results inour galactic evolution models. In fact, all three methods will selectthe same proxy cluster for all but ∼
15 of the 110 unmodeled GCsin our study. In what follows, we discuss models that use the scaleddistance in the fundamental plane to assign the proxy cluster.
The output of our GC dynamical simulations describes the proper-ties of the BH-binaries ejected from GCs. To model the present daypopulation of BH-LMXBs that are ejected from GCs, we use as in-puts for our galactic evolution models: the ejection time t ej , ejectionvelocity v ej , and the properties of the binary, the semi-major axis a , eccentricity e , the mass of the BH primary m , and the mass ofthe companion m . This is accomplished by constructing empiri-cal cumulative distribution functions (CDFs) of these quantities foreach of the 37 sets of parameters listed in Table 2, and then sam-pling these distributions in our Monte Carlo models. We assumethat the ejection time, ejection velocity, and binary properties areindependent and sample the marginal distributions of each. http://jobovy.github.io/galpy/ In the GC dynamical models, a , e , t ej , and v ej are treated ascontinuous variables. As such, we are able to sample the CDFs forthese quantities directly. We fit cubic splines to the empirical CDFsand invert the distributions by interpolation. The GC dynamicalmodels treat m and m as discrete quantities, which fall into themass bins shown in Table 1. In our galactic evolution models, how-ever, we want to consider continuous masses. To accomplish this,we first determine an object’s mass bin by sampling the discreteCDF output by the dynamical simulations. Next, we sample themass distribution within that bin using the evolved mass functiondescribed in section 2.2.1. Using these CDFs, we are able to gener-ate sample populations of the BH-binaries ejected by the 106 GCsin our galactic evolution simulations.During each realization, for each cluster, we first determinethe number of binaries that the cluster will eject during the 10 Gyrsimulation by sampling a Poisson distribution with rate parameter (cid:104) N ej (cid:105) (third column of Table 3). Once we have determined thenumber N bin of ejected binaries, we draw N bin samples from the a , e , m , m , t ej , and v ej distributions.Since the internal evolution of a binary is independent of its or-bit in the Galaxy, we separately compute the full internal evolution ofthe binary using the rapid binary population synthesis code BSE de-scribed in Hurley et al. (2002) with the updates described in Clausenet al. (2012) and Lamberts et al. (2016).
BSE combines interpolatedstellar evolution models with recipes for mass-transfer and otherbinary evolution processes to enable rapid modeling of a binarysystem’s lifetime. Binary population synthesis calculations employparameterized models to describe poorly understood processes inbinary evolution. In our
BSE runs, we assume that stable masstransfer is conservative. Additionally, we use a common-envelopeefficiency parameter of 1 . a , e , m , m as the initial conditions fora BSE run. Furthermore, we set the companion star’s metallicity tothat of its parent GC and its age to t ej . The latter has little effectbecause most of the ejected stars have lifetimes that exceed 10 Gyr.The binary is evolved for t evol =
10 Gyr − t ej , i.e., to the presentday. Systems are discarded if the companion star is not overflowingits Roche-lobe and transferring mass to the BH at the end of thesimulation. For each mass transferring binary, we determine theposition r GC and velocity v GC of its parent GC at t ej . We initializean orbit for the ejected binary at r GC and v GC + v ej , assumingthat the binaries are ejected isotropically. With the initial conditionsdetermined, we then evolve these binaries using galpy to determinetheir positions at the present day.Our galactic evolution models consider three BH-retention sce-narios. In the first, we assume that most BHs are ejected and use theresults from our GC dynamics models with N BH = M . We referto this set of models as MIN. In the second case, referred to as 200,we assume moderate BH retention, using the results from our GCdynamics models with N BH = N BH = N BH , we use theresults from the model with nearest N BH simulated for that samecluster. We compute 10 realizations for the MIN and 200 cases and5 × realizations for the MAX case. MNRAS , 1–27 (2017)
H-LMXBs from BH retaining GCs − − − − − − − log ( N BH c M c [ M − (cid:12) ]) − − − l o g ( h N e j i ) BH-NCBH-NSBH-BH
Figure 4.
Expected number of binary ejections (cid:104) N ej (cid:105) as a function of thenumber of retained BHs N BH , concentration c , and total cluster mass M c .The number of binaries ejected over the life of the cluster is well describedby the two characteristic variables of the fundamental plane, c and M c ,along with the number of BHs retained by the cluster. Our simulations of binary-single star interactions in GCs provideus with statistical properties of the ejected BH binaries they pro-duce including ejection time t ej , ejection velocity v ej , the orbitalproperties a and e , and the component masses m and m . Com-bining these results with the methods described in section 3.2, weobtain predictions for the distribution and properties of the galac-tic population of BH-LMXBs produced by GCs. Additionally, thesimulations allow us to explore merger events involving BHs suchas gravitational radiation driven mergers, both in the cluster andpost-ejection, as well as those mergers that occur during three-bodyencounters. We describe these results in detail below. We find that the number of ejected binaries and the properties ofthese binaries are strongly affected by the GC structure and thenumber of retained BHs. In Table 3, we list the expected number ofejected BH binaries over the life of each cluster, listed in order ofincreasing mass, including the exact number of BHs in each cluster.The ejected BH-binary expectation value is well described by thenumber of retained BHs N BH and the two characteristic variablesthat define the fundamental plane of GCs (see Figure 3), namelythe total cluster mass M c and the concentration c . In Figure 4, weplot the expected number of ejected BH binaries as a function ofthe three characteristic variables: N BH , M c , and c .The most important structural variable that impacts the ejectedbinary properties is the cluster mass. The total cluster mass enforcesa minimum energy needed to escape, which the binary must gainthrough repeated encounters. In order for a binary to escape fromthe cluster, it must acquire a recoil velocity from a final three-bodyencounter high enough to climb out of the cluster gravitationalpotential. In Figure 5, we show the distribution of the ejected binaryvelocities as a function of cluster mass, where the influence of themass of the cluster on the ejection velocity is apparent. The expected log ( M c [ M (cid:12) ]) − − l o g ( v e j [ k m s − ] ) d P / d l o g v e j Figure 5.
The distributions of ejection velocities v ej as a function of the totalcluster mass M c for the ejected binaries. Each vertical bar represents thedistribution of v ej for the corresponding mass M c and is normalized such thatthe integral over log v ej in each mass bin yields unity. The binary velocityfluctuates due to random encounters with other stars in the cluster untilthe binary acquires a high enough recoil velocity to exceed the minimumejection velocity, which is determined by the cluster mass. The increase inthe necessary velocity for escape is apparent in the increasing mean valueof each v ej distribution. number of ejections is then higher for lower-mass clusters due to thelower escape velocities associated with these clusters, as is visiblein Figure 4. To decouple this statement from the additional variablesin Figure 4, it can also be observed in Table 3 (which is orderedby increasing mass) that for a fixed number of retained BHs, theexpected number of ejections scales with the cluster mass.The mechanism through which the binary converts bindingenergy to kinetic energy is easiest to understand in the three-bodycenter of mass frame, where we perform our integration for encoun-ters. After an encounter, the final relative velocity at infinity is givenby¯ v ∞ = m ( m + m ) ¯ m ( ¯ m + ¯ m ) v ∞ + M ¯ m ( ¯ m + ¯ m ) ( U bin − ¯ U bin ) , (38)where U bin = − Gm m a is the binding energy of the binary and allunbarred quantities represent the initial binary before encountering m , while barred quantities represent the final binary and ¯ m is theejected mass. In the case of no exchange, and utilizing ∆ a ≡ ¯ a − a ,Equation 38 reduces to¯ v ∞ = v ∞ − M m m b (cid:18) Gm m ∆ aa (cid:19) . (39)In this frame, the binary velocity is related, through conservationof momentum, to the relative velocity simply by v b = m M v ∞ . Thechange in the kinetic energy, ∆ T ≡ ¯ T − T , of the binary is then ∆ T = − Gm m m M (cid:18) ∆ aa (cid:19) . (40)The amount by which the semi-major axis changes in an averageencounter, where the semi-major axis is reduced without exchange,is proportional to the semi-major axis, ∆ a ≈ − (cid:15) a , with (cid:15) in the range ∼ [ , . ] (Sigurdsson & Phinney 1993). Using this relation, and MNRAS000
The distributions of ejection velocities v ej as a function of the totalcluster mass M c for the ejected binaries. Each vertical bar represents thedistribution of v ej for the corresponding mass M c and is normalized such thatthe integral over log v ej in each mass bin yields unity. The binary velocityfluctuates due to random encounters with other stars in the cluster untilthe binary acquires a high enough recoil velocity to exceed the minimumejection velocity, which is determined by the cluster mass. The increase inthe necessary velocity for escape is apparent in the increasing mean valueof each v ej distribution. number of ejections is then higher for lower-mass clusters due to thelower escape velocities associated with these clusters, as is visiblein Figure 4. To decouple this statement from the additional variablesin Figure 4, it can also be observed in Table 3 (which is orderedby increasing mass) that for a fixed number of retained BHs, theexpected number of ejections scales with the cluster mass.The mechanism through which the binary converts bindingenergy to kinetic energy is easiest to understand in the three-bodycenter of mass frame, where we perform our integration for encoun-ters. After an encounter, the final relative velocity at infinity is givenby¯ v ∞ = m ( m + m ) ¯ m ( ¯ m + ¯ m ) v ∞ + M ¯ m ( ¯ m + ¯ m ) ( U bin − ¯ U bin ) , (38)where U bin = − Gm m a is the binding energy of the binary and allunbarred quantities represent the initial binary before encountering m , while barred quantities represent the final binary and ¯ m is theejected mass. In the case of no exchange, and utilizing ∆ a ≡ ¯ a − a ,Equation 38 reduces to¯ v ∞ = v ∞ − M m m b (cid:18) Gm m ∆ aa (cid:19) . (39)In this frame, the binary velocity is related, through conservationof momentum, to the relative velocity simply by v b = m M v ∞ . Thechange in the kinetic energy, ∆ T ≡ ¯ T − T , of the binary is then ∆ T = − Gm m m M (cid:18) ∆ aa (cid:19) . (40)The amount by which the semi-major axis changes in an averageencounter, where the semi-major axis is reduced without exchange,is proportional to the semi-major axis, ∆ a ≈ − (cid:15) a , with (cid:15) in the range ∼ [ , . ] (Sigurdsson & Phinney 1993). Using this relation, and MNRAS000 , 1–27 (2017) Giesler, Clausen, & Ott
Name N BH BH-NC BH-NS BH-BHPal 13 19.64 3.14 7 . × − . × − NGC 6838 20.61 6 . × − . × − . × − . × . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − NGC 6121 20.70 3 . × − . × − . × − . × − . × . × − . × − . × − . × . × NGC 6093 19.85 1 . × − . × − . × − . × − . × . × NGC 5286 12.29 6 . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × . × NGC 6205 20.10 6 . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − Table 3.
Expected number of binary ejections. For each cluster and number of retained BHs, we list the exact number of BHs in the cluster along with theexpected number of ejections over the cluster lifetime for three binary types: BH-NC, BH-NS, and BH-BH. The clusters follow the same order as Table 2, sortedaccording to increasing total cluster mass. The values of N BH are non-integer values as a consequence of modeling the population with a smooth distributionfunction. assuming a binary with constant m and m , Equation 40 reducesto ∆ T ∝ m M (cid:15) a , (41)yielding a simple relation that describes the gain in kinetic energy interms of the constant fractional change in the semi-major axis (cid:15) andthe ratio of the third body to the total mass of the three-body system.Additionally, Equation 41 shows that this change in kinetic energybecomes more efficient as the semi-major axis decreases, convert-ing more energy from binding to kinetic after each encounter thatshrinks the binary’s orbit. After repeated interactions, the increasein velocity due to the decrease in a becomes more substantial andthe binary can eventually reach the necessary velocity to escape.We can directly relate the necessary gain in kinetic energy to thechange in binding energy ∆ U = ¯ U bin − U bin , by simply rearrangingEquation 38 and assuming no exchange of masses, which yields ∆ T = − m M ∆ U . (42) In the process of the binary increasing its kinetic energy, the bindingenergy becomes more negative. Since the higher-mass clusters tendto hold on to the binaries longer, this strict minimum kinetic energyfor ejection is manifest in the more negative-valued binding energyof the binaries it ejects. It follows from this, that on average, thesemi-major axes of the binaries ejected from more massive clusterstend to be smaller. This is confirmed by Figure 6, which depicts thedistribution of orbital separations as a function of cluster mass.In addition to the increase in the expected number of ejectedbinaries in lower-mass clusters, the total number of expected ejec-tions also increases with an increase in the number of BHs. Whilethe number of ejections is expected to increase with the number ofBHs, interestingly, the fraction of ejected binaries that are BH-NCalso grows with the number of BHs (see Figure 4 and Table 3). Thisbehavior can be attributed to the fact that the BHs are not in energyequipartition with the rest of the cluster. Adding more BHs withoutaffecting the distribution of the luminous cluster members requiresthat the BHs are spread out farther from the core, where they have MNRAS , 1–27 (2017)
H-LMXBs from BH retaining GCs log ( M c [ M (cid:12) ]) − − l o g ( a [ AU ] ) d P / d l o g a Figure 6.
The distributions of the semi-major axes at ejection a as afunction of the total cluster mass M c for the ejected binaries. Each verticalbar represents the distribution of a for the corresponding mass M c andis normalized such that the integral over log a in each mass bin yieldsunity. High mass clusters require a high velocity for escape, which a binarymust acquire through three-body interactions in order to be ejected. Theenergy needed to escape is more easily gained once the orbital separationhas decreased sufficiently (see Equation 41). As a consequence, the meanvalue of a at ejection shifts to smaller separation with increasing clustermass M c . traditionally been expected to reside. Accordingly, the mean den-sity of BHs goes down, and they are less likely to interact with eachother. However, because they are well mixed with the stars at largerradii, the number of BH-NC binaries that form in three-body ex-changes grows. Additionally, since these binaries form farther fromthe core, they also have the benefit of a shallower potential to climbout of.Besides influencing the number of ejected binaries, the numberof retained BHs also affects the distribution of the semi-major axesof the ejected binaries. In Figure 7, we show the distribution ofsemi-major axes for the ejected BH-NC binaries in our cluster modelfor NGC 5694 for the three different choices of BHs retained. Wechoose this cluster since it is representative of the effect that thenumber of retained BHs has on the population of ejected BH-NCbinaries. Figure 7 displays an increase in the width of the distributionof semi-major axes for larger populations of BHs. This is againrelated to the necessary spreading of the BHs as we increase thenumber of BHs harbored by the cluster. Therefore, the BH-NCbinaries that form outside of the core, where the escape velocitydrops rapidly as a function of radius, can be ejected while theirbinding energies are of comparably lower magnitudes. Althoughthe more widely separated binaries are less likely to become mass-transferring systems, the simulations with large BH numbers tendto have much higher ejection rates. The higher ejection rates stillproduce enough tight binaries in the tail of distribution to outnumberthose produced with fewer BHs present.The remaining structural property of GCs that has a clear ef-fect on the population of ejected binaries is the cluster density. InFigure 8, we plot the distribution of ejection times as a functionof the luminous central density, which is related to the core den-sity as discussed in section 2.2.2. The distribution establishes thatthe cluster density has some impact on the time at which binaries − − − − − log ( a [ AU ]) d P / d l o g a N BH = N BH = N BH = Figure 7.
The probability distribution for the ejected BH-NC binary semi-major axes from NGC 5694, a representative case, with a population of 20,200, and 1000 BHs. An increase in the number of BHs requires spreadingthe BHs outside of the core, where they are more likely to form binarieswith NC objects. In the outskirts, the energy necessary to escape is muchsmaller, allowing the binary to escape before it has had sufficient time toharden. These binaries escape with comparatively low magnitude bindingenergy and wide orbital separations. log ( ρ L [ L (cid:12) pc − ]) l o g ( t e j [ y r ] ) d P / d l o g t e j Figure 8.
The distributions of time of ejection t ej as a function of the lumi-nous central density ρ L for the ejected binaries. Each vertical bar representsthe distribution of t ej for the corresponding core luminosity density ρ L andis normalized such that the integral over log t ej in each density bin yieldsunity. In higher density clusters, where encounters occur more frequently,many binaries are ejected after only a few Gyr, while in the lower densityclusters most ejections occur near the end of the 10 Gyr evolution. are ejected from their host GC. The time between binary-singleencounters can be approximated by t enc = Γ − = v m π G ( m b + ¯ m ) n o a , (43)where v m is the mean velocity of stars in the cluster, n o is its coredensity, and ¯ m is the mean mass. Combining this result with Equa- MNRAS000
The distributions of time of ejection t ej as a function of the lumi-nous central density ρ L for the ejected binaries. Each vertical bar representsthe distribution of t ej for the corresponding core luminosity density ρ L andis normalized such that the integral over log t ej in each density bin yieldsunity. In higher density clusters, where encounters occur more frequently,many binaries are ejected after only a few Gyr, while in the lower densityclusters most ejections occur near the end of the 10 Gyr evolution. are ejected from their host GC. The time between binary-singleencounters can be approximated by t enc = Γ − = v m π G ( m b + ¯ m ) n o a , (43)where v m is the mean velocity of stars in the cluster, n o is its coredensity, and ¯ m is the mean mass. Combining this result with Equa- MNRAS000 , 1–27 (2017) Giesler, Clausen, & Ott
14h 16h 18h 20h 22h 0h 2h 4h 6h 8h 10h − ◦ − ◦ − ◦ − ◦ − ◦ ◦ ◦ ◦ ◦ ◦ ◦ GCBH-LMXB − − − − − − − − d l o g P / d d e g Figure 9.
The spatial probability distribution of the simulated population of BH-LMXBs from GCs with N BH = h . m , − . ◦ , where the high density of objects explains the clustering ofBH-LMXBs and GCs. tion 40, we can obtain an approximation for the rate at which a binaryincreases its kinetic energy ∆ T / ∆ t . As encounters approximatelyoccur in increments of the encounter timescale, letting ∆ t = t enc ,we find that the rate at which the binary increases its kinetic energy, ∆ T ∆ t = (cid:18) π G m m m (cid:15) v m (cid:19) n o , (44)scales with the cluster core density. Therefore, the time it takes fora binary to acquire a high enough velocity to escape is reducedfor higher density clusters. As can be seen in Figure 8, in clustersof higher density, where encounters occur more frequently, mostBH-NC systems are ejected after only 3 Gyr of evolution whereasin lower density clusters most ejections take place near the end ofthe 10 Gyr simulation (i.e. the present day), Here we focus strictly on the population of the present-day mass-transferring systems that have successfully become BH-LMXBs.These results reflect the contribution to the BH-LMXB populationfrom the entire population of non-core collapsed Milky Way GCs.The production of BH-LMXBs is based on a subset of 15 simulatedGCs and the methods detailed in section 3.2. In the following sec-tion, we discuss the distribution and the properties of this populationof BH-LMXBs from GCs.As discussed at the end of section 3.1, some clusters requirechoosing a BH retention fraction of unity, f r BH =
1, in order toobtain the desired quantity of BHs. This occurs in the lowest-masscluster for each set of N BH , i.e. Pal 13 for N BH =
20, NGC 6838for N BH = N BH = f r BH < The number of mass transferring systems that develop from theBH-NC binaries that are ejected from our model clusters stronglydepends on the assumed BH retention in GCs. We employ the samenotation as in section 3.2.2 for BH retention: MIN refers to N BH =
20, 200 refers to N BH = N BH = + − mass-transferring BH low-mass systemsand the MAX case yields an expectation value of 156 + − ejectedBH-LMXBs, with the stated uncertainties bounding the 95 per centconfidence interval.The clusters that contribute the largest number of BH-LMXBsare those with the highest BH-NC ejection rates (see Table 3). Asis visible in Figure 4, the expected number of ejections can be ap-proximated as a function of the number of retained BHs N BH andthe two fundamental parameters describing the cluster: the concen-tration c and the total cluster mass M c . While the initial semi-majoraxis at ejection a , which is sensitive to the cluster mass (Figure 6),is an important factor in determining whether a BH-NC will leadto mass transfer, surprisingly, the fraction of BH-NCs that becomeBH-LMXBs appears nearly constant across clusters. Equivalentlystated, (cid:104) N BH − LMXB (cid:105) ∼ f LMXB (cid:104) N ej (cid:105) appears to hold true for the setof clusters modeled, where f LMXB ∼ .
25 represents the fraction ofejected BH-NC binaries that evolve into BH-LMXBs. Although thedistributions of most orbital parameters, which determine whethera system will evolve into a BH-LMXB, vary from cluster to cluster,the thermal eccentricity distribution shared by all clusters ensuresthat a roughly equal proportion of the ejected binaries will becomeBH-LMXBs. For clusters that tend to eject wider binaries, it is onlythe highly eccentric systems that become BH-LMXBs, and viceversa.
MNRAS , 1–27 (2017)
H-LMXBs from BH retaining GCs For a given BH retention, the number of successfully formedBH-LMXBs from GCs is potentially a function of the ejectiontime, initial separation, initial eccentricity, primary and companionmasses, and the complex internal evolution of the binary. Yet, sincewe find that the ejection properties are largely determined by thecluster properties, namely the quantities defining the fundamentalplane, the size of the BH-LMXB population from GCs is wellapproximated by the cluster properties alone.
As GCs generally have low escape velocities, the ejected BH-LMXBs typically escape with relatively low velocities. Due to this,the distribution of BH-LMXBs closely mimics the distribution ofGCs in the Milky Way galaxy. In Figure 9, we present the spatialprobability distribution of BH-LMXBs from GCs, for the MAXcase, on a Mollweide projection of the galactic map in an equa-torial coordinate system. Additionally, we include the distributionof galactic GCs and known BH-LMXBs from BlackCAT (Corral-Santana et al. 2016), a catalog of candidate BH-LMXBs, which weuse in all figures including an observed population, unless statedotherwise. Although the 200 case produces fewer BH-LMXBs, thedistribution is qualitatively similar to the MAX case. The highestprobability density region is near the galactic center, where themajority of GCs reside. However, as Figure 5 illustrates, the distri-butions of the ejection velocities have widths that span an order ofmagnitude or more. As a consequence, some fraction of the binarieshave ejection velocities that allow them to separate from their parentcluster. Additionally, the binaries that are ejected at an earlier timein the GC’s orbit have sufficient time to diverge from the host GCorbit. The higher density streaks in Figure 9 can be attributed tothese binaries that have drifted from the parent GC.As GCs primarily follow halo orbits that extend well out ofthe galactic plane, the GCs are easily able to populate this spacewith BH-LMXBs. In Figure 10, we provide the spatial probabilitydistribution for BH-LMXBs from the MAX case in the R − z plane.Again, we present only the MAX case, as the 200 case is simi-larly distributed but with a lower overall probability density. Themedian absolute distance from the galactic plane is | z | = .
63 kpcand the median distance from the galactic center in the plane is R = .
51 kpc. While it is clear from Figure 10 that many of theBH-LMXBs from GCs are located in the galactic disk, the distri-bution extends well out of the galactic plane into the lower densityregions above and below the disk. BH-LMXBs that form in the fieldwill generally reside in the high density galactic plane, unless theyreceive substantial kicks at birth, which might eject them into the‘high-z’ regions. However, the magnitude of BH-LMXB kicks isstill uncertain and the magnitude necessary to reach the highest ofBH-LMXBs from GCs is considered unlikely (see, e.g., Repetto &Nelemans 2015; Mandel 2016). In Figure 11, we show the cumula-tive distribution function of the absolute distance | z | perpendicularto the galactic plane for the MAX case, the 200 case, and the ob-served population of BH-LMXBs. The observed population termi-nates at a maximum | z | ∼ A typical BH-LMXB with a GC origin has an initial semi-major axisof 5 . R (cid:12) , initial BH mass of 8 . M (cid:12) , and an initial companion R [kpc] − − z [ k p c ] GCBH-LMXB − − − − − − − − d l o g P / d R d z Figure 10.
The spatial probability distribution of the simulated population ofBH-LMXBs from GCs with N BH = R − z plane. The coordinate z specifies the distance perpendicular to the galactic plane and R is thein-plane distance from the galactic center at the origin. The populations ofMilky Way GCs (marked by black circles) and known BH-LMXBs (markedby orange stars) are included for reference. While many of the BH-LMXBsfrom GCs populate the galactic disk, the distribution extends well out of thegalactic plane into the high- | z | region. | z | [kpc] F ( | z | ) Obs N BH = N BH = Figure 11.
The normalized cumulative distribution function of the absolutedistance perpendicular to the galactic plane | z | . The included distributionsare the BH-LMXBs produced in our GC simulations for the cases of N BH = N BH = N BH = mass of 0 . M (cid:12) . The median present-day period is 4 .
48 h and themedian present-day BH mass is 8 . M (cid:12) , which has increased abovethe initial median BH mass due to accretion from the companion.As discussed in section 3.2.2, the masses used in the Monte Carlomodels for the ejected binaries are sampled according to the EMFfrom the mass bin corresponding to the mass in the ejected BH-NC.This is done for both the primary BH mass M BH and the companionmass m to obtain the mass distributions, which we discuss below.In Figure 12, we show the distribution of the BH mass in the MNRAS000
48 h and themedian present-day BH mass is 8 . M (cid:12) , which has increased abovethe initial median BH mass due to accretion from the companion.As discussed in section 3.2.2, the masses used in the Monte Carlomodels for the ejected binaries are sampled according to the EMFfrom the mass bin corresponding to the mass in the ejected BH-NC.This is done for both the primary BH mass M BH and the companionmass m to obtain the mass distributions, which we discuss below.In Figure 12, we show the distribution of the BH mass in the MNRAS000 , 1–27 (2017) Giesler, Clausen, & Ott M BH [ M (cid:12) ] d P / d M B H Obs N BH = N BH = Figure 12.
The probability distributions of BH masses in BH-LMXBsfor the observed population (Özel et al. 2010) and for the BH-LMXBsproduced in our GC simulations for the cases of N BH = , N BH = M (cid:12) . population of BH-LMXBs from GCs for both cases that producemass transferring systems. Along with the BH mass distributionsfor the 200 and MAX cases, we include the inferred BH massdistribution from observations (Özel et al. 2010). Although theobserved mass distribution reaches down to ∼ M (cid:12) , our EMF doesnot produce BH masses in the range M BH < M (cid:12) . The BH primarymass is peaked at 7 . M (cid:12) and displays a preference for the lower-mass BHs. The lack of systems at high BH mass can be attributed totwo contributing factors. The leading contribution is the distributionof BH masses in the ejected BH-NCs, which is dominated by thetwo lowest BH mass bins (i.e. 8 . M (cid:12) and 20 . M (cid:12) ). Althoughthese are produced in nearly equal numbers, the preference for thelowest mass bin that arises in the BH-LMXBs is due to a secondaryeffect introduced during the binary stellar evolution. High massratio systems are prone to disrupting the companion star, endingthe possibility of evolving into a stable BH-LMXB. Despite thesebarriers to forming BH-LMXBs with high mass BHs, there remainsa small population of high mass present-day BH-LMXBs, with M BH > M (cid:12) , which accounts for ∼ m < . M (cid:12) , where the maximum mass is constrained by the MSturnoff-mass, m to = . M (cid:12) . The present-day companion massis a function of the mass-transfer rate and the time since the onset ofmass transfer. The majority of the companion masses are MS stars,however there exists a subpopulation of WD companion masseswhich account for ∼
10 per cent of the companions in the MAXcase and ∼
20 per cent in the 200 case. In Figure 13, we displaythe companion mass distribution for the MAX case, 200 case, andthe observed population of BH-LMXBs. The lack of lower-masscompanions in the 200 case relative to the MAX case is due to thehigher fraction of WDs, which have masses m WD (cid:38) . M (cid:12) . In theMAX case there is a larger number of BHs in the outskirts wherethe lowest masses reside, whereas the 200 case is more centrally m [ M (cid:12) ] d P / d m Obs N BH = N BH = Figure 13.
The probability distributions of the companion masses in BH-LMXBs for the cases N BH = N BH = concentrated where there is an increase in the probability of pickingup a higher mass companion and which includes a larger populationof WDs. The observed population in Figure 13 is generated fromthe observational data in the candidate BH-LMXB catalog Black-CAT. There are 18 confirmed BH-LMXBs in the catalog that have ameasurement of the BH mass M BH and the mass ratio q , which weuse to estimate the companion mass m = q M BH . The companionmasses in the observed population have large error bars due to theuncertainty in the measurements of the BH mass and the mass ratio.The initial eccentricity of the binaries follows a thermal dis-tribution, while the initial semi-major axis, as discussed in 4.2.1,is typically ( a / AU ) (cid:28)
1, due to their GC origin. The small initialseparation of the BH-NCs leads to a distribution of periods p where ∼
99 per cent of the BH-LMXBs have p (cid:46) . p (cid:46) . ∼
99 per cent of the populationhave p (cid:46) ∼ / MNRAS , 1–27 (2017)
H-LMXBs from BH retaining GCs Period p [ h ] d P / dp Obs N BH = N BH = N BH = N BH = Figure 14.
The probability distribution of orbital periods in the simulatedBH-LMXBs from GCs for the two stellar companion sub-populations: WDand MS. The periods for the observed population of BH-LMXBs that areless than 13 h are included for reference and are identified by orange tickmarks (18 of the 28 candidate BH-LMXBs from BlackCAT). To preserve therelative size of the MS and WD companion populations, each distributionis independently normalized and then multiplied by the factors N BH − MS / N and N BH − WD / N , respectively, with N = N BH − MS + N BH − WD . This nor-malization is applied to each N BH case independently. operating on the BH-LMXBs with a WD companion, however dueto the compact nature of WDs, the critical separation which leads toRoche lobe overflow occurs at smaller separations, hence the shorterorbital periods. The binary evolution for the BH-LMXBs from GCsis significantly different from the evolution of field binaries. In thestandard binary evolution picture, the companion evolves to over-fill its Roche lobe, which can lead to mass transfer at relativelylarge separations. The MS stars in BH-LMXBs from GCs have notevolved significantly within the cluster, but evolve on much longertimescales, preventing them from achieving mass transfer at wideseparations.In Figure 15, we provide a temperature-luminosity diagram forthe mass-transferring MS companions. We exclude the WD systemsfrom the diagram, since they are likely too faint for observation. TheMS companions have temperatures ∼ ∼ × − – 0 . (cid:12) , making these identifiable as K/M late-typeMS stars below the MS turnoff.A distinct characteristic of these systems are their kinematicproperties. In Figure 16, we show the distribution of the magnitudeof the velocity v of the BH-LMXBs from GCs. The velocity v iscomputed from the components of the space velocity in the heliocen-tric galactic coordinate system ( U , V , W ) , a right-handed coordinatesystem with U in the direction of the galactic center, V along thedirection of rotation, and W pointing toward the galactic north pole.The median values of the velocity components for the MAX caseare ( U , V , W ) = (− . , − . , − . ) km s − . The large neg-ative velocity in the V component is indicative of this populationnot participating in galactic rotation. The peculiar velocity — thevelocity of a source relative to a local standard of rest, obtained byremoving the contribution of galactic rotation at the source distancein the galactic plane R — is sometimes used to infer a ‘natal kick’ forBH-LMXBs. Although it is possible to convert the Galactic spacevelocity to a peculiar velocity, this inferred ‘kick velocity’ is only log ( T [ K ]) − − − − − − − l o g ( L [ L (cid:12) ] ) − − d l o g P / d l o g L d l o g T Figure 15.
Temperature-luminosity diagram for the BH-LMXB com-panion mass in the simulated population of BH-LMXBs from GCs with N BH = m < m to . v [ km s − ] d P / dv N BH = N BH = Figure 16.
The probability distributions for the space velocity v of thesimulated BH-LMXB population for the two BHs retention values N BH = N BH = v = v ej + v GC ,where v ej is the ejection velocity and v GC is the velocity of the host GC. Since v ej is approximately the GC escape velocity, the magnitude v is dominated bythe relatively large contribution from v GC . As such, the velocity distributionof BH-LMXBs is consistent with the velocity distribution of GCs, which isreflected in the high mean velocities. justified in assuming the source was born in the galactic disk, whereit participates in galactic rotation. For BH-LMXBs formed in thefield, which is most likely to occur in the disk, this is a reasonableassumption. However, the V component of the BH-LMXBs fromGCs indicate low rotational velocities, which is consistent with theparent GC halo orbits, which are typically non-circular and extend MNRAS000
The probability distributions for the space velocity v of thesimulated BH-LMXB population for the two BHs retention values N BH = N BH = v = v ej + v GC ,where v ej is the ejection velocity and v GC is the velocity of the host GC. Since v ej is approximately the GC escape velocity, the magnitude v is dominated bythe relatively large contribution from v GC . As such, the velocity distributionof BH-LMXBs is consistent with the velocity distribution of GCs, which isreflected in the high mean velocities. justified in assuming the source was born in the galactic disk, whereit participates in galactic rotation. For BH-LMXBs formed in thefield, which is most likely to occur in the disk, this is a reasonableassumption. However, the V component of the BH-LMXBs fromGCs indicate low rotational velocities, which is consistent with theparent GC halo orbits, which are typically non-circular and extend MNRAS000 , 1–27 (2017) Giesler, Clausen, & Ott . × − . × − . × − . × − . × − NGC 6838 20.61 8.27 1.02 8 . × − . × − . × − . × − . × − . × . × − . × − . × − . × − . × − . × − . × − . × − . × − . × . × − . × − . × − . × − . × − . × − . × − . × − . × . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − NGC 6121 20.70 1 . × . × − . × − . × − . × − . × . × . × − . × − . × . × . × . × − . × − . × . × NGC 5694 20.49 2 . × . × − . × − . × . × . × − . × − . × . × . × . × . × . × − . × . × . × NGC 6093 19.85 3 . × . × − . × − . × . × . × − . × . × . × . × . × . × − . × . × . × NGC 5286 12.29 1 . × . × − . × − . × − . × . × . × − . × − . × . × . × . × . × − . × . × . × NGC 6656 19.80 1 . × . × − . × − . × − . × . × . × − . × − . × . × . × . × . × NGC 1851 20.76 2 . × . × − . × . × . × − . × − . × . × . × . × . × . × − . × . × . × NGC 6205 20.10 1 . × . × − . × − . × − . × . × . × − . × − . × . × . × . × − . × . × NGC 6441 20.98 2 . × . × − . × − . × . × . × − . × . × . × . × . × . × . × . × NGC 104 22.49 2 . × . × − . × − . × . × . × − . × . × . × . × . × . × − . × . × . × NGC 5139 20.84 7.15 8 . × − . × − . × − . × − . × − . × . × − . × − . × . × − . × . × . × − . × − . × . × Table 4.
Expected number of mergers. For each cluster and number of retained BHs, we list the exact number of BHs in the cluster along with the expectednumber of mergers over the cluster lifetime. well out of the galactic plane. As the BH-LMXBs with GC originsare ejected at relatively low velocities along the GC’s orbit in thegalaxy, this population of BH-LMXBs has a velocity distributionconsistent with the high-velocity halo orbits of GCs. As these sys-tems have high apparent peculiar velocities, due to their halo orbitsand the lack of participation in galactic rotation, attempting to infera ‘natal kick’ from the peculiar velocity in such a case is ill-posedand leads to the conclusion of a large required ‘natal kick.’
As briefly discussed in section 2.3, we allow for gravitational radia-tion driven mergers between compact objects. Since all of our ‘testbinaries’ contain at least one BH, the allowable set of GW mergerpairs is limited to BH-NS, BH-WD, and BH-BH. In addition tothose binaries that merge during their evolution within the cluster,binaries of these types can also be ejected from the cluster. In the case of the ejection of a compact pair, we calculate the expectedmerger time t d using the ejected binary parameters and refer tothese as post-ejection mergers if t ej + t d < t H , where t H = yr isapproximately the Hubble time. The total merger rate includes thesepost-ejection mergers in addition to the in-cluster mergers. Here wepresent an estimate of the merger rates averaged over the 10 yrsimulations for different BH retention values.For notational convenience, we refer to a parameter set as x i ,where the index i runs over the 39 parameter sets which make upeach row of Table 2 and corresponds to a specific GC and valueof N BH . We compute the expected number of mergers for eachparameter set by considering the probability of a BH being involvedin a merger, defined simply by P m ( x i ) = N mergers ( x i ) N runs ( x i ) , multiplied bythe BH population (cid:104) N m (cid:105) i = P m ( x i ) N BH ( x i ) . (45)In the case of a merger involving two BHs, the expectation value iscalculated using N BH ( x i )/ MNRAS , 1–27 (2017)
H-LMXBs from BH retaining GCs (cid:104) R ( N BH )(cid:105) BH-BH BH-NS BH-WD (cid:104) R ( )(cid:105) . × − . × − . × − (cid:104) R ( )(cid:105) . × − . × − . × − (cid:104) R ( )(cid:105) Table 5.
The contribution to the compact merger rate density from all GCsin the universe, stated in Gpc − yr − . Each row corresponds to the mergerrate contribution from GCs with the simulated BH population specified by N BH in (cid:104) R ( N BH )(cid:105) . The merger rate densities are averaged over the lifeof the cluster, weighted by the GC mass function to account for the non-uniform mass distribution of GCs, and assumes a GC spatial density of ρ GC = .
77 Mpc − . driven compact object mergers over the lifetime of each cluster fora given BH population. The number of BH-BH mergers is stronglycorrelated with the GC core density n o . Each population of BHshas a merger expectation value that follows a power-law in the coredensity with exponent ∼ .
58. Since we do not include primordialbinaries, exchange encounters are the only means to forming BH-BH binaries that can later merge. The average rate of encounters isdirectly proportional to the density, with the highest density clustersproviding the largest number of opportunities to successfully formBH-BH binaries. There are additional correlated variables, such asthe concentration c and velocity dispersion σ , however these aresecondary to the density n o and likely due to their own correlationwith n o .Given the expected number of mergers for each cluster, we de-termine a weighted average using the GC mass function, since thetotal cluster mass of GCs is not uniformly distributed (McLaughlin& Pudritz 1996). We do this individually for each group of simu-lations belonging to the sets N BH = { , , } , utilizing theGC mass spectrum dN ( M c )/ dM c of McLaughlin & Pudritz (1996).For each simulated cluster, we assign a weight w i = N ( M c ( x i )) andcompute the expected number of mergers per cluster in the MilkyWay from (cid:104) N m ( N BH )(cid:105) = (cid:205) i w i (cid:104) N m (cid:105) i (cid:205) i w i . (46)For clarity, to obtain the expected number of mergers for N BH = N BH =
20. Theresulting expected number of BH-BH mergers over the life of acluster for each choice of N BH are (cid:104) N m ( )(cid:105) = . (cid:104) N m ( )(cid:105) = .
08, and (cid:104) N m ( )(cid:105) = . t GC = yrs old, and that the spatial densityof GCs in the universe is ρ GC = .
77 Mpc − (see supplementalmaterials of Rodriguez et al. 2015). Using the weighted averagescomputed above as our ‘typical’ cluster merger values and assigningthis value to each GC in the volume, we obtain the merger ratedensity due to all GCs in the universe, (cid:104) R ( N BH )(cid:105) = (cid:104) N m ( N BH )(cid:105) t GC ρ GC . (47)In Table 5, we provide the computed estimated merger rate densitiesfor compact object mergers due to GCs for the three populations of N BH we consider. Although there is an increased interest in theBH-mass spectrum for BH-BH mergers in GCs, stimulated by thelarger than expected BH masses recently detected by aLIGO (Abbottet al. 2016d), the use of just three discrete BH masses precludes thepossibility of such an analysis. Since BH-BH mergers from GCs only partially contribute tothe total merger rate, with the remaining mergers coming fromthe field, the rates due to GCs should not exceed the upper boundof the total estimated merger rate. The most recent observationalevidence constrains the BH-BH merger rate density to lie in therange 12 −
213 Gpc − yr − (Abbott et al. 2017). The BH-BH mergerrate densities given in Table 5 for the three different BH retentionscenarios are well below the upper bound, presenting no conflictwith the observed rate. However, since the fraction of the totalmergers attributable to GCs is still largely uncertain, none of theof GC rates presented here can be ruled out based on the currentobserved rate of BH-BH mergers.The bounds of our merger rates, which span a wide range ofuncertainty in BH retention, are consistent with previous studies thatprovide estimates of the BH-BH merger rate from GCs (O’Learyet al. 2006; Sadowski et al. 2008; Downing et al. 2011; Morscheret al. 2015, Rodriguez et al. 2016a). However, we find that only ∼
10 per cent of the BH-BH mergers occur outside of the clusterboundaries, which differs from a subset of these previous stud-ies. In Downing et al. (2011) no mergers occur in-cluster, whilein Morscher et al. (2015), ∼
85 per cent of BH-BH mergers occurpost-ejection, and Rodriguez et al. (2016a) find that ∼
90 per centmerge outside the cluster. In contrast to the small number of BH-BH binaries these studies find merging in cluster, O’Leary et al.(2006) finds that only ∼ −
72 per cent of the BH-BH mergers arepost-ejection. Finally, Sadowski et al. (2008) is most closely alignedwith our results, with ∼
10 per cent of mergers occurring out of thecluster.This discrepancy in merger location can be attributed to thedistribution of the BHs in the cluster and their interactions with thelower-mass components. In models with centrally clustered BHs,the BHs are segregated from the remainder of the cluster, formingan isolated and decoupled system. These self-interacting BHs ef-ficiently form BH binaries. Strong binary-binary interactions caneject these binary BHs from the cluster, where they might latermerge in isolation. In addition to the efficient removal of BH bina-ries from the core, binary-single interactions are equally efficientat ejecting single BHs from the cluster. Furthermore, these strongencounters are likely to interrupt potential mergers of eccentric BHbinaries which would merge in-cluster if uninterrupted. This chan-nel leads to a majority of BH-BH mergers outside of the cluster andeventually depletes the GC of BHs (e.g., O’Leary et al. 2006; Baner-jee et al. 2010; Downing et al. 2011). We assume that in order forGCs to retain significant BH populations, the BHs must avoid segre-gating in the core, which we accomplish through a modified velocitydispersion for the BHs, as discussed in section 2.2.4. This modifiedvelocity dispersion spreads the BHs throughout the cluster, wherethey can interact with the lower-mass stars. This supposition is sim-ilar to the assumptions made in Sadowski et al. (2008) and producesqualitatively similar results.In our simulations, a key channel for producing BH-BH bina-ries is through the formation of a binary composed of a BH anda non-BH outside of the core, which eventually drift to the centerwhere there is a high density of BHs. The non-BH will be prefer-entially exchanged with one of the more massive BHs in the core,producing a BH-BH binary that will realize one of three outcomes:(1) the BH-BH binary will be dismantled in the high density region,(2) given a sufficiently large eccentricity (hence a shorter orbital de-cay time), will eventually merge in the core, or (3) will harden andbe ejected from the cluster. This formation channel is similar to thatdescribed in Sadowski et al. (2008). As discussed in section 2.3.5,we allow for single BHs to exchange into existing binaries. The ma-
MNRAS000
MNRAS000 , 1–27 (2017) Giesler, Clausen, & Ott − − − − log ( − e ) d P / d l o g ( − e ) BH-BH m : e ( t o ) BH-BH m : e ( t m ) BH-BH ej : e ( t o ) BH-BH ej : e ( t ej ) f ( e ) = e Figure 17.
The probability distributions of eccentricity for two populationsof BH-BH binaries in GCs: BH-BH binaries which form and merge in cluster(BH-BH m , black lines) and the BH-BH binaries which form and are ejectedfrom the cluster (BH-BH ej , blue lines). For each population, we show theeccentricity distribution at the time the binary forms, e ( t o ) (solid lines), andthe distribution of eccentricities at the binary’s final state (dashed lines).The final state of the in-cluster mergers is at a time t m , the time at whichthe computed merger time is less than the cluster timestep. The final statefor the ejected binaries is the time of ejection t ej . A thermal eccentricitydistribution, with probability density f ( e ) = e , is included for reference. jority of binaries that a single BH encounters are binaries composedof two low-mass stars. Successful exchanges of a more massive BHfor one of the lower-mass stars tend to produce high-eccentricityBH–non-BH binaries following the relation (cid:104) e (cid:105) ≈ − . (cid:0) m non − BH m BH (cid:1) , (48)which is independent of the initial eccentricity and applicable when m non − BH (cid:28) m BH (Sigurdsson & Phinney 1993). For the three BHmasses considered, M BH = { . , . , . } M (cid:12) , and a clusternon-BH star with an average mass of (cid:104) m non − BH (cid:105) ≈ . (cid:12) , thisleads to mean initial eccentricities of (cid:104) e (cid:105) ≈ { . , . , . } .Once the binary makes it to the core, the non-BH is easily exchangedfor one of the many massive BHs, yielding a highly eccentric BH-BH binary according to Equation 48. In Figure 17, we display theeccentricity distributions for the BH-BH binaries at formation andat merger or ejection for those binaries that have end states (2) and(3), as described above, respectively. Some fraction of the eccentricbinaries that form through this channel are driven to high enougheccentricities that they can merge in-cluster in-between encounters.The remainder are subject to further encounters that drive theireccentricities toward a thermalized distribution, are hardened in theprocess, and are eventually ejected.The eccentricity distribution of merging BH-BH binaries isimportant for the detection of the resulting gravitational waves. Theeccentricity tends to zero as the orbit shrinks, however moderndetectors are sensitive to the GW signal at frequencies when thebinary is still in the inspiral phase and the eccentricity is finite.The aLIGO (LIGO Scientific Collaboration et al. 2015) detectorsare sensitive to ∼
10 Hz, at design sensitivity, while the futurespace-based detector LISA (Amaro-Seoane et al. 2013) will besensitive to much lower frequencies ∼ − − − − − − − − log e d P / d l o g e BH-BH in : f GW =
10 Hz : aLIGOBH-BH in : f GW = out : f GW =
10 Hz : aLIGOBH-BH out : f GW = Figure 18.
The eccentricity probability distributions for two populations ofBH-BH mergers from GCs for the two detectors aLIGO and LISA. The twopopulations correspond to the BH-BH mergers occurring in-cluster (solidlines) and those that merge outside of the cluster, post-ejection (dashed lines).The black lines correspond to the eccentricity of each population when itreaches a corresponding gravitational wave frequency of f GW =
10 Hz, thelower bound frequency of the aLIGO band at design sensitivity. The bluelines represent the eccentricity distribution at f GW = centricity at a specific frequency by evolving a o and e o , accordingto (cid:104) de / da (cid:105) (Peters 1964), up until some target value a associatedwith the frequency in consideration. In Figure 18, we display theresidual eccentricity of the inspiraling BH-BH binaries as they firstenter the design-sensitivity frequency bands for aLIGO and LISA.It is apparent that for aLIGO, both the ejected mergers and the ini-tially high-eccentricity in-cluster mergers have residual eccentricitydistributions below 10 − , which has a negligible effect on detec-tions using circularized templates. However, in the case of LISA,while the ejected mergers result in a small eccentricity at 1 mHz,the initially highly eccentric in-cluster merger population remainssignificantly eccentric at this frequency.Utilizing (cid:104) de / da (cid:105) to determine the evolved eccentricity as-sumes that the binary evolves in isolation. For the in-cluster merg-ers, we classify a BH-BH binary as merged once the orbital decaytime has fallen below the cluster timestep. However, this couldleave significant time for further dynamics to modify the eccentric-ity such that the binary will not in fact merge in cluster (Banerjeeet al. 2010). To account for this possibility, the in-cluster mergersin Figure 18 only include those mergers which satisfy the addi-tional constraint t dec < (cid:104) t enc (cid:105) , which is satisfied for ∼
70 per centof in-cluster mergers. Here, the average encounter time is approx-imated by (cid:104) t enc (cid:105) = t bin / N enc with t bin corresponding to the timesince the binary’s formation and N enc is the number of three-bodyencounters the binary has been subject to during the time t bin . Theremaining ∼
30 per cent of mergers are uncertain and are not furtherevolved; they may be broken up, ejected, or merge after subsequentinteractions.
MNRAS , 1–27 (2017)
H-LMXBs from BH retaining GCs (cid:104) R ( N BH )(cid:105) BH-NC BH-WD BH-NS BH-BH (cid:104) R ( )(cid:105) . × − . × − . × − . × − (cid:104) R ( )(cid:105) .
08 1 . × − . × − . × − (cid:104) R ( )(cid:105) .
27 2 .
14 3 . × − . × − Table 6.
The rate of three-body mergers in GCs computed for the Milky Waygalaxy and stated in MWEG − Myr − . Each row corresponds to the three-body merger rate in Milky Way GCs with the simulated BH populationspecified by N BH in (cid:104) R ( N BH )(cid:105) . The merger rates are averaged over the lifeof the cluster, weighted by the GC mass function to account for the non-uniform mass distribution of GCs, and assumes N GC (cid:39)
150 for the numberof GCs in the galaxy.
In addition to the GW-driven mergers, we also calculate the rate oftidally driven mergers or ‘collisions’ that occur during three-bodyencounters. The merger criteria are based on a minimum separationbetween bodies, as discussed in section 2.3.4. We compute theexpected number of three-body merger events only for those thatinvolve a BH. Although we track the number of three-body mergersfor all object types, including NS-NS, MS-WD, etc., we are missinga significant fraction of these mergers by only tracking single BHsor binaries with at least one BH. We compute the expected numberof mergers in a manner similar to the computation of GW mergersabove.The left columns of Table 4 list the expected number of mergersinvolving a BH that occur during three-body encounters over thelifetime of each cluster for a given BH population. These three-bodymergers are computed using Equation 45 to obtain an expected valuefor each cluster in the set. As the majority of these events will onlybe observationally relevant locally, we provide these rates solely forthe Milky Way galaxy. Using the computed values from Table 4we construct a cluster weighted average with Equation 46. Fromthis we use a modified version of Equation 47, with N GC (cid:39) ρ GC to obtain the final approximate rate for each event: (cid:104) R ( N BH )(cid:105) = (cid:104) N m ( N BH )(cid:105) t GC N GC . These computed rates for BH-BH, BH-NS, BH-WD and BH-NC are shown in Table 6, stated in terms of the numberof expected events per Milky Way equivalent galaxy (MWEG) perMyr. The BH-NC merger rate includes the three-body mergers ofboth BH-RG and BH-MS.These rates are included to ensure that a large population ofretained BHs in GCs does not lead to a conflict with observations.Even in the case of maximal BH retention, the occurrence of theseevents is relatively infrequent. The most commonly occurring three-body collision is that between a BH and a NC star. The interactionof a NC object with a BH, commonly referred to as a tidal disruptionevent (TDE), is often studied in the context of supermassive BHsrather than stellar-mass BHs. However, there is some interest in GC-relevant NC collisions with stellar-mass BHs, which are referred toas micro-TDEs (Perets et al. 2016). These events lead to full orpartial tidal disruption of the NC star and are accompanied by long-duration energetic flares. There is large uncertainty in the signalsassociated with these events as the strength and duration of the signaldepends heavily on the details of the encounter (see, e.g., Perets et al.2016).The signals associated with the compact mergers are likely toappear as head-on mergers due to the criteria associated with cat-egorizing mergers during three-body encounters; the exclusion ofhigher order corrections to Newtonian gravity in our three-body cal- culations requires extremely close-encounters due to the relativelysmall size of the compact objects involved. Despite the uncertaintyin the observables produced in three-body collisions, the rate of oc-currence is low enough that our model does not generate a conflictwith present observations. In our simulations, GCs produce a population of BH-LMXBs witha unique set of characteristic properties. These properties providesome constraints on the likelihood of a BH-LMXB having a GCorigin. In this section, we identify the key characteristics of BH-LMXBs from GCs and determine which of the currently knownBH-LMXBs are consistent with this population.As discussed in section 4.2.3 and visible in Figure 12, the spec-trum of BH masses in BH-LMXBs from GCs in our simulationsis roughly consistent with the observed population of BH masses.This makes the BH mass a poor candidate for differentiating be-tween field-formed BH-LMXBs and those with a GC origin. As aconsequence of the age of GCs, the companions are typically un-evolved MS stars, with masses necessarily below the turnoff-mass m to = . M (cid:12) . Additionally, they reside on a tightly confinedbranch of a temperature-luminosity diagram (see Figure 15). Thisprovides the first distinctive characteristic of BH-LMXBs formed inGCs: a companion mass of m (cid:46) . M (cid:12) and a spectral class con-sistent with late-type K/M stars. BlackCAT (Corral-Santana et al.2016) currently contains 18 observed BH-LMXB systems with theproper information to compute an estimate of the companion mass.Of the 18 systems, six BH-LMXBs have companion masses exceed-ing the maximum companion mass in our population of BH-LMXBsfrom GCs. Two of these six are near the edge of the distribution withwith m (cid:38) . M (cid:12) , while the other four have m (cid:62) . M (cid:12) , sug-gesting these are more consistent with a field-formation scenario.A second property of a BH-LMXB with a GC origin is acharacteristically short period. As shown in Figure 14, there is asharp limit in the distribution confining GC-origin BH-LMXBs toperiods shorter than p ∼ . p > | z | perpendicular to the galactic plane (seeFigure 11), the overall distribution of the BH-LMXBs from GCsdoes not provide a strict criterion for discerning between GC ori-gin and field origin. Figure 10 illustrates that while the simulatedpopulation extends much farther out of the galactic plane than theobserved distribution, there is still a significant population of GC-origin BH-LMXBs that reside in the plane, overlapping the regionwhere field-formed binaries are expected to have the highest den-sity. This makes discerning a potential origin for BH-LMXBs inthis region difficult. Additionally, for the many systems clusterednear the galactic center or those that reside in the plane, the highdensity of objects and dust make these systems equally difficult toobserve optically. Although a number of BH-LMXB candidates aredetectable in these regions through X-ray, the detailed propertiesof these systems remain unknown due to current optical limita-tions. The spatial distribution of BH-LMXBs from GCs, in general,makes observations of the population difficult, even for those outof the plane. Observation and confirmation of BH-LMXBs rely ona dynamical measurement of the BH mass through optical spec-troscopy, introducing a bias toward sources at distances D <
10 kpc
MNRAS000
MNRAS000 , 1–27 (2017) Giesler, Clausen, & Ott from the Sun (Repetto & Nelemans 2015). For the population ofBH-LMXBs from our model GCs, the MAX and 200 cases bothproduce a median distance of D = . ( U , V , W ) introduced in section 4.2.3, Mirabel et al. (2001) found aspace-velocity for this system of ( U = − ± , V = − ± , W = − ± ) km s − . The large magnitude v ∼
145 km s − and the largenegative V component are consistent with a high-velocity halo orbitand a lower than average rotational velocity about the galactic cen-ter. This description is consistent with the velocity distribution ofour population of BH-LMXBs from GCs, which inherit the high-velocity halo-orbits when they are ejected from the GC. As a con-sequence of the high-velocity halo orbit, which manifests itself as ahigh computed peculiar velocity, this system is commonly invokedto support large natal kicks (Gualandris et al. 2005; Fragos et al.2009; Repetto et al. 2012; Repetto & Nelemans 2015). Confidentlyidentifying an origin for this system could help to shed some light onthe issue. The relatively low-metallicity environments of GCs pro-vides an additional constraint on properly categorizing BH-LMXBsas originating in GCs versus in the field. Although all of the previouscharacteristics point to a GC origin, perhaps one of the strongestarguments against a GC origin for this system is the supersolarabundance of elements in the secondary star found by GonzálezHernández et al. (2006), which is consistent with a metal-rich pro-genitor and makes a GC origin highly unlikely. However, thereexist a conflicting claim presented by Frontera et al. (2001), wherethrough broad-band X-ray spectroscopy, it was concluded that thecompanion has a metallicity of Z / Z (cid:12) ∼ .
1, consistent with thelow metallicities expected of systems at large | z | or those with a GCorigin. Given that metallicity provides a strong constraint on theorigin of a BH-LMXB, additional observations appear necessary toreduce the uncertainty of this case.To our knowledge, there are no known velocity measurementsor metallicity measurements for the four other BH-LMXBs withpossible GC origins. Although an increasing number of BH-LMXBcandidates are being discovered in X-rays, only a few have beenconfirmed and characterized with detailed optical follow-up obser-vations. Over time, more data will become available, better con-straining the properties of the galactic BH-LMXB population. If even a single BH-LMXB could be confidently attributed to a GCorigin this would provide a strong argument in favor of BH retentionin GCs. There is growing observational evidence and theoretical support fora sizable BH population in present-day galactic GCs. These BHscan acquire low-mass companions through dynamical interactionswithin the GC. Those binaries that are ejected from the GC canevolve into BH-LMXBs and can populate a large region of spaceabove and below the galactic plane. These binaries could potentiallyexplain observed BH-LMXBs at large distances from the planewithout a need for large BH birth kicks.In this study, we have presented a population of Milky WayBH-LMXBs formed through dynamical interactions in GCs. Toexplore the BH-LMXB population dependence on BH retentionin GCs, we performed simulations for retained BH populationsof 20, 200, and 1000 BHs. The simulated GCs broadly cover theparameter space and represent a realistic subset of Milky Way GCs.We generated a large number of binary evolution realizations foreach set of initial GC parameters and number of retained BHs. Thisallowed us to derive statistical distributions for the number of ejectedbinaries and their relevant properties. Using the statistics from theGC simulations, we performed Monte Carlo simulations to obtaina present day population of BH-LMXBs ejected from GCs.We find that in the case of minimal BH retention ( N BH =
20) no observable BH-LMXBs are produced, while the N BH =
200 and N BH = + − and 156 + − BH-LMXBs,respectively. Here, the uncertainties represent the bounds of the 95per cent confidence interval. As there is no observable populationfor minimal BH retention, this suggests that finding any BH-LMXBof GC origin would imply that GCs retain sizable BH populationsof more than a few tens of BHs.Aside from the difference in the size of the population, theproperties and distributions of BH-LMXBs are qualitatively similarfor the two cases that produce BH-LMXBs, 200 and MAX. Wefind that BH-LMXBs from GCs have velocity distributions inher-ited from their host clusters that are consistent with stars on high-velocity halo orbits. Additionally, the ejected BH-LMXBs have aspatial distribution that is also similarly aligned with the GC galacticdistribution. This shared distribution is described by a high densityin the galactic plane and near the galactic center, with a significantfraction distributed well above and below the galactic plane. Thetypical binary is located at an absolute distance of R = . | z | = . D = .
74 kpc from the Sun. The presenceof a large population of BH-LMXBs at large distances from theplane is characteristic of BH-LMXBs from GCs, as field formedBH-LMXBs must be subject to large kicks in order to access thisregion. The average present-day BH-LMXB ejected from a GC iscomposed of a 8 . M (cid:12) BH and a 0 . M (cid:12) K/M late-type MS starbelow the turnoff-mass, with a characteristically short orbital periodof p = .
186 h. These properties and their associated distributionsare key observable characteristics of this predicted population ofBH-LMXBs formed in GCs.Comparing our BH-LMXB systems with the ensemble of ob-served BH-LMXBs, we find that five of these are candidates for hav-ing a GC origin. There are a total of 27 confirmed BH-LMXBs, butjust 18 of these have sufficient observations for comparing measured
MNRAS , 1–27 (2017)
H-LMXBs from BH retaining GCs Name M BH [ M (cid:12) ] m [ M (cid:12) ] p [ h ] | z | [ kpc ] ReferencesMAXI J1659-152 5 . ± . . ± .
05 2 . ± × − . ± .
05 [0,1]SWIFT J1357.2-0933 > . > .
33 2 . ± × − > .
75 [2,3]SWIFT J1753.5-0127 > . ± . (cid:62) . ± .
03 3 . ± × − . ± . . ± .
65 0 . ± .
083 4 . ± × − . ± .
09 [8-11]GRO J0422+32 8 . ± . . ± .
31 5 . ± × − . ± .
06 [12-15]
Table 7.
Properties of the five observed systems that are consistent with the properties of our simulated population of BH-LMXBs with GC origins. Thecolumns refer to the primary BH mass M BH , the companion mass m , the orbital period p , and the absolute distance perpendicular to the galactic plane | z | .[0] Yamaoka et al. (2012), [1] Kuulkers et al. (2013), [2] Mata Sánchez et al. (2015), [3] Corral-Santana et al. (2013), [4] Shaw et al. (2016), [5] Neustroevet al. (2014), [6] Zurita et al. (2008), [7] Cadolle Bel et al. (2007), [8] Khargharia et al. (2013), [9] Calvelo et al. (2009), [10] Torres et al. (2004), [11] Gelinoet al. (2006), [12] Casares et al. (1995), [13] Beekman et al. (1997), [14] Webb et al. (2000), [15] Gelino & Harrison (2003) properties against our results. The five systems that are compatiblewith our simulated population of BH-LMXBs from GCs are SWIFTJ1357.2-0933, SWIFT J1753.5-0127, XTE J1118+480, and GROJ0422+32. XTE J1118+480 is one of the rare systems with a mea-sured space velocity and it is atypically large for a system formedin the galactic disk, with v ∼
145 km s − . This system is commonlydiscussed in the context of formation kicks, since a high-velocitykick is required to explain the large distance from the galactic plane, | z | ∼ .
52 kpc, under the assumption that it originated in the plane.However, if XTE J1118+480 comes from a GC, which producesBH-LMXBs at a median distance of | z | ∼ . N BH = .
81 Gpc − yr − , while in the case ofminimal retention, N BH =
20, the rate is as low as 3 . × − Gpc − yr − . This rate represents an average over the cluster life-times and assumes a spatial density of GCs throughout the universeof ρ GC = .
77 Mpc − . Our maximum retention rate is consistentwith previous estimates of the GC merger rate contribution andis compatible with the recent observations by aLIGO. Althoughour model produces rates in good agreement with previous studies,our simulations result in a larger than average fraction of merg-ers occurring in-cluster, as opposed to post-ejection. We attributethe discrepancy to the increased interaction between the BHs andthe lower mass stars as a consequence of our cluster BH distribu-tion. The BH-BH binaries that merge in-cluster are a consequenceof the large eccentricities, acquired through dynamical formation, leading to significantly shortened orbital decay times. The dynami-cally formed BH-BH binaries that merge in-cluster are formed withan average eccentricity of e ∼ .
96. At the time of merger in theaLIGO band, the residual eccentricities are small and in the range10 − (cid:46) e (cid:46) − . However, we find that when passing through theLISA band years before merger, they still have eccentricities in therange 10 − (cid:46) e (cid:46)
1. Models in which the BHs are confined to asubcluster at the core of GCs produce mergers with substantiallysmaller eccentricities. As the merger formation channels are suffi-ciently different for a BH subcluster model, LISA might be able tohelp distinguish how a population of retained BHs is distributed inGCs by observing the distribution of eccentricities.The present study provides new insights into the populationand properties of BH-LMXBs of GC origin. However, there are anumber of important limitations that should be kept in mind wheninterpreting our results. While there is mounting evidence to supportthat present-day GCs are BH retaining, how GCs are able to retaina significant population of BHs and how those BHs are distributedis still uncertain. Our choice of distributing the BHs throughout thecluster is motivated by preserving the observed structural proper-ties of each modeled GC in the presence of a large BH population.However, this spreading leads to an increase in interaction betweenthe BHs and the lower-mass stars, which is typically a rare occur-rence if the BHs remain clustered in the core. If GCs are able toretain a significant population of BHs that remain centrally clus-tered, formation of BH-NC binaries will likely be suppressed. Thereduced formation of BH-NC binaries would significantly reducethe number of ejected BH-NCs, directly diminishing the numberof BH-LMXBs from GCs. Future studies regarding the impact ofthe BH distribution within BH-retaining GCs are necessary to fullyunderstand the consequences of this limitation. Furthermore, theresults presented here rely on the outcomes of many independentrealizations. Since we perform each simulation independently in astatic cluster background, we are neglecting the change in the BHpopulation and its impact on the cluster as single BHs and BHbinaries are ejected over the cluster lifetime. Additionally, we donot account for binary-binary interactions, which have the poten-tial to disrupt existing binaries or possibly aid in ejecting them.Models which account for these limitations are necessary to betterunderstand the impact of ignoring these processes. While N -bodysimulations and Monte Carlo based models can resolve some ofthese issues, the computational expense remains a limiting factor inperforming many realizations. However, as the computational tech-niques and resources continue to improve, it will soon be possibleto produce many high-accuracy GC simulations that address theselimitations. MNRAS000
1. Models in which the BHs are confined to asubcluster at the core of GCs produce mergers with substantiallysmaller eccentricities. As the merger formation channels are suffi-ciently different for a BH subcluster model, LISA might be able tohelp distinguish how a population of retained BHs is distributed inGCs by observing the distribution of eccentricities.The present study provides new insights into the populationand properties of BH-LMXBs of GC origin. However, there are anumber of important limitations that should be kept in mind wheninterpreting our results. While there is mounting evidence to supportthat present-day GCs are BH retaining, how GCs are able to retaina significant population of BHs and how those BHs are distributedis still uncertain. Our choice of distributing the BHs throughout thecluster is motivated by preserving the observed structural proper-ties of each modeled GC in the presence of a large BH population.However, this spreading leads to an increase in interaction betweenthe BHs and the lower-mass stars, which is typically a rare occur-rence if the BHs remain clustered in the core. If GCs are able toretain a significant population of BHs that remain centrally clus-tered, formation of BH-NC binaries will likely be suppressed. Thereduced formation of BH-NC binaries would significantly reducethe number of ejected BH-NCs, directly diminishing the numberof BH-LMXBs from GCs. Future studies regarding the impact ofthe BH distribution within BH-retaining GCs are necessary to fullyunderstand the consequences of this limitation. Furthermore, theresults presented here rely on the outcomes of many independentrealizations. Since we perform each simulation independently in astatic cluster background, we are neglecting the change in the BHpopulation and its impact on the cluster as single BHs and BHbinaries are ejected over the cluster lifetime. Additionally, we donot account for binary-binary interactions, which have the poten-tial to disrupt existing binaries or possibly aid in ejecting them.Models which account for these limitations are necessary to betterunderstand the impact of ignoring these processes. While N -bodysimulations and Monte Carlo based models can resolve some ofthese issues, the computational expense remains a limiting factor inperforming many realizations. However, as the computational tech-niques and resources continue to improve, it will soon be possibleto produce many high-accuracy GC simulations that address theselimitations. MNRAS000 , 1–27 (2017) Giesler, Clausen, & Ott
ACKNOWLEDGMENTS
The authors thank Sterl Phinney, Steinn Sigurdsson, and SaulTeukolsky for valuable discussions. This work is partially supportedby the Sherman Fairchild Foundation and by NSF under awardNo. CAREER PHY-1151197. The simulations were carried out onNSF/NCSA Blue Waters under PRAC award no. ACI-1440083 andon the Caltech cluster Zwicky, supported by the Sherman FairchildFoundation and NSF award No. PHY-0960291.
REFERENCES
Aarseth S. J., 2012, MNRAS, 422, 841Abbott B. P., et al., 2016a, Phys. Rev. Lett., 116, 061102Abbott B. P., et al., 2016b, Phys. Rev. Lett., 116, 241103Abbott B. P., et al., 2016c, ApJ, 818, L22Abbott B. P., et al., 2016d, ApJ, 818, L22Abbott B. P., et al., 2017, Phys. Rev. Lett., 118, 221101Amaro-Seoane P., et al., 2013, GW Notes, Vol. 6, p. 4-110, 6, 4Bahcall N. A., Hausman M. A., 1977, ApJ, 213, 93Banerjee S., Baumgardt H., Kroupa P., 2010, MNRAS, 402, 371Beekman G., Shahbaz T., Naylor T., Charles P. A., Wagner R. M., MartiniP., 1997, MNRAS, 290, 303Belczynski K., Sadowski A., Rasio F. A., Bulik T., 2006, ApJ, 650, 303Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition. PrincetonUniversity PressBlaauw A., 1961, Bull. Astron. Inst. Netherlands, 15, 265Bovy J., 2015, ApJS, 216, 29Breen P. G., Heggie D. C., 2013, MNRAS, 432, 2779Cadolle Bel M., et al., 2007, ApJ, 659, 549Calvelo D. E., Vrtilek S. D., Steeghs D., Torres M. A. P., Neilsen J., Filip-penko A. V., González Hernández J. I., 2009, MNRAS, 399, 539Campanelli M., Lousto C., Zlochower Y., Merritt D., 2007, ApJ, 659, L5Casares J., Jonker P. G., 2014, Space Sci. Rev., 183, 223Casares J., Martin A. C., Charles P. A., Martin E. L., Rebolo R., HarlaftisE. T., Castro-Tirado A. J., 1995, MNRAS, 276, L35Catalán S., Isern J., García-Berro E., Ribas I., 2008, MNRAS, 387, 1693Chernoff D. F., Weinberg M. D., 1990, ApJ, 351, 121Clausen D., Wade R. A., Kopparapu R. K., O’Shaughnessy R., 2012, ApJ,746, 186Clausen D., Sigurdsson S., Chernoff D. F., 2013, MNRAS, 428, 3618Cohn H., 1979, ApJ, 234, 1036Corral-Santana J. M., Casares J., Muñoz-Darias T., Rodríguez-Gil P., Shah-baz T., Torres M. A. P., Zurita C., Tyndall A. A., 2013, Science, 339,1048Corral-Santana J. M., Casares J., Muñoz-Darias T., Bauer F. E., Martínez-Pais I. G., Russell D. M., 2016, A&A, 587, A61Downing J. M. B., Benacquista M. J., Giersz M., Spurzem R., 2011, MN-RAS, 416, 133Duric N., 2004, Advanced astrophysicsFragos T., Willems B., Kalogera V., Ivanova N., Rockefeller G., Fryer C. L.,Young P. A., 2009, ApJ, 697, 1057Fregeau J. M., Ivanova N., Rasio F. A., 2009, ApJ, 707, 1533Frontera F., et al., 2001, ApJ, 561, 1006Gelino D. M., Harrison T. E., 2003, ApJ, 599, 1254Gelino D. M., Balman Ş., Kızıloˇglu Ü., Yılmaz A., Kalemci E., TomsickJ. A., 2006, ApJ, 642, 438González Hernández J. I., Rebolo R., Israelian G., Harlaftis E. T., FilippenkoA. V., Chornock R., 2006, ApJ, 644, L49Gualandris A., Colpi M., Portegies Zwart S., Possenti A., 2005, ApJ, 618,845Harris W. E., 1996, AJ, 112, 1487Heggie D. C., Hut P., McMillan S. L. W., 1996, ApJ, 467, 359Hénon M. H., 1971, Ap&SS, 14, 151Hurley J. R., Tout C. A., Pols O. R., 2002, MNRAS, 329, 897Hut P., Bahcall J. N., 1983, ApJ, 268, 319 Ivanova N., Heinke C. O., Rasio F. A., Belczynski K., Fregeau J. M., 2008,MNRAS, 386, 553Ivanova N., Chaichenets S., Fregeau J., Heinke C. O., Lombardi Jr. J. C.,Woods T. E., 2010, ApJ, 717, 948Janka H.-T., 2013, MNRAS, 434, 1355Janka H.-T., 2017, ApJ, 837, 84Jeans J. H., 1919, MNRAS, 79, 408Jonker P. G., Nelemans G., 2004, MNRAS, 354, 355Kharchenko N. V., Piskunov A. E., Schilbach E., Röser S., Scholz R.-D.,2013, A&A, 558, A53Khargharia J., Froning C. S., Robinson E. L., Gelino D. M., 2013, AJ, 145,21King I. R., 1966, AJ, 71, 64Krauss L. M., Chaboyer B., 2003, Science, 299, 65Kroupa P., 2001, MNRAS, 322, 231Kulkarni S. R., Hut P., McMillan S., 1993, Nature, 364, 421Kuulkers E., et al., 2013, A&A, 552, A32LIGO Scientific Collaboration et al., 2015, Class. Quantum Grav., 32,074001Lamberts A., Garrison-Kimmel S., Clausen D. R., Hopkins P. F., 2016,MNRAS, 463, L31Lee H. M., Ostriker J. P., 1986, ApJ, 310, 176Lyne A. G., Lorimer D. R., 1994, Nature, 369, 127Mandel I., 2016, MNRAS, 456, 578Mata Sánchez D., Muñoz-Darias T., Casares J., Corral-Santana J. M., Shah-baz T., 2015, MNRAS, 454, 2199McClintock J. E., Remillard R. A., 2006, Black hole binaries. pp 157–213McLaughlin D. E., 2000, ApJ, 539, 618McLaughlin D. E., Pudritz R. E., 1996, ApJ, 457, 578Meylan G., Heggie D. C., 1997, A&ARv, 8, 1Milone A. P., et al., 2012, A&A, 540, A16Mirabel I. F., Dhawan V., Mignani R. P., Rodrigues I., Guglielmetti F., 2001,Nature, 413, 139Moreno E., Pichardo B., Velázquez H., 2014, ApJ, 793, 110Morozova V., Piro A. L., Renzo M., Ott C. D., Clausen D., Couch S. M.,Ellis J., Roberts L. F., 2015, ApJ, 814, 63Morscher M., Umbreit S., Farr W. M., Rasio F. A., 2013, ApJ, 763, L15Morscher M., Pattabiraman B., Rodriguez C., Rasio F. A., Umbreit S., 2015,ApJ, 800, 9Neustroev V. V., Veledina A., Poutanen J., Zharikov S. V., Tsygankov S. S.,Sjoberg G., Kajava J. J. E., 2014, MNRAS, 445, 2424O’Leary R. M., Rasio F. A., Fregeau J. M., Ivanova N., O’Shaughnessy R.,2006, ApJ, 637, 937Özel F., Psaltis D., Narayan R., McClintock J. E., 2010, ApJ, 725, 1918Paxton B., Bildsten L., Dotter A., Herwig F., Lesaffre P., Timmes F., 2011,ApJS, 192, 3Perets H. B., Li Z., Lombardi Jr. J. C., Milcarek Jr. S. R., 2016, ApJ, 823,113Peters P. C., 1964, Physical Review, 136, 1224Peterson C. J., King I. R., 1975, AJ, 80, 427Pfahl E., Rappaport S., Podsiadlowski P., 2002, ApJ, 573, 283Repetto S., Nelemans G., 2015, MNRAS, 453, 3341Repetto S., Davies M. B., Sigurdsson S., 2012, MNRAS, 425, 2799Rodriguez C. L., Morscher M., Pattabiraman B., Chatterjee S., Haster C.-J.,Rasio F. A., 2015, Phys. Rev. Lett., 115, 051101Rodriguez C. L., Chatterjee S., Rasio F. A., 2016a, Phys. Rev. D, 93, 084029Rodriguez C. L., Morscher M., Wang L., Chatterjee S., Rasio F. A., SpurzemR., 2016b, MNRAS, 463, 2109Sadowski A., Belczynski K., Bulik T., Ivanova N., Rasio F. A.,O’Shaughnessy R., 2008, ApJ, 676, 1162Salpeter E. E., 1955, ApJ, 121, 161Sana H., et al., 2012, Science, 337, 444Shaw A. W., Charles P. A., Casares J., Hernández Santisteban J. V., 2016,MNRAS, 463, 1314Sigurdsson S., Hernquist L., 1993, Nature, 364, 423Sigurdsson S., Phinney E. S., 1993, ApJ, 415, 631Sigurdsson S., Phinney E. S., 1995, ApJS, 99, 609Sippel A. C., Hurley J. R., 2013, MNRAS, 430, L30MNRAS , 1–27 (2017)
H-LMXBs from BH retaining GCs Spitzer L., 1987, Dynamical Evolution of Globular Clusters. Princeton Uni-versity Press, Princeton, NJSpitzer Jr. L., Hart M. H., 1971, ApJ, 164, 399Strader J., Chomiuk L., Maccarone T. J., Miller-Jones J. C. A., Seth A. C.,2012, Nature, 490, 71Torres M. A. P., Callanan P. J., Garcia M. R., Zhao P., Laycock S., KongA. K. H., 2004, ApJ, 612, 1026Wang C., Jia K., Li X.-D., 2016, MNRAS, 457, 1015Webb N. A., Naylor T., Ioannou Z., Charles P. A., Shahbaz T., 2000, MNRAS,317, 528Yamaoka K., et al., 2012, PASJ, 64, 32Zonoozi A. H., Küpper A. H. W., Baumgardt H., Haghi H., Kroupa P., HilkerM., 2011, MNRAS, 411, 1989Zurita C., Durant M., Torres M. A. P., Shahbaz T., Casares J., Steeghs D.,2008, ApJ, 681, 1458de Mink S. E., Sana H., Langer N., Izzard R. G., Schneider F. R. N., 2014,ApJ, 782, 7This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000