LISA Detection of Binary Black Holes in the Milky Way Galaxy

TL;DR

This paper estimates the number of tightly bound binary black holes in the Milky Way and discusses their potential detection by LISA, including possible X-ray signatures, based on merger rates from LIGO.

Contribution

It provides the first estimate of the abundance of close binary black holes in the Milky Way and assesses their detectability with LISA.

Findings

LISA could detect binary black holes in the Milky Way.

Tightly bound binaries likely originate from larger separations.

Possible X-ray signatures of such binaries are identified.

Abstract

Using the black hole merger rate inferred from LIGO, we calculate the abundance of tightly bound binary black holes in the Milky Way galaxy. Binaries with a small semimajor axis ($\lesssim 10 R_\odot$) originate at larger separations through conventional formation mechanisms and evolve as a result of gravitational wave emission. We find that LISA could detect them in the Milky Way. We also identify possible X-ray signatures of such binaries.

Equations

$$dN(a,t)=\rho (a,t)da.$$

$$\frac{\partial \rho (a,t)}{\partial t}-\frac{\partial }{\partial a}\left[K(M_1,M_2)\frac{\rho (a,t)}{a^3}\right]=S(a,t),$$

$$K(M_1,M_2)\equiv \frac{64}{5}\frac{G^3}{c^5}(M_1{M}_2)(M_1+M_2).$$

$$\rho (a,t)=a^{3}F\left[\frac{a^4}{4K}+t\right],$$

$$\frac{\partial }{\partial a}\left[K\frac{\rho (a)}{a^3}\right]=0.$$

$$\rho (a)=C\frac{a^3}{K},$$

$$\rho (a_m)\frac{K}{a_m^3}=R.$$

$$\rho (a)=\rho (a_m){\left[\frac{a}{a_m}\right]}^3=\frac{R}{K}{a}^3,$$

$$N(\le a)\approx 3\times 10^{-2}{R}_{100}{\left[\frac{a}{{\mathrm{R}}_{\odot }}\right]}^4,$$

$$\rho (a,t)=-\frac{1}{K}{a}^3\left[A\Theta (a-a_0){\int }_{a_0}^a\tilde{S}da\right]+a^{3}F\left[\frac{a^4}{4K}+t\right],$$

$$\rho (a,t)=-\frac{1}{K}{a}^3\left[A\frac{\Theta (a-a_0)}{(n+1)}(a^{n+1}-a_0^{n+1})\right]+a^{3}F\left[\frac{a^4}{4K}+t\right].$$

$$\rho (a,t)=-\frac{1}{K}{a}^{3}A\Theta (a-a_0)+a^{3}F\left[\frac{a^4}{4K}+t\right].$$

$$S(a)\propto a^n\exp \left[-\frac{a}{a_c}\right],$$

$$\frac{\partial \rho (a,t)}{\partial t}-\frac{\partial }{\partial a}\left[K(M_1,M_2)\frac{\rho (a,t)}{a^3}\right]=K_2{a}^n\exp \left[-\frac{a}{a_c}\right],$$

$$\rho (a,t)=a^{3}F\left[\frac{a^4}{4K}+t\right]+\frac{a^3{a}_c^{n+1}{K}_2}{K}\Gamma \left[1+n,\frac{a}{a_c}\right],$$

$$\rho (a)=a^{3}C+\frac{a^3{a}_c^{n+1}{K}_2}{K}\Gamma \left[1+n,\frac{a}{a_c}\right].$$

$$\rho (a_m)\frac{K}{a_m^3}+\sigma =R,$$

$$\sigma$$

$$=K_2{a}_c^{n+1}\left[\Gamma (1+n)-\Gamma \left(1+n,\frac{a_m}{a_c}\right)\right].$$

$$C=\frac{R-K_2{a}_c^{n+1}\Gamma \left(1+n,a_m/a_c\right)}{K}.$$

$$\rho (a)$$

$$\approx a^3\frac{R}{K},$$

$$\rho (a)={\int }_{K_{min}}^{\infty }{f}_K{\rho }_K(a,K)dK,$$

$${\int }_{K_{min}}^{\infty }{f}_{K}dK=1,$$

