Probing the Black Hole Merger History in Clusters using Stellar Tidal Disruptions

TL;DR

This paper proposes using electromagnetic signals from stellar tidal disruptions caused by binary black holes in star clusters as a new method to trace their merger history, offering a novel multi-messenger approach.

Contribution

It introduces a new observable, the variability in tidal disruption events, to distinguish black hole merger pathways in dense star clusters.

Findings

Orbital period distribution peaks on days in globular clusters.

Tidal disruption variability can map BBH orbital periods.

Method can differentiate formation scenarios of BBH mergers.

Abstract

The dynamical assembly of binary black holes (BBHs) in dense star clusters (SCs) is one of the most promising pathways for producing observable gravitational wave (GW) sources, however several other formation scenarios likely operate as well. One of the current outstanding questions is how these different pathways may be distinguished apart. In this paper we suggest a new multi-messenger observable that can be used to constrain the formation of BBH mergers originating from SCs: the electromagnetic signal from tidal disruptions (TDs) of stars by BBHs. Such TDs will show variability in their light curve from the orbital motion of the disruptive BBHs, and can therefore be used to map the BBH orbital period distribution, and thereby also the dynamical mechanisms that eventually drive the BBHs to merger. Using an analytical approach including General Relativistic effects, we find that the orbital period distribution of BBHs within globular clusters peaks on timescales of days, which we argue is unique to this assembly pathway. We propose that the search for variable TDs in current and future EM transient surveys might be used to constrain the merger history of BBHs in SCs.

Equations

$$a_{\mathrm{HB}}=\frac{3}{2}\frac{G{m}_{\mathrm{BH}}}{v_{\mathrm{dis}}^2},$$

$$a_{\mathrm{n}}=a_{\mathrm{HB}}\times {\delta }^n.$$

$$\frac{d{\Gamma }_{\mathrm{TD}}}{d\log T}=\frac{d{N}_{\mathrm{BBH}}}{dn}\frac{d{\Gamma }_{\mathrm{TD}}}{d{N}_{\mathrm{BBH}}}\frac{dn}{d\log T},$$

$$\frac{d{N}_{\mathrm{BBH}}}{dn}\propto {\mathcal{W}}_n\times \mathcal{P}(n),$$

$${\mathcal{W}}_n\propto t_{\mathrm{bs},\mathrm{n}}\approx {\left({\eta }_{\mathrm{BH}}{\sigma }_{\mathrm{bs},\mathrm{n}}{v}_{\mathrm{dis}}\right)}^{-1}\propto {\delta }^{-n},$$

$$\mathcal{P}(n)=\prod _{i=0}^n{\mathcal{P}}_i=\prod _{i=0}^n(1-\widetilde{{\mathcal{P}}_i}),$$

$$\ln \mathcal{P}(n)=\sum _{i=0}^n\ln (1-\widetilde{{\mathcal{P}}_i})\approx -\sum _{i=0}^n\widetilde{{\mathcal{P}}_i},$$

$$\widetilde{{\mathcal{P}}_i}\approx {\left(t_{\mathrm{bs},i}/t_{\mathrm{cm},i}\right)}^{2/7}={\tilde{\mathcal{P}}}_{\mathrm{HB}}\times {\epsilon }^i,$$

$${\tilde{\mathcal{P}}}_{\mathrm{HB}}\approx {\left(\frac{8192}{3645\pi }\frac{1}{G^3{c}^5}\frac{v_{\mathrm{dis}}^{11}}{m_{\mathrm{BH}}^3{\eta }_{\mathrm{BH}}}\right)}^{2/7},$$

$$\ln \mathcal{P}(n)\approx -{\tilde{\mathcal{P}}}_{\mathrm{HB}}\sum _{i=0}^n{\epsilon }^i.$$

$$\mathcal{P}(n)\approx \exp \left(-{\tilde{\mathcal{P}}}_{\mathrm{HB}}\times \frac{1-{\epsilon }^{n+1}}{1-\epsilon }\right).$$

$$n_{\max }\approx -\frac{\log {\tilde{\mathcal{P}}}_{\mathrm{HB}}}{\log \epsilon }.$$

$$\frac{d{\Gamma }_{\mathrm{TD}}}{d{N}_{\mathrm{BBH}}}\approx \frac{d{\Gamma }_{\mathrm{s}}}{d{N}_{\mathrm{BBH}}}\times P_{\mathrm{TD}},$$

$$\frac{d{\Gamma }_{\mathrm{TD}}}{d{N}_{\mathrm{BBH}}}\propto {\eta }_{\mathrm{s}}{\sigma }_{\mathrm{s}}{v}_{\mathrm{dis}}\times \frac{R_{\mathrm{TD}}}{a}\propto \frac{{\eta }_{\mathrm{s}}{m}_{\mathrm{BH}}{R}_{\mathrm{TD}}}{v_{\mathrm{dis}}}.$$

$$\begin{array}{rl}\frac{d{\Gamma }_{\mathrm{TD}}}{d\log T}& \propto \frac{d{N}_{\mathrm{BBH}}}{dn}\propto \frac{d{N}_{\mathrm{BBH}}}{d\log T}\\ & \propto T^{-2/3}\exp \left(-{\tilde{\mathcal{P}}}_{\mathrm{HB}}\times \frac{1-\epsilon (T/T_{\mathrm{HB}}{)}^{-20/21}}{1-\epsilon }\right),\end{array}$$

$$T_{\mathrm{HB}}=\frac{\pi \sqrt{27}}{2}\frac{G{m}_{\mathrm{BH}}}{v_{\mathrm{dis}}^3},$$

$$\begin{array}{rl}{T}_{\min ,\tau }& \approx T_{\mathrm{HB}}\times {\left(t_{\mathrm{bs},\mathrm{HB}}/\tau \right)}^{3/2}{\left(\frac{\delta }{1-\delta }\right)}^{3/2}\\ & \propto m_{\mathrm{BH}}^{-2}{v}_{\mathrm{dis}}^{3/2}{\eta }_{\mathrm{BH}}^{-3/2},\end{array}$$

$$\begin{array}{rl}{T}_{\min ,\mathrm{m}}& \approx T_{\mathrm{HB}}\times {\tilde{\mathcal{P}}}_{\mathrm{HB}}^{21/20}\\ & \propto m_{\mathrm{BH}}^{1/10}{v}_{\mathrm{dis}}^{3/10}{\eta }_{\mathrm{BH}}^{-3/10},\end{array}$$

$$\begin{array}{rl}{T}_{\min ,\mathrm{ej}}& \approx T_{\mathrm{HB}}\times (v_{\mathrm{esc}}/v_{\mathrm{dis}}{)}^{-3}{\left(\frac{1-\delta }{9\delta }\right)}^{3/2}\\ & \propto m_{\mathrm{BH}}{v}_{\mathrm{dis}}^{-3},\end{array}$$