$$K_3=\frac{|m+1|}{K_{min}^{m+1}}=|m+1|K_{min}^{|m+1|}.$$

$${\int }_{K_{min}}^{\infty }{f}_K{\rho }_K(a_m,K)\frac{K}{a_m^3}dK=R.$$

$$R$$

$$=\frac{C{K}_3}{(m+1)}{K^{m+1}|}_{K_{min}}^{\infty }.$$

$$R=-\frac{C{K}_3}{(m+1)}{K}_{min}^{m+1},$$

$${\rho }_K(a,K)=-\frac{R(m+1)}{K_3{K}_{min}^{m+1}}\frac{a^3}{K}.$$

$$\rho (a)$$

$$=\frac{(m+1)}{m}\frac{R{a}^3}{K_{min}}.$$

Full text

1 Introduction

The Laser Interferometer Gravitation-Wave Observatory (LIGO) discovered gravitational waves from binary black holes (Abbott et al., 2016a, c), composed of black holes with masses $\gtrsim 10M_{\odot}$.

Two stars in an isolated binary can evolve to produce the progenitors of the LIGO sources. To possess the observed parameters of the LIGO sources, these binaries must have progenitors with high masses ($M\sim 40-100M_{\odot}$) and low metallicities (Belczynski et al., 2016). The binary evolution could be initially affected by mass transfer through a common envelope phase.

However, in the chemically homogeneous evolution model (de Mink and Mandel, 2016; Mandel and de Mink, 2016), two massive stars in a near contact binary spin rapidly due to tidal spin-orbit coupling. The rapid rotation of a star mixes its interior, allowing transport of hydrogen from the envelope to the core and metals from the core to the envelope (Maeder, 1987). In contrast to the standard binary evolution model, the stars do not follow a common envelope phase, due to their contraction within their Roche lobes.

Regardless of the formation mechanism, none of the progenitor systems could directly produce two black holes at arbitrarily small separations. Here we focus on tight binary black holes with a semimajor axis smaller than could feasibly be created directly by conventional stellar evolution mechanisms. These binary black holes were born at a larger semimajor axis through standard evolution, and then migrated to smaller separations via gravitational wave emission. Most of these binary black holes reside in an intermediate regime, where their coalescence time is shorter than the Hubble time but longer than the LIGO operation lifetime before they become detectable in the LIGO frequency band.

We focus our analysis on binary black holes within the Milky Way galaxy. While most of our equations could be applied to binaries at arbitrary distances, the signatures of the systems under consideration are not observable outside of the Milky Way. We will assume circular orbits, as the detected LIGO binaries are constrained to possess low eccentricities (Abbott et al., 2016a). This means that we neglect binary black hole production via many body encounters and strong interactions in globular clusters (Sigurdsson and Hernquist, 1993; Rodriguez et al., 2015). Aside from these assumptions, we will remain agnostic as to the specific mechanism producing the binaries.

Recently, Seto (2016) used an estimate of the population of Galactic binary black holes to predict that LISA will have the sensitivities required to detect binaries like GW150914. In this article, we used a more sophisticated population analysis to show that the result can be generalized to more complicated population models. Furthermore, we extend the calculation to black hole populations with varying masses and take into account the black hole mass function.

In section 2 we perform a population analysis of tight binary black holes in the Milky Way. In section 3 we show that the gravitational wave signatures of these tight binaries are observable by LISA. In section 4 we examine possible X-ray signatures of these binaries. Finally, section 5 summarizes our conclusions.

2 Population analysis

The number of binary black holes at a given time $t$ with semimajor axis between $a$ and $a+da$ can be written as,

[TABLE]

Assuming that their dynamical evolution is dominated by the emission of gravitational radiation, $\rho(a,t)$ obeys a simple advection equation,

[TABLE]

where $M_{1}$ and $M_{2}$ are the masses of the two black holes, $S(a,t)$ is a source term that parameterizes the production of binaries at a semimajor axis $a$, and $K$ is given by,

[TABLE]

In the simple case of $S=0$, the solution of equation (2) is given by

[TABLE]

where $F$ is some arbitrary function. The simplest solution can be obtained in a steady state, where equation (2) reduces to

[TABLE]

Integrating this equation gives

[TABLE]

where $C$ is an arbitrary constant. This is a special case of the solution in equation (3), with $F=C/K$.

The inferred merger rate from LIGO for two $30M_{\odot}$ black holes is between $2$ and $600$ Gpc*-3* yr*-1* (Abbott et al., 2016b) in comoving units. We can estimate the Galactic merger rate by adopting the number of Milky Way-like galaxies to be $10^{-2}$ per comoving Mpc3 (Montero-Dorta and Prada, 2009). Adopting $100$ Gpc*-3* yr*-1* as a fiducial LIGO inferred merger rate gives the LIGO Galactic merger rate to be $\sim 10^{-5}R_{100}$ mergers per galaxy per year.

A merger is detected by LIGO when the semimajor axis of the binary is small enough that it enters the LIGO frequency band. Denoting this critical semimajor axis as $a_{m}$, the merger rate $R$ is equal to the flux in $a$-space at $a=a_{m}$,

[TABLE]

Therefore,

[TABLE]

where we have substituted $\rho(a_{m})$ from equation (6). Using the inferred Galactic merger rate and specializing to binary black holes of mass $M_{1}=M_{2}=30M_{\odot}$, we can integrate over $a$ to obtain the number of Galactic binary black holes with semimajor axis $\leq a$,

[TABLE]

where $R_{\odot}$ is the solar radius and $R_{100}$ is the rate in units of $100\;\rm Gpc^{-3}\;yr^{-1}$.

2.1 Source functions with a minimum injection scale

Next, we generalize the population analysis to cases with a source function. Most formation mechanisms cannot produce binaries that are very tight. We model this situation with a source function that is proportional to a step function, $S=\tilde{S}(a)\Theta(a-a_{0})$, where $a_{0}$ is the minimum binary separation for the formation mechanism, $\Theta$ is the Heavyside step function, and $\tilde{S}(a)$ is an arbitrary function.

Given this source function, the general solution is,

[TABLE]

as long as the function $\tilde{S}(a)$ is not singular at $a=a_{0}$. For example, a power law source with normalization $A$ and index $n$, $S=Aa^{n}\Theta(a-a_{0})$, admits the general solution,

[TABLE]

Note also that in the case of a delta function source injected at $a=a_{0}$, the solution is given by,

[TABLE]

The most important feature of equation (9) is that the solution below the injection point, $a<a_{0}$ is unchanged from the sourceless case. This implies that as long as we restrict our analysis to $a\leq a_{0}$, we can simply use the sourceless solution. Our results will therefore be robust due to its insensitivity to the particular binary black hole production mechanism.

2.2 Power-law source functions with a maximum injection scale

For the sake of generality, we also consider a source function without a minimum scale, namely $S(a)\propto a^{n}$ with positive $n$ that extends all the way to $a=0$. In principle, $n$ is related to the power law index of the binary separation of massive stars (Sana et al., 2012). However, as not all massive binaries evolve into binary black holes, the mapping between the two indices is unknown. Since binaries are not produced up to arbitrarily high semimajor axis, we truncate our source function at high values of $a$ by an exponential factor,

[TABLE]

where $a_{c}$, is some large semimajor axis above which binary black holes are rarely produced. This source function corresponds to a mechanism that produces binaries over a broad range of $a$, where instead of a minimum injection scale, $a_{0}$, we now have a maximum injection scale, $a_{c}$. At large semimajor axes, this distribution corresponds to the end states of binary star evolution, whereas at small $a$ it corresponds to more exotic processes such as direct collapse (Loeb, 2016).

The binary population with this source term obeys

[TABLE]

where $K_{2}$ and $n$ are constants that in principe can be constrained by observations. The general solution to this equation is given by,

[TABLE]

where $F$ is an arbitrary function and $\Gamma$ is the incomplete Gamma function. As before, we can look for a steady state solution by setting the function $F$ to be a constant $C$, giving

[TABLE]

In this case, the merger rate $R$ equals the flux in $a$-space at $a=a_{m}$ plus a term corresponding to the source,