$$\begin{array}{rl}{T}_{\mathrm{peak}}& =T_{\mathrm{HB}}\times {\tilde{\mathcal{P}}}_{\mathrm{HB}}^{21/20}{\left(\frac{10}{7}\frac{\epsilon }{\epsilon -1}\right)}^{21/20}\\ & \propto m_{\mathrm{BH}}^{1/10}{v}_{\mathrm{dis}}^{3/10}{\eta }_{\mathrm{BH}}^{-3/10}\end{array}$$

$$\frac{d{N}_{\mathrm{BBH}}}{d\log T}\propto T^{-2/3}\mathcal{P}(T)\text{(dynamics)},$$

$$\frac{d{N}_{\mathrm{BBH}}}{d\log T}\propto \text{const.}\text{(Opik’s Law)}.$$

$$\frac{d{N}_{\mathrm{BBH}}}{d\log T}\propto T^{8/3}\text{(e=0, GW driven)}.$$

$$\frac{{\Gamma }_{\mathrm{TD}}^{\mathrm{BBH}}}{\text{gal.}}\approx 10^{-7}{\text{yr}}^{-1}\left(\frac{{\eta }_{\mathrm{s}}}{10^4{\text{pc}}^{-3}}\right){\left(\frac{m_{\mathrm{BH}}}{30M_{\odot }}\right)}^{4/3}\left(\frac{12{\text{kms}}^{-1}}{v_{\mathrm{dis}}}\right),$$

Full text

I Introduction

During the first and second observing runs of the advanced gravitational-wave detector network, LIGO and VIRGO detected ten merging binary black holes (BBHs) (2016PhRvL.116f1102A, ; 2016PhRvL.116x1103A, ; 2016PhRvX…6d1015A, ; 2017PhRvL.118v1101A, ; 2017PhRvL.119n1101A, ; 2018arXiv181112907T, ). Upgraded detectors are expected to detect hundreds of such mergers within the next few years (see e.g., 2016ApJ…833L…1A, ). The origin and distribution of these BBHs is a key scientific question that remains to be addressed (e.g., 2018arXiv180605195B, ).

Several formation channels have been suggested in the literature: dense stellar clusters (2000ApJ…528L..17P, ; 2010MNRAS.402..371B, ; 2013MNRAS.435.1358T, ; 2014MNRAS.440.2714B, ; 2015PhRvL.115e1101R, ; 2016PhRvD..93h4029R, ; 2016ApJ…824L…8R, ; 2016ApJ…824L…8R, ; 2017MNRAS.464L..36A, ; 2017MNRAS.469.4665P, ), field binaries (2012ApJ…759…52D, ; 2013ApJ…779…72D, ; 2015ApJ…806..263D, ; 2016ApJ…819..108B, ; 2016Natur.534..512B, ; 2017ApJ…836…39S, ; 2017ApJ…845..173M, ; 2018ApJ…863….7R, ; 2018ApJ…862L…3S, ), active galactic nuclei discs (2017ApJ…835..165B, ; 2017MNRAS.464..946S, ; 2017arXiv170207818M, ), galactic nuclei (2009MNRAS.395.2127O, ; 2015MNRAS.448..754H, ; 2016ApJ…828…77V, ; 2016ApJ…831..187A, ; 2016MNRAS.460.3494S, ; 2017arXiv170609896H, ), single-single GW captures of primordial black holes (BHs) (2016PhRvL.116t1301B, ; 2016PhRvD..94h4013C, ; 2016PhRvL.117f1101S, ; 2016PhRvD..94h3504C, ), and very massive stellar mergers (Loeb:2016, ; Woosley:2016, ; Janiuk+2017, ; DOrazioLoeb:2017, ). There are significant uncertainties in the predicted merger rates in these scenarios due to the complexity of the underlying astrophysical environments. The mergers observed by LIGO are in the gravitational-wave dominated regime, and are hence described only by a few intrinsic parameters, such as the masses and the spins of the black holes. Given a population of detected mergers, the distributions of these intrinsic parameters can be used to disentangle these channels.

Binary black holes formed in dynamical environments, such as globular clusters (GCs), galactic nuclei (GN), and nuclei star clusters (NSCs), can be eccentric in both LIGO (2009MNRAS.395.2127O, ; Kocsis:2012ja, ; 2014ApJ…784…71S, ; 2016ApJ…824L..12O, ; 2017ApJ…840L..14S, ; 2018MNRAS.476.1548S, ; 2018ApJ…853..140S, ; 2019MNRAS.482…30S, ; 2018PhRvD..97j3014S, ; 2018ApJ…855..124S, ; 2018arXiv181000901Z, ; 2018arXiv181104926R, ; 2018ApJ…860….5G, ) and LISA bands (e.g., 2018MNRAS.tmp.2223S, ; 2018arXiv181111812K, ), and have randomly oriented spins when exchange interactions are effective (e.g., 2016ApJ…832L…2R, ). In comparison, isolated field binaries are likely to have correlated spins and merge on near circular orbits (e.g. 2018ApJ…862L…3S, ). Exceptions to this standard picture include field BBHs undergoing Lidov-Kozai (LK, (e.g., 1962AJ…..67..591K, )) oscillations that tend to re-orient spins (e.g., 2017ApJ…846L..11L, ; 2018MNRAS.480L..58A, ; 2018ApJ…863…68L, ) and produce eccentric mergers (e.g., 2017ApJ…836…39S, ; 2018arXiv180506458A, ), and dynamically formed BBHs that have their spins re-oriented and possibly aligned through mass accretion following tidal disruptions in dense clusters (e.g. 2018arXiv181201118L, ). Binaries of primordial BHs that are formed via GW-driven single-single capture are also expected to show signs of eccentricity (e.g., 2016PhRvD..94h4013C, ). As a non-zero eccentricity would indicate a dynamical origin of BBH mergers, major effort is therefore being undertaken to observationally probe eccentric mergers (e.g., 2017PhRvD..95b4038H, ; 2018PhRvD..97b4031H, ; 2018arXiv180502716H, ; 2018arXiv180900672G, ). Various methods have been suggested for distinguishing between dynamical channels, including correlating with the BBH masses, cosmic merger evolution (e.g., 2018ApJ…866L…5R, ), and using future multi-band GW observations (e.g., 2017ApJ…842L…2C, ; 2018MNRAS.tmp.2223S, ; 2018MNRAS.481.4775D, ).

In this paper, we discuss a new electromagnetic window into the assembly of BBH mergers in dense stellar clusters (SCs): tidal disruptions (TDs) of stars by BBHs. Previous theoretical work on these events suggests that their lightcurves are interrupted on characteristic timescales related to the orbital period of the disrupting BBHs (e.g., Liu:2009fl, ; 2014ApJ…786..103L, ; 2016MNRAS.458.1712R, ; 2017MNRAS.465.3840C, ; 2018arXiv181201118L, ). We propose that populations of such events can be used to probe the orbital period distribution of BBHs inside SCs. We explore in detail how the dynamics that drive the BBHs to merger (2006ApJ…637..937O, ; 2018PhRvL.120o1101R, ; 2018PhRvD..97j3014S, ; 2018MNRAS.tmp.2223S, ; 2018arXiv181104926R, ) also shape their orbital period distribution. We use these insights to derive, for the first time, a simple analytical model of the BBH period distribution with the inclusion of General Relativistic (GR) effects, and relate the TD and BBH merger rates from SCs.