[TABLE]

where $\sigma$ is the rate of binary black holes created with $a\leq a_{m}$, given by,

[TABLE]

Since $a_{m}\ll a_{c}$, we get,

[TABLE]

The number density of binary black holes is then given by,

[TABLE]

where in the second equality we used the fact that both $a/a_{c}$ and $a_{m}/a_{c}$ are much smaller than unity. This result shows that in the case of a power law source function, the scaling of the sourceless solution $\rho\propto a^{3}$ remains valid even if the source function does not have a minimum scale as long as there is a maximum injection scale. The case of a power law with neither a maximum or minimum injection scale is treated in the Appendix.

2.3 Population analysis for varying black hole masses

In the previous sections, the abundance of binaries was derived for a given value of $K(M_{1},M_{2})$. In reality, the black hole population spans a range of masses with a probability given by the black hole mass function. Thus, the binaries possess varying values of $K$. In this section we incorporate this diversity of $K$ values. Note that there is a single advection equation for every value of $K$, i.e. there are many copies of equation (2), each for a different value of $K$. To make explicit its dependence on $K$, we label the density $\rho$ in equation (2) as $\rho_{K}(a,K)$. The total number of binary black holes per semimajor axis is given in terms of $\rho_{K}(a,K)$ by,

[TABLE]

where $f_{K}$ is the probability of finding a binary black hole with a particular value of $K$. As $K=K(M_{1},M_{2})$ is a function of the two black hole masses, its probability distribution is determined by the mass functions of the first and second black holes, $f_{M_{1}}$ and $f_{M_{2}}$, respectively. For simplicity, we adopt $f_{M_{1}}=f_{M_{2}}=\Phi_{M}$, and the distribution $f_{K}$ can be obtained from $\Phi_{M}$ by a series of convolutions. Adopting a phenomenological power-law relation, $f_{K}=K_{3}K^{m}$ with an index $m<-1$ and a normalization,

[TABLE]

we can derive the normalization constant $K_{3}$ in terms of $K_{min}$. Theoretically, the minimum $K$ value is obtained when both black holes are at the limit imposed by the Tolman-Oppenheimer-Volkoff equation of a Chandrasekhar-Landau mass ($\sim 3M_{\odot}$) each. However, there is evidence that there exists a mass gap under $\sim 5M_{\odot}$ (Ozel et al, 2010; Farr et al, 2011; Kreidberg et al., 2012). We thereby chose the minimum $K$ to be that when both black holes are $\sim 5M_{\odot}$ each.

For our phenomenological model with $m<-1$, the integral converges and can be solved to give,

[TABLE]

The merger rate is given by the sum over all masses of the fluxes in $a$-space,

[TABLE]

For the sourceless steady-state solution, $\rho(a)=Ca^{3}/K$, we can find $C$ in the phenomenological $f_{K}=K_{3}K^{m}$ model by noting that,

[TABLE]

For $(m+1)<0$, this integral converges to

[TABLE]

which implies that for $(m+1)<0$,

[TABLE]

The number of binary black holes per unit semimajor axis is therefore given by,

[TABLE]

The number of Galactic binary black holes with semimajor axis $\leq a$ is given by

[TABLE]

Aside from numerical factors of order unity the only change from the single $K$ case is that $K_{min}$ appears in the denominator in place of $K$. Since $K_{min}$ is $\sim 200$ times smaller than the $K$ for two $30$ solar mass black holes, this number is of order $200$ larger than in equation (8).

2.4 Population analysis for varying black hole masses: Chabrier/Kroupa IMF

Next, we proceed beyond the phenomenological toy model for $f_{K}$ assuming that the black hole mass function follows the power-law dependence of massive stars, $\Phi_{M}=kM^{-2.3}$ (Chabrier, 2003), where $k$ is a constant. We define the quantity,

[TABLE]

so that the respective distribution $f_{\tilde{K}}$ follows,

[TABLE]

If the masses $M_{1}$ and $M_{2}$ are independently distributed, $f_{K}$ can be derived from $\Phi_{M}$ as follows. We first switch from the random variables $M_{1}$ and $M_{2}$ to $\tilde{K}$ and $W$, where $W\equiv M_{1}$ and,