We find that the orbital period distribution of BBHs within standard GCs peaks on timescales of days, with only a weak dependence on BH mass. This period distribution is entirely different from other theoretical distributions derived under simple considerations, such as those driven purely by GW emission (2017MNRAS.469..930C, ), and the empirical Opik’s Law distribution for field binaries 1924PTarO..25f…1O . This emboldens us to suggest that if current and future transient surveys observe a population of interrupted TDs, the inferred orbital periods will directly inform us about the dynamical assembly of BBHs in dense SCs. Provided we properly account for the observational biases of the surveys and understand the factors governing the detectability of the TDs, such detections can also constrain the contribution of SCs to the BBH merger rate in the universe.

The paper is structured as follows. In Section II we describe our GR dynamical model and derive an analytical expression for the rate distribution of star-BBH TDs as a function of the BBH orbital time. This is followed by Section III, where we discuss different astrophysical limits on the BBH orbital distribution. In Section IV we summarize our main results from our dynamical models, where in Section V we compare to numerical data derived using Monte Carlo (MC) techniques. In Section VI we discuss prospects and challenges related to using interrupted TDs to probe the BBH orbital distribution. We conclude our study in Section VII.

II Stellar Disruptions by Binary Black Holes

In this section we derive the rate of stars disrupted by BBHs in SCs per BBH orbital period, as a function of the BH masses and SC properties. For this, we use our analytical model presented in 2018PhRvD..97j3014S ; 2018MNRAS.tmp.2223S for describing the GR dynamical evolution of the BBHs, that we throughout this work assume are equal mass. This is a reasonable approximation as both mass-segregation and few-body dynamics tend to keep equal mass objects together (e.g., Heggie:1975uy, ; 2016PhRvD..93h4029R, ; 2017MNRAS.469.4665P, ). Our derivation is divided into a few separate steps, as indicated by the following subsections. Below we start by describing our dynamical model.

II.1 Dynamical Model

In our model we assume that BBHs form at a steady rate in the SC core, through the interaction of three initially unbound BHs (e.g., Heggie:1975uy, ; 1976A&A….53..259A, ; Hut:1983js, ), with their semi-major axes (SMA) equal to the hard-binary (HB) value (Hut:1983js, ),

[TABLE]

where $m_{\rm BH}$ denotes the mass of one of the (assumed equal mass) BHs, and $v_{\rm dis}$ is the SC velocity dispersion. After formation the dynamics of the BBHs is governed by binary-single interactions involving other (equal mass) BHs. Each of these interactions leads to a decrease in the SMA of the interacting BBH from SMA $a$ to $\delta a$, a process referred to as ‘hardening’, where we assume a constant value $\delta=7/9$, which is the mean value for equal mass interactions (2018PhRvD..97j3014S, ). Following this approximation, a given BBH will therefore after $n$ interactions have a SMA equal to

[TABLE]

In this paper we will occasionally refer to a BBH with SMA $a_{\rm n}$ to be in ‘state $n$’ or at ‘hardening step $n$’. This binary-single hardening process continues until either the release of binding energy from a binary-single interaction recoils the BBH out the cluster (e.g., 2000ApJ…528L..17P, ), or that the interacting BBH undergoes a GW merger during (e.g., 2006ApJ…640..156G, ; 2014ApJ…784…71S, ; 2018arXiv181000901Z, ; 2018arXiv181104926R, ) or in-between (e.g., 2018PhRvL.120o1101R, ; 2018PhRvD..97j3014S, ; 2018arXiv181104926R, ) its binary-single interactions. For further descriptions of this model we refer the reader to 2018MNRAS.tmp.2223S . Information is also provided in Figure 1.

II.2 Formulation

We start by considering an ensemble of SCs each with a constant production rate of BBHs that undergo binary-single hardening according to our model described in the above Section II.1. We consider this distribution of BBHs at a random point in time – the time of observation – at which we observe a set of stars disrupted by the BBHs. The question is, what is the distribution of BBH orbital periods for the set of BBHs disrupting the stars? For answering this, we start by writing the differential tidal disruption rate per $\log$ BBH orbital time as a product of the following terms,

[TABLE]

where $T$ is the BBH orbital period, $\Gamma_{\rm TD}$ is the rate of stars disrupted by BBHs, $n$ is the $n$’th hardening step, and $N_{\rm BBH}$ is the number of BBHs present at state $n$ at the time of observation. The last term ${dn}/{d\log{T}}$ is a simple Jacobian factor that can be calculated from first using that $T_{\rm n}\propto a_{\rm n}^{3/2}\propto\delta^{3n/2}$, where $T_{\rm n}$ denotes the BBH orbital period at state $n$, which follows from Kepler’s laws and Eq. (2). This implies that $\log T_{\rm n}\propto n(3/2)\log\delta\propto n$, from which we conclude that ${dn}/{d\log{T}}$ simply equals a constant. In the following subsections we derive the two remaining terms ${dN_{\rm BBH}}/{dn}$ and ${d\Gamma_{\rm TD}}/{dN_{\rm BBH}}$, where our final expression for ${d\Gamma_{\rm TD}}/{d\log{T}}$ is presented in Section II.5.

II.3 Distribution of Binary Black Holes

For calculating the term ${dN_{\rm BBH}}/{dn}$, i.e. the number of BBHs per state $n$ at the time of observation, we start by factorizing it as,

[TABLE]

where $\mathscr{W}_{n}$ is a weight factor that essentially equals the probability that the BBH is in state $n$ at the time of observation given that it is not able to undergo any type of merger inside the cluster (Newtonian term), and $\mathscr{P}(n)$ is the probability that the BBH do not merge before and during state $n$ (GR term). Note that we do not consider the merging population in our calculations, i.e. we only account for the BBHs in Eq. (4) that also ‘survives’ state $n$, although there is a small probability for the merging BBHs to disrupt stars during their GW inspiral: This represents an interesting case, but is beyond what we can study with our semi-analytical models. In the following two sections we derive the ‘observational weight factor’ $\mathscr{W}_{n}$ and the term describing the ‘depletion from GW Mergers’ $\mathscr{P}(n)$, respectively.

II.3.1 Observational Weight Factor

The factor $\mathscr{W}_{n}$ is proportional to the time the BBH spends at state $n$, which in our model equals the time between binary-single interactions at state $n$, a time we denote $t_{\rm bs,n}$. The factor $\mathscr{W}_{n}$ can therefore be approximated by (e.g., 2018arXiv180708864S, ),

[TABLE]