[TABLE]

where the positive root was chosen since $W$, $M_{2}$, and $\tilde{K}$ are positive definite. The distribution function $f_{\tilde{K}W}$ is then given by,

[TABLE]

where the determinant of the Jacobian of the transformation,

[TABLE]

The marginal distribution $f_{\tilde{K}}$ is therefore,

[TABLE]

which for the assumed mass function is given by,

[TABLE]

Here, $W_{\rm min/max}$ corresponds to the minimum and maximum black hole masses; in particular, $W_{\rm min}$ is again chosen to be $\sim 5M_{\odot}$ and $W_{\rm max}\sim 100M_{\odot}$, respectively.

Although the integral in equation (35) could only be solved numerically, the distribution is well represented by a power law form between $W_{\rm min}$ and $W_{\rm max}$ with a negative index $-\alpha\sim-2.1$ (see Figure 1). In order to simplify the analysis, we will therefore adopt,

[TABLE]

where $A=159$. This gives

[TABLE]

We proceed analogously to the previous section, where the main difference is that we now have an upper cutoff on black hole masses at $W_{\rm max}$.

For the steady state solution, $\rho(a)=Ca^{3}/K$, the rate of binary black holes entering the LIGO band is given by,

[TABLE]

In equation (2.4), the integration limits refer to the maximum and minimum black hole masses that LIGO is sensitive to. The assumption that LIGO is sensitive to all black hole masses available in the black hole mass function translates to substituting $K_{max}$ and $K_{min}$ as these limits. Note that as the IMF is dominated by low mass black holes, our result is only weakly dependent on the exact value of $K_{max}$. Using $R$ to eliminate the constant $C$ yields,

[TABLE]

The abundance of binaries is therefore given by,

[TABLE]

The number of binary black holes with semimajor axis $\leq a$ is then,

[TABLE]

While this result is derived by assuming that the black hole mass function follows the power-law mass function of massive stars, simulations indicate that not all massive stars form black holes, and that stars above $\sim 50M_{\odot}$ blow off too much of their mass to produce black holes (Sukhbold, 2016). These complications could lower the abundance of black holes relative to that predicted by equation (42). This means that our prediction should be interpreted as an upper bound to the number of binary black holes in the Galaxy.

2.5 Comparison with population synthesis models

In order to compare our result to population synthesis models, we first transform our variable from $a$ to the frequency $f$. The number of binary black holes with semimajor axis $\leq a_{\rm max}$ is equal to the number of binary black holes with frequency $\geq f_{\rm min}$, where $f_{\rm min}$ is the orbital period of the smallest black holes in the population with separation $a_{\rm max}$.

For the systems under consideration, $a_{\rm max}=10R_{\odot}$ and the smallest black hole mass is $5M_{\odot}$, resulting in a minimum frequency of $f_{\rm min}=2\times 10^{-5}$Hz. Equation (42) therefore predicts $\approx 3\times 10^{4}R_{100}$ binary black holes in the Galaxy with frequency greater than $2\times 10^{-5}$Hz.

The population synthesis model A of Belczynski et al. (2010) predicts $\sim$thousands of binary black holes with frequencies greater than $2\times 10^{-5}$Hz. Noting that the LIGO rate $R_{100}$ ranges from $0.02$ to $6$, this is consistent with our result. An earlier calculation by Nelemans et al. (2001) predicts a number that is an order of magnitude larger than model A of Belczynski et al. (2010), which is still consistent with our result.

3 Gravitational wave signal from Milky Way binaries

Most of the $3\times 10^{4}R_{100}$ Galactic binary black holes with orbital separation $a\leq 10R_{\odot}$, will not enter the LIGO bandpass in a short enough time for them to be observed by LIGO. For example, the timescale for two $30M_{\odot}$ black holes to coalesce from $a\sim$ a few $R_{\odot}$ is thousands of years. These binaries, however, will be observable by LISA111http://www.elisascience.org which is sensitive to lower frequencies than LIGO.