where $\eta_{\rm BH}$ is the number density of single BHs in the SC core, and $\sigma_{\rm bs,n}\propto m_{\rm BH}a_{\rm n}/v_{\rm dis}^{2}$ is the cross section for a strong binary-single interaction at state $n$ (e.g., 2018ApJ…853..140S, ).

II.3.2 Depletion from Gravitational Wave Mergers

To derive $\mathscr{P}(n)$ we first express it as the following product,

[TABLE]

where $\mathscr{P}_{i}$ denotes the probability for that the BBH do not merge at state $i$, and $\tilde{\mathscr{P}_{i}}$ that the BBH do merge at state $i$. To evaluate this product we first rewrite it as,

[TABLE]

where for the last term we have assumed that $\tilde{\mathscr{P}_{i}}\ll 1$. The term $\tilde{\mathscr{P}_{i}}$ denotes the total probability for the BBH to merge at state $i$, which in our model either can happen in the time-span before it undergoes its next binary-single interaction (e.g., 2018MNRAS.tmp.2223S, ) or during the binary-single interaction itself through a two-body GW capture (2014ApJ…784…71S, ). Normally, the probability for undergoing a merger in-between binary-single interactions is far larger than doing a binary-single interaction, therefore, to leading order $\tilde{\mathscr{P}_{i}}$ is well approximated by the probability that the BBH at state $i$ has a GW inspiral life time, $t_{{\rm m},i}$, less than its binary-single encounter time, $t_{{\rm bs},i}$. Using that the GW inspiral life time can be written as $t_{{\rm m},i}\approx t_{{\rm cm},i}\times(1-e^{2})^{7/2}$, where $t_{{\rm cm},i}$ denotes the circular GW inspiral life time for which the BBH eccentricity $e=0$ (Peters:1964bc, ), and assuming that the BBH eccentricity distribution follows that of a so-called thermal distribution $P(e)=2e$ (Heggie:1975uy, ), the probability $\tilde{\mathscr{P}_{i}}$ can then be written as (e.g., 2018ApJ…853..140S, ; 2018PhRvD..97j3014S, ),

[TABLE]

where $\varepsilon\equiv\delta^{-10/7}$, and $\tilde{\mathscr{P}}_{\rm HB}$ is the probability that the BBH merges before its next binary-single interaction at the HB limit when its SMA $=a_{\rm HB}$. This probability can be expressed as,

[TABLE]

which follows directly from Eq. (5), Eq. (8), and the relation $t_{{\rm cm},i}\propto a_{i}^{4}/m_{\rm BH}^{3}$ given by Peters:1964bc . Now using the relation from the above Eq. (8) we can write Eq. (7) as,

[TABLE]

This sum can be evaluated analytically using the solution $\sum_{i=0}^{n}x^{i}=(1-x^{n+1})/(1-x)$, which follows from the standard literature on ‘geometrical series’. With this, we finally find,

[TABLE]

As seen, the effect from including the possibility for a BBH to merge in-between its binary-single interactions gives rise to a BBH depletion that takes the form of a super-exponential decay. We note here that the BBHs that merge in-between their binary-single interactions, i.e. the population giving rise to the decay, have been shown to constitute $\approx 30-50\%$ of the observable BBH merger rate from GCs (e.g., 2018PhRvL.120o1101R, ; 2018MNRAS.tmp.2223S, ; 2018arXiv181104926R, ); the function $\mathscr{P}(n)$ can therefore be directly linked to LIGO/Virgo events. One should note that the above expressions are only valid for $n<n_{\rm max}$, where $n_{\rm max}$ is the value of $n$ for which the corresponding $\tilde{\mathscr{P}_{n}}=1$. Using Eq. (8) we find,

[TABLE]

Fig. 2 illustrates the terms $\mathscr{W}_{n}$, $\mathscr{P}(n)$, and the factor ${dN_{\rm BBH}}/{dn}\propto\mathscr{W}_{n}\times\mathscr{P}(n)$ given by Eq. (4), as a function of $n$ for values describing a dense SC system (see figure caption). The maximum value $n_{\rm max}$ for these values is $n_{\rm max}\approx 31$, which indeed is where $\mathscr{P}(n)$ starts to rapidly fall off. The limit for $n_{\rm max}$ is not naturally build into the equations, as we assumed in Eq. (7) that $\tilde{\mathscr{P}_{i}}\ll 1$, which helped us evaluating the sum. Our main results are not affected by these small differences. Other upper limits on $n$ exists as well, some of which will be described in Section III.

Finally, we note that in 2018PhRvD..97j3014S the total probability for a BBH to merge in-between its binary-single interactions during hardening from the HB limit to step $n$ was derived in the integral-limit using slightly different approximations (Section III.D.2 in 2018PhRvD..97j3014S ), from which it was found that $\tilde{\mathscr{P}}(n)\approx(7/10)\tilde{\mathscr{P}}_{n}/(1-\delta)$. In certain limits our derived expression from Eq. (11) should be identical to that previous result. After some mathematical manipulations, one indeed finds that $\tilde{\mathscr{P}}(n)=1-\exp(-\tilde{\mathscr{P}}_{\rm HB}{(1-{\varepsilon}^{n+1})}/{(1-{\varepsilon})})$ is $=(7/10)\tilde{\mathscr{P}}_{n}/(1-\delta)$ in the limit where $\mathscr{P}(n)\gg 0$, $\tilde{\mathscr{P}}_{n}\gg\tilde{\mathscr{P}}_{\rm HB}$, and $\delta\approx 1$, all of which were assumed in 2018PhRvD..97j3014S . This serves as an excellent check of the different approaches that so far have been employed to estimate the probability for a BBH to undergo a GW merger during hardening.

II.4 Rate of Tidal Disruptions

To estimate the term ${d\Gamma_{\rm TD}}/{dN_{\rm BBH}}$, i.e. the stellar tidal disruption rate per BBH, we first factorize it as,

[TABLE]

where $\Gamma_{\rm s}$ denotes the rate of star-BBH interactions, and $P_{\rm TD}$ the probability for that a single star-BBH interaction results in a stellar tidal disruption. The first term $d\Gamma_{\rm s}/dN_{\rm BBH}$ is $\propto\eta_{\rm s}\sigma_{\rm s}v_{\rm dis}$, where $\eta_{\rm s}$ is the number density of stars, and $\sigma_{\rm s}\propto m_{\rm BH}a/v_{\rm dis}^{2}$ is the cross section for a star-BBH interaction (e.g., 2017ApJ…846…36S, ). The second term $P_{\rm TD}$ is proportional to the probability for that the incoming star undergoes a pericenter passage w.r.t. to one of the two BHs during the star-BBH interaction that is $\lesssim R_{\rm TD}$, where $R_{\rm TD}$ is the stellar disruption distance. In the equal mass case estimating this probability is straight forward and can e.g. be done by describing the temporary three-body state as a binary with a bound single (e.g., 2018PhRvD..97j3014S, ); however, the relative small mass of the star in our case will rarely allow it to form a bound system with any of the BHs. Simulations show that the light star instead will undergo chaotic single encounters with each BH on orbits with a wide range of energies and eccentricities (e.g., 2017MNRAS.465.3840C, ; 2018MNRAS.477.4009D, ; 2018ApJ…863L..24C, ; 2018arXiv181201118L, ). Especially, one finds that the star is highly likely to disrupt on orbits that, relative to the disrupting BH, are unbound, which is possible since each BH evolves with a relative high velocity around their common center-of-mass (e.g. 2018arXiv181201118L, ). Although one might be able to develop an analytical description of this limit, we here make use of the recent work by 2018MNRAS.477.4009D , which indicates the probability $P_{\rm TD}$ is $\propto R_{\rm TD}/a$ (we note that this exact same scaling is found in the equal mass case (e.g., 2017ApJ…846…36S, ; 2018PhRvD..97j3014S, )). Using this scaling one now finds,