Focusing on the case of two $\sim 30M_{\odot}$ black holes, we find from equation (8) that the tightest binary black hole in our Galaxy has $a\sim 2.5R_{\odot}$. For such a binary consisting of two $30M_{\odot}$ black holes, the gravitational wave frequency is $f\sim 3\times 10^{-4}$ Hz, which is within the LISA bandpass (Farmer and Phinney, 2003). The angular-averaged gravitational wave strain for the $n=2$ mode is given by (Peters and Mathews, 1963; Seto, 2016),

[TABLE]

where $M_{c}\equiv(M_{1}M_{2})^{3/5}(M_{1}+M_{2})^{-1/5}$ is the chirp mass. Integrating the signal over an observational period $\tau$, the signal to noise ratio becomes (Seto, 2016),

[TABLE]

where $h(f)$ is the LISA instrumental noise, with a value of $\sim 3\times 10^{-18}\;\rm{Hz^{-1/2}}$ at $0.5$ mHz (Amaro-Seoane et al., 2012). Scaling the noise with frequency as the power law $h(f)\propto f^{-2}$ (Seto, 2016), the signal to noise ratio becomes,

[TABLE]

For observations across the Milky Way with $d=20$ kpc, we find,

[TABLE]

Figure 2 shows the expected number of such Milky Way binaries as a function of their SNR.

3.1 Confusion with cosmological sources

A supermassive binary black hole at cosmological distances can possess similar values of strain amplitude and frequency to a Galactic binary black hole, thus masquerading as a Galactic source. However, this confusion can be eliminated by measuring the change in gravitational wave frequency as a function of time, $\dot{f}$.

The time derivative of the gravitational wave frequency is given by (Cutler and Flanagan, 1994)

[TABLE]

Supermassive binary black holes possess chirp masses that are much greater than that of Galactic binaries. As such, their frequency changes at a much faster pace than Galactic binaries.

4 Electromagnetic flag

4.1 Binary black hole accretion in a hierarchical triple system

A binary black hole could accrete gas if it resides in a hierarchical triple system, where the third object is a main sequence star. The wind of the third star would lead to accretion at the Bondi-Hoyle-Lyttleton rate (Hoyle and Lyttleton, 1939; Bondi and Hoyle, 1944),

[TABLE]

where $M_{tot}=M_{1}+M_{2}$, $\rho_{w}$ the mass density of the stellar wind, $v_{w}$ the wind speed, and $c_{w}\equiv(5kT_{w}/3m_{p})$ is the sound speed, with $T_{w}$ be the wind temperature. Scaling the wind parameters to solar values at a distance of $1$ AU and assuming an efficiency of $0.1$ for converting rest mass into radiation, the luminosity produced by the binary where both black holes are $\sim 30$ solar masses is,

[TABLE]

where $l$ is the separation of the star from the binary. The maximum luminosity is given by the Eddington limit,

[TABLE]

Due to the orbital motion of the black holes around the center of mass, the observed flux would be modulated by Doppler beaming. Assuming that the emitted flux, $F_{\nu 0}$ scales with frequency as $F_{\nu 0}\propto\nu^{\beta}$, the Doppler modulation is given by

[TABLE]

where $F_{\nu}$ is the observed flux, and $D$ is the Doppler factor. To first order, the flux modulation is given by (D’Orazio, Haiman, and Schiminovich, 2015),

[TABLE]

where $\phi$ is the orbital phase and $i$ the inclination.

At the Eddington luminosity, the amplitude of the flux modulation of two $30M_{\odot}$ binary, is given by,

[TABLE]

where $d$ is the distance to the object and we have assumed $\beta\sim 1$ (Sobolewska et al., 2011; D’Orazio, Haiman, and Schiminovich, 2015). Given the flux sensitivity of XRM-Newton of $2\times 10^{-15}\;\rm erg\;cm^{-2}\;s^{-1}$, for a semimajor axis of $a\sim 10R_{\odot}$, the flux modulation of these objects will be observable out to $d\sim 10$ Mpc.

Realistically, it is unlikely for such binaries to emit at the Eddington luminosity. For general Bondi-Hoyle-Lyttleton accretion, the flux modulation depends on the binary and stellar wind parameters,

[TABLE]