[TABLE]

From this one concludes that ${d\Gamma_{\rm TD}}/{dN_{\rm BBH}}$ is independent of the BBH orbital period $T$ for fixed $m_{\rm BH}$. The term ${d\Gamma_{\rm TD}}/{dN_{\rm BBH}}$ acts therefore only as a normalization factor, and plays therefore no role in the calculation of the functional form for ${d\Gamma_{\rm TD}}/{d\log T}$. Finally, we note that one do expect corrections to this simple picture when the stellar radius becomes comparable to the BBH SMA $a$ (e.g. 2018arXiv181201118L, ).

II.5 Final Solution

Using that $\delta^{n}\propto T_{\rm n}^{2/3}$ one can now finally write the solution to Eq. (3) in terms of $T$ as,

[TABLE]

where $T_{\rm HB}$ denotes the BBH orbital period at the HB limit,

[TABLE]

which follows from Eq. (1). Note here that in the above Eq. (15) we have further indicated that the derived distribution is proportional to the term ${dN_{\rm BBH}}/{d\log T}$. This is a highly useful quantity that will be relevant for other and future studies. One should also note that the above Eq. (15) does not include possible lower limits on $T$, but this will be discussed in further detail in Section III, before we explore our solution in Section IV.

Finally, to check our derivations, we compared our expression from Eq. (15) with results derived using our semi-analytical MC code described in 2018MNRAS.tmp.2223S . We find excellent agreement, which serve as a good validation of our results so far. In Section V we show comparison of our models to data derived using a Hénon-type MC code.

III Limits on the Binary Black Hole Orbital Period

The upper limit on the BBH orbital period $T$, denoted by $T_{\rm max}$, is in our model always set by the HB value given by Eq. (16); however, different dynamical and astrophysical effects will determine the lower limit, $T_{\rm min}$. In the following we discuss the three limits: ‘Time Limit’ (${T_{\rm min,\tau}}$), ‘Merger Limit’ ($T_{\rm min,m}$), and ‘Ejection Limit’ (${T_{\rm min,ej}}$), so that $T_{\rm min}=\max([T_{\rm min,\tau},T_{\rm min,m},T_{\rm min,ej}])$. In the sections below we quote numbers based on the following ‘fiducial values’: $m_{\rm BH}=30M_{\odot}$, $v_{\rm dis}=12\ \text{kms}^{-1}$, $v_{\rm esc}/v_{\rm dis}=5$, and $\eta_{\rm BH}=10^{4}\ \text{pc}^{-3}$.

III.0.1 Time Limit

It takes a finite time $\tau$ for a BBH to harden from its initial orbital time $T_{\rm HB}$ to some smaller time $T_{\tau}$. An imposed limit on $\tau$ will therefore map to a corresponding lower value on $T$, that we here denote $T_{\rm min,\tau}$. The absolute upper limit on $\tau$ is the Hubble time $t_{\rm H}$. For deriving $T_{\rm min,\tau}$, we start by writing out the time it takes for a BBH to reach state $n$, $\tau_{\rm n}=\sum_{i=0}^{i=n-1}t_{\rm bs,i}=t_{\rm bs,HB}(1-(T_{\rm n}/T_{\rm HB})^{-2/3})/(1-1/\delta)$, where we have used that $t_{{\rm bs},i}=t_{\rm bs,HB}\times\delta^{-i}$ which follows from Eq. (5), and $t_{\rm bs,HB}$ refers to the time between binary-single interactions at the HB limit. This expression can be rearranged from which we now find,

[TABLE]

where $T_{\rm min,\tau}\approx 0.3\ \text{days}$ for our chosen fiducial values and $\tau=t_{\rm H}$. Note here the strong dependence on $m_{\rm BH}$.

III.0.2 Merger Limit

In our model it is impossible for a BBH to pass state $n_{\rm max}$ given by Eq. (12) through binary-single interactions alone. Since $T_{\rm n}\propto(9/7)^{-3n/2}$ then the upper limit $n_{\rm max}$ directly maps to a lower limit on $T$, referred to as $T_{\rm min,m}$. This limit can be found by the use of Eq. (12), or by solving for the period $T$ that satisfies $t_{\rm bs}=t_{\rm cm}$, from which follows,

[TABLE]

where $T_{\rm min,m}\approx 1.7\ \text{days}$ for our chosen fiducial values. Note that this limit is not automatically included in Eq. (15) as we assumed that $\tilde{\mathscr{P}_{i}}\ll 1$ in Eq. (7) to evaluate the summation; however, Eq. (15) does rapidly decline just around that limit as we later illustrate in Fig. 3.

III.0.3 Ejection Limit

A BBH undergoing binary-single interactions inside a stellar cluster receives in each interaction (if no merger takes place) a dynamical velocity kick $v_{\rm kick}$ that increases with increasing $n$, i.e. decreasing SMA $a_{\rm n}$ (e.g., 1992ApJ…389..527H, ; 2018PhRvD..97j3014S, ). The BBH in its hardening process will therefore eventually receive a $v_{\rm kick}$ that exceeds the escape velocity $v_{\rm esc}$ of the cluster. As shown in 2018PhRvD..97j3014S , the value of the BBH SMA for which $v_{\rm kick}=v_{\rm esc}$ can be written as $a_{\rm ej}=({1}/{6})\left({1}/{\delta}-1\right){m_{\rm BH}}/{v_{\rm esc}^{2}}$. This marks a lower limit on the BBH SMA, as no BBH is able to pass $a_{\rm ej}$ without being ejected from the cluster (see Fig. 1). The corresponding minimum orbital period $T_{\rm min,ej}$ is from Kepler’s laws given by,

[TABLE]

where $T_{\rm min,ej}\approx 9.8\ \text{days}$ for our chosen fiducial values. Note here the strong dependence on $v_{\rm dis}$. Comparing with the other two considered limits, $T_{\rm min,ej}$ takes the largest value and represents therefore the actual lower limit for our fiducial values. However, our three considered limits scale differently with $m_{\rm BH},v_{\rm dis},\eta_{\rm BH}$, and the two other limits might therefore dominate in other cases. We comment on this below, where we discuss our results including the derived limits.