Substituting the solar wind parameters at $d=1$ AU and the XRM-Newton sensitivity for $\Delta F$ yields for the observer distance of $d\sim 300$ pc at which a pair of $30M_{\odot}$ black holes will be detectable. The X-ray surveyor222http://wwwastro.msfc.nasa.gov/xrs/, a next generation x-ray observatory with a flux sensitivity of $\sim 10^{-19}\;\rm erg\;cm^{-2}\;s^{-1}$, will be able to detect this flux modulation out to a distance of $\sim 30$ kpc, which allows their detection throughout the entirety of the Milky Way galaxy.

Owing to the fact that most massive stars are in multiple systems (Raghavan et al., 2010; Tokovinin, 2014), and that there is precedent for X-ray binaries in triple systems (Grindlay et al, 1988; Thorsett et al, 1999; Chou and Grindlay, 2001; Zdziarski et al., 2007; Prodan and Murray, 2015), most binary black holes will likely possess a third companion. Further evidence of this comes from observations of superorbital modulations in high-mass X-ray binary systems, which could be caused by a third companion (Farrell, Sood, and O’Neill, 2006; Corbet and Krimm, 2013). However, only a fraction of binary black holes in a hierarchical triplet would host companions in the relevant mass range for accretion to be efficient. Since only companions with masses $\sim 1M_{\odot}$ and above generates sufficient luminosity to be observed throughout the Milky Way, the percentage of observable systems is

[TABLE]

where $\Phi_{M}$ is the stellar IMF, $T_{H}\sim 1.4\times 10^{10}$ yr is the Hubble time, and $T_{c}$ is the main sequence lifetime of the companion star, given by the broken power law (Laughlin, Bodenheimer and Adams, 1997; Salaris and Cassisi, 2006; Loeb, Batista, and Sloan, 2016)

[TABLE]

where $(\tau,\Gamma)=(10^{10}\;\rm{years},\;2.5)$ for stars less massive than $3M_{\odot}$ and $(\tau,\Gamma)=(7.6\times 10^{9}\;\rm{years},\;3.5)$ for more massive stars. The ratio in equation (56) takes into account both the stellar IMF and the fact that more massive stars live a shorter amount of time. Substituting the Chabrier IMF for $\Phi_{M}$, we obtain,

[TABLE]

This implies that out of the $\sim 300$ black holes with $a\leq 10R_{\odot}$ predicted by equation (8), only a few systems would host the appropriate companion to be observable throughout the entirety of the Milky Way. This number is further diminished by the fact that only a fraction of all triple systems have the stars at a close enough distance, as the signal scales as $\rho_{w}\propto l^{-2}$. We therefore conclude that the most efficient method to detect these binaries is through their gravitational wave emission with LISA.

4.2 Tidal disruption flares from planets and asteroids

Another source of electromagnetic activity could be associated with the tidal disruption of planets and asteroids by the black holes. White dwarfs and neutron stars are known to host rocky debris around them (Vanderburg et al., 2015; Zuckerman et al., 2010; Farihi et al, 2009; Koester et al., 2014; Wolszczan and Frail, 1992; Podsiadlowski, 1993; Phinney and Hansen, 1993). Orbits around binary black holes can be chaotic and are subject to Kozai-Lidov oscillations, leading to an enhanced tidal disruption rate (Ivanov et al., 2005; Chen et al., 2009; Li et al., 2015).

For a planet or asteroid of mass $m_{p}$ and radius $r_{p}$ being tidally disrupted by a black hole of mass $M$, the length of the flare is defined as the time it takes for the emission to drop under the Eddington limit. This is given by (Ulmer, 1998),

[TABLE]

where $l_{p}$ the pericenter distance and $l_{t}$ the tidal radius. Assuming a Neptune-like planet with $m_{p}\sim 10^{29}$ g and $r_{p}\sim 10^{29}$ cm, the flare time becomes $t_{f}\sim 1.4$ years.