IV Results

The distribution ${d\Gamma_{\rm TD}}/{d\log{T}}$ given by Eq. (15) is plotted in Fig. 3 as a function of log $T$ for individual BH masses $m_{\rm BH}[M_{\odot}]=5,10,20,30$. The solid lines show our derived solution that includes the possibility for GW merger during hardening, where the dotted lines show the Newtonian solution $\propto T^{-2/3}$ obtained from setting $\tilde{\mathscr{P}}_{\rm HB}=0$ in Eq. (15). As seen, the four distributions all peak near the same orbital time $T$. Furthermore, the peak locations are neither located near the minimum ($T_{\rm min}$) nor the maximum orbital time ($T_{\rm HB}$), but somewhere in between. The explanation follows from the two terms in Eq. (4): BBHs with high values of $T$ are less likely to be observed due to their relative small $t_{\rm bs}$ (Eq. (5)), where the number of BBHs with low values of $T$ is greatly reduced due to depletion through mergers during the hardening (Eq. (11)). This is in contrast to the Newtonian solution, which predicts that the $T$ distributions always will peak near their lowest possible value.

To see how the location of the peak value, $T_{\rm peak}$, scales with $m_{\rm BH}$ and the cluster parameters, we can solve for the peak location using the standard procedure by differentiating ${d\Gamma_{\rm TD}}/{d\log{T}}$ w.r.t. $\log T$, equal the expression to zero, and then solve for $\log T$. From this we find that $T_{\rm peak}$ can be written as,

[TABLE]

where $T_{\rm peak}\approx 8.9\ \text{days}$ for our chosen fiducial values. As seen, the location is indeed almost independent of $m_{\rm BH}$ due to the suppression $1/10$ in the exponent, and only weakly dependent on $v_{\rm dis}$ and $\eta_{\rm BH}$.

In Fig. 3 we have shown with vertical lines the lower limits $T_{\rm min,t_{\rm H}}$, $T_{\rm min,m}$, and $T_{\rm min,ej}$, defined in Section III. As seen, for our chosen values, $T_{\rm min}$ for BH masses $m_{\rm BH}[M_{\odot}]=5,10,20,30$ is $T_{\rm min,\tau}(>T_{\rm peak})$, $T_{\rm min,ej}(<T_{\rm peak})$, $T_{\rm min,ej}(<T_{\rm peak})$, $T_{\rm min,ej}(>T_{\rm peak})$, respectively. It is clear that $T_{\rm min,ej}$ is likely to be the most relevant limit for merging BBHs in classical GCs; however, the value of $T_{\rm min,ej}$ depends strongly on $v_{\rm dis}$ as seen in Eq. (19). This implies, for example, that for our chosen cluster values a dispersion $v_{\rm dis}\gtrsim 20\ \text{kms}^{-1}$ the lower limit $T_{\rm min}$ will instead be $T_{\rm min,m}$ for all our considered values of $m_{\rm BH}$. In that case ${T_{\rm peak}}/{T_{\rm min,m}}=\left(({10}/{7}){{\varepsilon}}/({{\varepsilon}-1})\right)^{21/20}\approx 5$, i.e., the peak $T_{\rm peak}$ is always visible in the high velocity dispersion limit, and no sources are expected to be found with $T\lesssim T_{\rm peak}/5$. This limit further implies and confirms that in NSCs all BBHs will merge inside the system before dynamical ejection takes place (e.g., 1995MNRAS.272..605L, ; 2016ApJ…831..187A, ), which greatly increases the possibility for repeated mergers and possible build-up of heavy central objects (e.g., 2016ApJ…831..187A, ; 2018arXiv181103640A, ).

V Comparing to MC Simulations

In Figure 4 we compare our analytical derived solutions from Section II.5 to the BBH orbital period distribution, $dN_{\rm BBH}/d\log T$, extracted from full-scale GC models developed using the MC code, CMC. CMC is a Hénon-type MC code used to model the long-term evolution of GCs, incorporating various evolutionary effects including two-body relaxation, single and binary-star evolution, and small-$N$ resonant encounters with GR corrections (for a review and summary of the latest updates to CMC, see (Henon1971a, ; Henon1971b, ; Joshi2000, ; Joshi2001, ; Fregeau2003, ; Pattabiraman2013, ; Chatterjee2010, ; Chatterjee2013, ; Rodriguez2015a, ; Rodriguez2018a, ; Rodriguez2018b, )). The data we compare to is taken from the CMC models used in (2018arXiv181111812K, ) and (Rodriguez2018b, ). We consider only those CMC models with $N=2\times 10^{6}$ particles, which most closely match our chosen fiducial cluster parameters ($v_{\rm dis}=12\ \text{kms}^{-1}$, $v_{\rm esc}/v_{\rm dis}=5$, and $\eta_{\rm BH}=10^{4}\ \text{pc}^{-3}$; see Section III). In order to generate a GC population representative of the clusters found in the local universe at present, we use a scheme similar to that of (2018arXiv181111812K, ) where the cluster models are assigned weights based on the present-day GC mass function and metallicity distributions of (El-Badry2018, ). Additionally, cluster ages are drawn from the metallicity-dependent age distributions of (El-Badry2018, ). We perform $10^{4}$ weighted cluster draws from the CMC models and then extract the orbital parameters of all BBHs found in each cluster draw at the present-day in order to generate the CMC orbital period distribution in Figure 4. The two black curves in the figure show our solution given by Eq. (15), with lower cutoffs set by $T_{\rm min,ej}$ derived using Eq. (19), for two different set of cluster parameters (‘params. 1’ and ‘params. 2’), as specified in the caption.

The overlap between the CMC data and our analytical solution is surprisingly good. The distribution $dN_{\rm BBH}/d\log T$ derived from CMC follows exactly our main predictions including an initial rise $\propto T^{-2/3}$ with a lower cutoff set by a combination of dynamical ejections and relativistic BBH merger depletion. As seen, the ‘params. 2’ values seem to fit the data slightly better than our fiducial values. It would be very interesting to explore how a MC derived distribution $dN_{\rm BBH}/d\log T$ changes with a wider range of GC parameters, metallicities, and age; however, that is out of the scope of this paper. With this excellent agreement between numerical and analytical techniques, we proceed below by discussing our results.

VI Discussion

In the previous sections, we derived the distribution of the orbital periods of BBHs in SCs. Here we start by comparing this result to the distributions expected in other astrophysical environments. We finish by discussing the prospects for and challenges involved in using interrupted TDs as probes of the BBH orbital time distribution, and comment on other probes.

VI.1 Orbital Distributions

The orbital period distribution of BBHs in SCs is shaped by the same GR dynamics that also drive the BBHs to merger. A measure of the orbital period distribution can therefore be used to indirectly probe the corresponding BBH merger history. However, this is only possible if the distribution from SCs is different from those found in other scenarios. In the paragraphs below we compare our derived distribution given by Eq. (15),