Since this timescale is much longer than the orbital timescale, lightcurves from these events will possess amplitude modulation due to the Doppler effect, as described in section 4.1. The Doppler modulation of such systems is bright enough to be detectable from throughout the Milky Way by existing telescopes such as the XMM-Newton and the Chandra X-ray Observatory. As the number of sources in the sky per unit time will depend on the duty cycle of such flares, a monitoring campaign of a large patch of the sky is required to find such flaring sources. Confusion with other sources would make the identification of such binaries difficult.

5 Conclusions

By calibrating the population of binary black holes based on the merger rate infered by LIGO, we have found that LISA could detect a handful of such binaries in the Milky Way galaxy (see Figure 2). A lack of detections will set constraints on the binary production mechanisms.

We also considered electromagnetic flags of these tight binary black holes in the Milky way, and found them to have weaker observational prospects. Gravitational wave signals could be leveraged to provide both the system’s masses and semimajor axis.

6 Acknowledgements

The authors thank Josh Grindlay for comments on the manuscript. This work was supported in part by the Black Hole Initiative at Harvard University, funded by the John Templeton Foundation.

7 Appendix

7.1 Power law source functions without a minimum or maximum scale

Consider a nonzero, time-independent power-law source term $S(a)\propto a^{n}$ that extends all the way from $a=0$ to $\infty$. Unlike the example considered in the main text, we consider a case without a smallest injection scale or an upper cutoff scale. The binary black hole population with this source term obeys

[TABLE]

where $K_{2}$ and $m$ are constants that in principe can be either derived or estimated from observations. The general solution to this equation is given by

[TABLE]

where $F$ is an arbitrary function. As before, we search for a steady state solution by setting the function $F$ to be the constant $C$, giving

[TABLE]

In this case, the merger rate $R$ is equal to the flux in $a$-space at $a=a_{m}$ plus a term corresponding to the source term,

[TABLE]

where $\sigma$ is the rate of binary black holes created per year with semimajor axis $a\leq a_{m}$, given by

[TABLE]

In this case,

[TABLE]

which allows us to deduce that $C=R/K$. Therefore, the number of Galactic binary black holes with semimajor axis $\leq a$ is given by

[TABLE]

where $N_{\rm hom}(\leq a)$ is the number of Galactic binary black holes in the source-less case. Note that when we set the source term to zero by using $K_{2}=0$, we will recover the sourceless solution. The effect of such a source term with $n>0$, i.e. where there is a higher rate of binary black hole production at large semimajor axis, is paradoxically a suppression of the number of binary black holes in the steady state solution due to the requirement that the rate $R$ be kept unchanged. Since a source function that does not have a maximum injection scale is unphysical, the solution we obtain is also unphysical. In this case, there is a scale, $a_{\rm crit}$, above which $N(\leq a)$ becomes negative. This calculation should be viewed only as a pedagogical toy model to illustrate an example of the ramifications of modifying one of our assumptions.

7.2 Population analysis for varying black hole masses: Power law source

We can repeat the analysis for the case of a power law source term. In this case we have,

[TABLE]

where, inspired by our previous solutions, we have explicitly written out the $1/K$ dependency of $C$ so that it is now a constant with respect to $K$. The merger rate is therefore given by,

[TABLE]

which converges when $(m+1)<0$ to

[TABLE]

Here $\sigma$ is again the rate of binary black holes created per year with semimajor axis $a\leq a_{m}$. However, in order to be consistent with our choice of the black hole mass function, we need to take into account the fact that different amounts of binary black holes are created for different black hole masses. Equivalently, binary black holes with different $K$’s are produced at different abundances. As a result, $\sigma$ has to include an extra integral over $K$,

[TABLE]

When $m+1<0$, this integral converges to,

[TABLE]

Following through, this gives the number of binary black holes per unit semimajor axis to be,

[TABLE]

when $m<0$. Note that when we set the source term to zero through $K_{2}=0$, we will recover the sourceless solution. The number of binary black holes per unit semimajor axis is therefore,

[TABLE]

where $\rho_{\rm{hom}}(a)$ is the solution in the sourceless case and the number of binary black holes with semimajor axis $\leq a$ is given by,

[TABLE]

where as before $N_{\rm{hom}}$ is the sourceless solution. This solution is again pathological due to the presence of a critical semimajor axis above which $N(\leq a)$ becomes negative.