[TABLE]

to BBHs following Opik’s Law and BBHs that evolve through GW emission only.

The ‘Opik’s Law distribution’ suggests that the distribution of binary SMA $a$ follows $dN_{\rm BBH}/da\propto 1/a$ in the field 1924PTarO..25f…1O . This relation corresponds to an orbital distribution scaling as,

[TABLE]

Opik’s distribution corresponds therefore to a flat distribution in $\log T$. This is clearly different from our derived SC distribution which decreases for increasing $\log T$.

For BBHs evolving by GW radiation only, the steady state solution in the $e=0$ limit can be found from solving a homogeneous advection equation (2016PhRvL.116w1102S, ; 2017MNRAS.469..930C, ). The solution to this is $dN_{\rm BBH}/dT\propto T^{5/3}$, which corresponds to,

[TABLE]

A population of BBHs entirely driven by GW radiation is therefore expected to have its period distribution increasing with increasing $\log T$. This behavior is opposite to what is found in our dynamical driven SC case.

To conclude, our derived $dN_{\rm BBH}/d\log T$ distribution of BBHs in SCs is clearly different from what is expected of BBHs not driven by dynamics. Therefore, if a distribution is measured that favors short period BBHs with a peak $T_{\rm peak}\sim$ days, then this would suggest that BBHs in SCs are present in the nearby universe and contribute to the observed GW merger rate. However, since the limit $T_{\rm min,ej}$ depends strongly on $v_{\rm dis}$, it is not entirely clear what the observed average value and dispersion is of the lower limit $T_{\rm min}$ for a representative GC ensemble. More detailed GC simulations must be used to provide further insight into this.

VI.2 Probing the BBH Orbital Time

Several observational probes exist that can be used to infer the orbital time distribution of BBHs in SCs. In the sections below we start by discussing prospects and challenges for using EM signals following a stellar TD event, after which we sketch out a few other possibilities.

VI.2.1 Stellar Tidal Disruptions

A star disrupted by a BBH is likely to give rise to a EM signal that initially follows the standard $t^{-5/3}$ decline (Rees1988, ; 2013ApJ…767…25G, ), after which the orbital motion of the two BHs gives rise to interruptions after a time $t_{\rm itr}$ that is comparable to the BBH orbital time $T$. This interestingly suggests that one is able to infer $T$ from observations of the interruption time $t_{\rm itr}$ (e.g., Liu:2009fl, ; 2014ApJ…786..103L, ; 2018arXiv181201118L, ). This approach was applied in 2014ApJ…786..103L , who studied the TDE J1201+30 which started out as a canonical TD after which clear interruptions appeared $\approx 27$ days after discovery. Using a simple $N$-body code for resolving the stream dynamics (Liu:2009fl, ), the authors estimated the orbital time and eccentricity of the disrupting super-massive BBH (SMBBH) to be $\approx 150$ days and $\approx 0.3$, respectively. The scenario of star-BBH TDs has been studied for both TDs by SMBBHs ($m_{\rm BH}\sim 10^{7}M_{\odot}$) (e.g. Liu:2009fl, ; 2016MNRAS.458.1712R, ; 2017MNRAS.465.3840C, ; 2018MNRAS.477.4009D, ) and recently by stellar mass BBHs ($m_{\rm BH}\sim 15M_{\odot}$) (2018arXiv181201118L, ), and generally it is found that the rate of accretion indeed is linked to the BBH’s orbital parameters. However, it is far from clear to what degree the EM light-curve actually traces the rate of mass fallback $\dot{M}$, which often is taken as a proxy for the EM emission from the TDE (2018arXiv180110180S, ). If the emission is beamed through a relativistic jet, one also expects light-curve variations from both precession and orbital motion caused by the companion BH and possible also the spin of the disrupting BH (e.g. 2011Natur.476..421B, ; 2013PhRvD..87h4053S, ). As further discussed in 2018arXiv181201118L and briefly below, we are likely only to observe jetted TDEs, and the jet dynamics for a disrupting stellar mass BBHs is therefore important to understand. This is a non-trivial problem and we therefore keep that for a future study.

We now turn to the astrophysical rates of star-BBH TDs expected from GCs. This rate depends both on the number density of stars, the absolute number of BBHs, but not on the the distribution of SMA a of the BBHs (see Section II.4). A simple ‘$n\sigma v$’ estimate gives,

[TABLE]

where this rate is per galaxy ($5$ BBHs per GC, and $200$ GCs per galaxy) derived for solar type stars ($1M_{\odot},1R_{\odot}$) interacting with BBHs of equal mass. We have further assumed that $P_{\rm TD}=2R_{\rm TD}/a$. This rate is on the lower side compared to both the derived rate of star-SMBH TDs (e.g., 2014ApJ…794….9M, ; 2017ApJ…841..132L, ; 2018MNRAS.480.5060S, ) and long gamma ray burst (GRB) (e.g. 2010MNRAS.406.1944W, ; 2018arXiv181201118L, ), but not unreasonably low for being observationally interesting. However, our derived rate is the rate of disruptions and not necessarily the rate of observable events. The observable rate depends further on the TD luminosity ($L$) and the energy extraction mechanism, which still are unsolved problems (e.g. 2014ApJ…783…23G, ; 2018arXiv180110180S, ). If $L$ is set by the mass fall back rate of the tidal stream, $\dot{M}$, then $L\propto 1/\sqrt{m_{\rm BH}}$, which suggests that lighter BHs should result in a higher $L$ than heavier BHs. However, if $L$ is Eddington limited, then $L\propto m_{\rm BH}$, which immediately leads to rates that are far below observable limits (2018arXiv181201118L, ). One proposed model for how radiation is processed and escapes involves an Eddington limited accretion disk with a relativistic launched jet (e.g. 2009ApJ…697L..77R, ; 2011MNRAS.416.2102G, ; 2012ApJ…760..103D, ). In this picture the energy flux carried by the jet can easily exceed the Eddington limit and will likely scale as $\propto\dot{M}c^{2}$, which would make such events visible at cosmological distances. However, in this case the observer has to be located near the cone of the jet opening, which significantly reduces the probability for detection. Taking this effect into account, the observed rate is likely $1/50-1/100$ times smaller (e.g. 2010MNRAS.406.1944W, ; 2017ApJ…849L..29K, ) than the disruption rate given by Eq. (24). However, the jet might cover a larger area on the sky if precession and movement from the BBH orbital motion is significant, which would lead to a higher probability for detection (e.g. 2013PhRvD..87h4053S, ). All in all, compared to long GRBs, which have been observed extensively by Swift, our derived rate of star-BBH TDs from GCs is about two-orders of magnitude lower. The rate could be higher if a subpopulation of the SCs have central stellar densities that are $>{10^{4}\text{pc}^{-3}}$, which could be the case for NSCs and GCs without a dominating BH core (e.g., Kremer2018d, ). The average number of GCs per galaxy could also be $>200$. Furthermore, active star formation in NSCs leads to a population of stars with large radii which also would lead to enhanced TD rates.

Finally, if both star-BH and star-BBH TDs are observed, then their number ratio can be used to constrain the fraction of BBHs in GCs, which play a key role in determining the BBH merger rate observable by LIGO/Virgo. To see how, we start by writing the rate of star-BBH TDs from a single GC as $\Gamma_{\rm TD}^{\rm BBH}\propto N_{\rm BBH}\times 4\eta_{\rm s}m_{\rm BH}R_{\rm TD}/v_{\rm dis}$, which follows from the arguments given in Section II.4, and where we have assumed that $P_{\rm TD}=2R_{\rm TD}/a$. Similarly, the total rate of star-BH TDs can be written as $\Gamma_{\rm TD}^{\rm BH}\propto N_{\rm BH}\times\eta_{\rm s}m_{\rm BH}R_{\rm TD}/v_{\rm dis}$ (2016ApJ…823..113P, ). The ratio between the two rates is therefore given by $\Gamma_{\rm TD}^{\rm BBH}/\Gamma_{\rm TD}^{\rm BH}\approx 4N_{\rm BBH}/N_{\rm BH}=4f_{\rm BBH}$, where $f_{\rm BBH}$ is the number of BBHs compared to single BHs. There are many uncertainties here, such as single BHs and BBHs do not necessarily cluster in the same way throughout the core region, and not all star-BBH TDs will be observationally different from star-BH TDs, but our simple relation does suggests that important astrophysical properties may be extracted from BH and BBH stellar TDs.

VI.2.2 Alternative Probes

Other ways of observing the BBH orbital distribution includes micro-lensing (e.g., 2002ApJ…579..639B, ; 2015ApJ…810L..20M, ; 2016MNRAS.458.3012W, ), periodic self-lensing if the BBHs accrete matter (e.g., 2018MNRAS.474.2975D, ; DOrazDiStef_InPrep, ), and GW emission (e.g. 2018arXiv180708864S, ; 2018arXiv181111812K, ). The possibility for detecting accreting BBHs through self-lensing is especially interesting, as this allows one to observe the remnant system of a star-BBH TD over a timespan comparable to the life time of the accretion disks, which is orders of magnitude longer than the TD flares and jets we have discussed so far. The probability for observing such systems is proportional to their life time, and accreting BBHs will therefore in principle have a much higher probability for being detected; however, for self-lensing to be efficient one has to observe near the orbital plane of the BBH, which puts strong constrains on the geometry in a similar way as for the jetted TDs. It is not yet clear if such requirements will reduce the number of sources below observable limits. Other approaches for indirectly probing a possible central population of BHs is to look for correlations with the GC stellar properties (e.g., 2018MNRAS.478.1844A, ; Kremer2018b, ; Kremer2018d, ; Weatherford2018, ), for which machine learning algorithms might be particular useful (e.g., 2018arXiv181106473A, ).

VII Conclusions

The orbital time distribution of BBHs in SCs is tightly linked to the GR dynamics that also drive the BBHs to merger. Therefore, if the BBH orbital distribution can be probed, then one will in principle be able to indirectly constrain the current BBH merger rate from SCs observable by LIGO/Virgo, as well as the dynamical environment of the BBHs. In this paper we propose that the BBH orbital time distribution might be indirectly probed using interrupted stellar TDs (e.g., Liu:2009fl, ). The main idea behind this ‘multi-messenger’ observable is that light-curves from star-BBH TDs are expected to show interruptions after a time-scale comparable to the BBH orbital time due to deflections of the tidal stream by the BHs.

Using a novel analytical framework for describing the dynamical evolution of BBHs in SCs, including the possibility for GW merger during hardening, we have shown that the corresponding BBH orbital time distribution have a near universal shape that peaks at a time-scale of days (Section IV). The distribution differs significantly from other astrophysical scenarios not driven by dynamics, including the Opik’s Law distribution (flat in $\log T$), purely GW driven BBHs (increasing with $\log T$), and SMBBH TDEs (distribute at time-scales of $10^{2}$ days). This suggests that if interrupted TDs are observed with breaks in their light-curves appearing after a time of hours to days (the time of interruption is usually found to be $\ll T$ (Liu:2009fl, )), then this would indicate a nearby dynamical driven population of BBHs embedded in a halo of stars, such as a GC. Actively searching for interrupted TD light-curves in current and upcoming EM surveys might therefore help us indirectly constrain the dynamical formation of BBH mergers in nearby SCs. In addition, the fraction of BBHs to single BHs can also be determined by probing the relative rate of star-BBH to star-BH TDs (Section VI).

We acknowledge that there are several uncertainties related to our proposed method, where the greatest is likely related to the TD luminosity and thereby the observable rate. As discussed in Section VI, super Eddington accretion is clearly needed to create rates that are within observable limits (see also: 2016ApJ…823..113P, ; 2018arXiv181201118L, ); however, exactly how a TD is powered even for single BH TDs is still an open question (e.g., 2015ApJ…806..164P, ). Despite this, the encouraging results from both the SMBBH limit (e.g., 2014ApJ…786..103L, ) and the stellar mass BBH limit (2018arXiv181201118L, ) make our proposed method an interesting possibility, and we do hope that our work will inspire groups to perform complementary studies of both the rate of star-BBH interactions in dense SCs as well of the accretion process itself.

Finally, this paper also presents the first GR analytical solution to the orbital time distribution of BBHs in dense SCs (Section II.5). We have compared our analytical solution to both a semi-analytical MC approach (Section II.5) and to data extracted from the Hénon-type MC code CMC (Section V). We find excellent agreement. As shown in Section IV, our inclusion of GR effects (representing the possibility for GW merger inside the cluster before ejection) gives rise to significant changes compared to the purely Newtonian limit (Fig. 3). This clearly illustrates the need for GR numerical methods for describing the distribution of BBHs even in classical GC systems (e.g., 2018arXiv181104926R, ). We hope our derived distribution will turn out to be useful for other studies related to BBHs in dense SCs and their detectability.

Acknowledgments. — The authors thank Carl L. Rodriguez, S. Schrøder, T. Fragos, B. Mockler, S. I. Mandel, W. Farr, C. Miller, D. J. D’Orazio, K. Hotokezaka, M. Gaspari and A. Askar for stimulating discussions. JS acknowledges support from the Lyman Spitzer Fellowship. TV acknowledges the support of the Friends of the Institute for Advanced Study. MLJR acknowledges that all praise and thanks belongs to Allah (any benefit is due to God and any shortcomings are my own). ERR acknowledge support from the DNRF (Niels Bohr Professor) and NSF grant AST-1615881. KK acknowledges support by the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE-1324585. The authors further thanks the Niels Bohr Institute for its hospitality while part of this work was completed, and the Kavli Foundation and the DNRF for supporting the 2017 Kavli Summer Program.