Gravitational-wave sources from mergers of binary black-holes catalyzed by fly-bys interactions in the field

TL;DR

This paper introduces a new mechanism for black hole binary mergers driven by stellar fly-bys that excite eccentricities in wide binaries, leading to gravitational-wave emissions, with a rate comparable to observations.

Contribution

It presents an analytic model for wide binary black hole mergers catalyzed by stellar fly-bys, a novel formation channel supported by numerical verification.

Findings

Merger rate estimated at ~10 times the fraction of wide binaries per Gpc^3 per year.

Predicted observational signatures include spin-orbit misalignment and high mass-ratio preferences.

The model's merger rate aligns with current observational estimates.

Abstract

Several scenarios were suggested for the origins of gravitational-wave (GW) sources from mergers of stellar binary black holes (BBHs). Here we propose a novel origin through catalyzed formation of GW-sources from ultra-wide binaries in the field. Such binaries experience perturbations from random stellar fly-bys which excite their eccentricities. Once a wide-binary is driven to a sufficiently small peri-center approach, GW-emission becomes significant, and the binary inspirals and merges. We derive an analytic model and verify it with numerical calculation to compute the merger rate to be $\sim10{\rm \times f_{wide}\,{\rm Gpc}^{-3}yr^{-1}}$ ($f_{{\rm wide}}$ is the fraction of wide BH-binaries), which is comparable to the observationally inferred rate. The observational signatures from this channel include spin-orbit misalignment; preference for high mass-ratio BBH; preference for high velocity-dispersion host-galaxies; and a uniform delay-time distribution.

Tables

Table 1. Table 1: This table summarize the numerical results. In all runs we chose b max = 5 × 10 4 ​ AU subscript 𝑏 max 5 superscript 10 4 AU b_{{\rm max}}=5\times 10^{4}{\rm AU} ,o ensure we simulate the impulsive regime. The stellar number density is n ∗ = 0.1 ​ pc − 3 subscript 𝑛 0.1 superscript pc 3 n_{*}=0.1{\rm pc^{-3}} and the velocity dispersion is ⟨ v enc ⟩ = 50 ​ k ​ m ​ s − 1 , delimited-⟨⟩ subscript 𝑣 enc 50 k m superscript s 1 \left\langle v_{{\rm enc}}\right\rangle=50{\rm kms^{-1}}, the perturber mass is m p = 0.6 ​ M ⊙ . subscript 𝑚 𝑝 0.6 subscript 𝑀 direct-product m_{p}=0.6M_{\odot}. n 𝑛 n - represents the number of mergers due to GW out of N 𝑁 N simulations with the same initial setup, namely a 𝑎 a and e 𝑒 e . The ratio of n / N 𝑛 𝑁 n/N is the estimate of the probability of merger for a specific a 𝑎 a and e 𝑒 e .

$a\left[{10}^4\mathrm{AU}\right]$		$e=0.8$	$e=0.9$	$e=0.95$	$e=0.99$	$e=0.995$	$e=0.999$	$e=0.9995$	$e=0.9999$
$a=0.5$	$n$	-	-	-	$0$	$34$	$470$	$722$	$1515$
	$N$	-	-	-	$24927$	$11348$	$24587$	$24682$	$24740$
$a=1$	$n$	-	$0$	$5$	$135$	$472$	$1801$	$158$	$3759$
	$N$	-	$100000$	$250430$	$99434$	$100367$	$94240$	$7491$	$125087$
$a=2$	$n$	$1$	$3$	$96$	$686$	$126$	$207$	$791$	$183$
	$N$	$595746$	$249568$	$248101$	$235553$	$24495$	$12474$	$35122$	$4999$
$a=3$	$n$	-	$0$	$16$	$204$	$393$	$943$	$327$	$56$
	$N$	-	$99950$	$99588$	$96550$	$94873$	$94915$	$24643$	$2494$
$a=4$	$n$	-	-	$2$	$62$	$85$	$207$	$64$	$61$
	$N$	-	-	$50000$	$48898$	$24522$	$29813$	$7488$	$4979$

Equations

$$t_{\mathrm{merger}}\approx \frac{a^4}{{\beta }^{\prime }}\times {\left(1-e^2\right)}^{7/2}$$

$${\beta }^{\prime }=\frac{85}{3}\frac{G^3{m}_1{m}_2\left(m_1+m_2\right)}{c^5},$$

$$e_c={\left[1-{\left(\frac{{\beta }^{\prime }T}{a^4}\right)}^{2/7}\right]}^{1/2}.$$

$$F_q={\int }_{e_c}^{1}2ede=1-e_c^2={\left(\frac{{\beta }^{\prime }T}{a^4}\right)}^{2/7}.$$

$$\theta =\frac{⟨\Delta v⟩}{v_k}$$

$$v_k={\left(\frac{G{m}_b}{3a}\right)}^{1/2}.$$

$$⟨\Delta v⟩\simeq \frac{3Ga{m}_p}{v_{\mathrm{enc}}{b}^2}$$

$${\theta }^2=\frac{9G^2{a}^2{m}_p^2}{{\left(v_{\mathrm{enc}}{b}^2\right)}^2}\frac{3a}{G{m}_b}=\frac{27G{a}^3{m}_p^2}{m_b{\left(v_{\mathrm{enc}}{b}^2\right)}^2}$$

$$F_s=\frac{\pi {\theta }^2}{4\pi }=\frac{27}{4}{\left(\frac{m_p}{m_b}\right)}^2\left(\frac{G{m}_b}{a{v}_{\mathrm{enc}}^2}\right){\left(\frac{a}{b}\right)}^4.$$

$$\frac{F_s}{F_q}=\frac{27}{4}{\left(\frac{m_p}{m_b}\right)}^2\left(\frac{G{m}_b}{a{v}_{\mathrm{enc}}^2}\right){\left(\frac{a}{b}\right)}^4{\left(\frac{a^4}{{\beta }^{\prime }T}\right)}^{2/7}$$

$${\dot{L}}_q=\frac{F_q}{T}.$$

$$f=n_{\ast }\sigma v_{\mathrm{enc}}.$$

$$n_{\ast }\pi b^2{v}_{\mathrm{enc}}=\frac{1}{T}.$$

$$\frac{F_s}{F_q}=\frac{27}{4}{\left(\frac{m_p}{m_b}\right)}^2\left(\frac{G{m}_b}{a{v}_{\mathrm{enc}}^2}\right){\left(\frac{a}{b}\right)}^4{\left(\frac{a^4}{{\beta }^{\prime }T}\right)}^{2/7}=1.$$

$${\left(v_{\mathrm{enc}}{b}^2\right)}^2\le \frac{27}{4}\frac{G{m}_p^2{a}^{29/7}}{m_b{\left({\beta }^{\prime }T\right)}^{2/7}}$$

$$a_{\mathrm{crit}}={\left[\frac{4}{27}\frac{m_b{\beta }^{\prime 2/7}{T}^{-12/7}}{G{m}_p^2{n}_{\ast }^2{\pi }^2}\right]}^{7/29}.$$

$$L_{a<a_{\mathrm{crit}}}=1-{\left(1-F_q\right)}^{t\cdot f}$$

$$L_{a<a_{\mathrm{crit}}}=tf{F}_q.$$

$$f=n_{\ast }\pi {\left(\frac{27}{4}\frac{G{m}_p^2{a}^{29/7}}{m_b{\left({\beta }^{\prime }T\right)}^{2/7}}\right)}^{1/2}$$

$$L_{a<a_{\mathrm{crit}}}=t{n}_{\ast }{m}_{\ast }{a}^{13/14}\times {\left(\frac{27G{\left({\beta }^{\prime }T\right)}^{2/7}}{4m_b}\right)}^{1/2}$$

$$L_{a>a_{\mathrm{crit}}}=1-{\left(1-F_q\right)}^{t/T}.$$

$$L_{a>a_{\mathrm{crit}}}=t{F}_q\frac{1}{T}=t{\left(\frac{{\beta }^{\prime }T}{a^4}\right)}^{2/7}\frac{1}{T}=t{a}^{-8/7}\times T^{-5/7}{\beta }^{\prime 2/7}.$$

$$t_{1/2}=0.00233\frac{v_{\mathrm{enc}}}{G{m}_p{n}_{\ast }a}.$$

$$L_{a<a_{\mathrm{crit}}}=\tau n_{\ast }{m}_{\ast }{a}^{13/14}\times {\left(\frac{27G{\left({\beta }^{\prime }T\right)}^{2/7}}{4m_b}\right)}^{1/2}\left(1-e^{-t/\tau }\right)$$

$$L_{a>a_{\mathrm{crit}}}=\tau a^{-8/7}\times T^{-5/7}{\beta }^{\prime 2/7}\left(1-e^{-t/\tau }\right),$$

$$n_{\ast }\left(r\right)=n_0{e}^{\left(-\left(r-r_{\odot }\right)/R_l\right)},$$

$$f_{\mathrm{BBH}}\approx f_{\mathrm{primary}}\times f_{\mathrm{secondary}}\times f_{\mathrm{wide}}\approx$$

$$2\times 10^{-4}\left(\frac{f_{\mathrm{primary}}}{10^{-3}}\right)\left(\frac{f_{\mathrm{seoncdary}}}{0.4}\right)\left(\frac{f_{\mathrm{wide}}}{0.5}\right).$$

$$\Gamma =\int \int L\left(a,r\right)\times f_{\mathrm{BBH}}\times \left(\frac{1Myr}{10\mathrm{Gyr}}\right)\mathrm{dadr}\approx 0.42My{\mathrm{r}}^{-1}.$$

$$\mathcal{R}=10^3{\rho }_{\mathrm{gal}}\times \Gamma \approx$$

Full text

I Introduction

Extensive theoretical studies over the past few decades have proposed the existence of gravitational-wave (GW) sources arising from the mergers of two compact objects, and provided a wide range of predicted production rates of such sources (e.g. Belczynski et al., 2002, 2004, 2007, 2008, 2016; de Mink and Belczynski, 2015; Dominik et al., 2012; Antonini and Perets, 2012; Antonini et al., 2014; Antognini et al., 2014; Petrovich and Antonini, 2017, and more). Observationally, eleven confirmed GW mergers have been detected by aLIGO and VIRGO since their initial operation. These include $9$ mergers of binary black-holes (BBHs) and a single merger from a binary neutron-star (NS)The LIGO Scientific Collaboration et al. (2018). The currently inferred BBH-merger rate from these observations (in the local Universe) is $\mathscr{R_{{\rm BBH}}}=9.7-101{\rm Gpc^{-3}yr^{-1}}$; while the merger rates of binary neutron-star is $\mathscr{R_{{\rm BNS}}}=110-3840{\rm Gpc^{-3}yr^{-1}}$; and the upper limit of BH-NS merger is $\mathscr{R_{{\rm BHNS}}}<600{\rm Gpc^{-3}yr^{-1}}$.

Three main evolutionary channels were proposed in the context of GW mergers. The first deals with merger in dense environments such as galactic centers or globular clusters (e.g. Rodriguez et al., 2016, 2018; Leigh et al., 2018), where binary mergers are catalyzed by strong interactions with stars in these dense environment. In such environments, strong three-body interactions lead to harden compact binaries (drive them to shorter periods) and excite their eccentricities. Such models predict GW-production rates in the range of $\sim 2-20{\rm Gpc^{-3}yr^{-1}}.$ The second evolutionary channel deals with the isolated evolution of initially massive close binary stars (e.g. Belczynski et al., 2008, 2016; Dominik et al., 2012, 2015; Mandel and de Mink, 2016). Some of the massive close binaries strongly interact through one or two common envelope phases (e.g. Dominik et al., 2012) in which the interaction of a star with the envelope of an evolved companion leads to its inspiral in the envelope and the production of a short period binary. A fraction of the post-CE binaries are sufficiently close to merge via GW emission within Hubble time. The large uncertainties in the initial conditions of the binaries, the evolution in the common-envelope phase, the natal-kick experienced by NS/BHs at birth; and the mass-loss processes of massive stars give rise to a wide range of expected GW-sources production rates in the range $\sim 10^{-2}-10^{3}{\rm Gpc^{-3}yr^{-1}}.$ The third evolutionary channel deals with mergers induced by secular evolution of triple systems either in the field (Antonini et al., 2017) or in dense environments (e.g. Samsing and D’Orazio, 2018; Samsing, 2018; Petrovich and Antonini, 2017; Antonini et al., 2018; Antonini and Perets, 2012). In this channel the secular perturbations by a third companion (Lidov-Kozai evolution (Lidov, 1962; Kozai, 1962) ) can drive BBHs into high eccentricities such that they merge within a Hubble-time; the rates expected in this channel are $\sim 0.5-15{\rm Gpc^{-3}yr^{-1}.}$

Here we present a fourth channel of binary evolution, in which we focus of wide (SMA >$1000{\rm AU})$ BBHs in the field perturbed by random fly-by interactions of field stars in their host galaxy. Kaib and Raymond (2014) and Michaely and Perets (2016) showed that although evolution of stars and binaries in the low-density environment in the field is typically thought to be collisionless, wide binary systems can be significantly affected by fly-by interactions of field stars stars, and effectively experience a collisional evolution. In particular, Kaib and Raymond (2014) calculated the probability of a head-on collision between two main-sequence stars in the Milky-Way Galaxy due to interaction with random stellar perturbers. Michaely and Perets (2016) followed these directions and suggested a novel formation scenario for low-mass X-ray binaries from wide-binaries in the field.

Here we show that collisional evolution in the field could also be highly important for the formation of GW-sources. We analyze the evolution of wide-orbit BBHs and show that a fraction of these can be driven into high eccentricities and close pericenter distances due to interaction with stellar perturbers. In cases where the pericenter distance of a given binary is driven into a sufficiently small distance GW-emission becomes significant and the binary rapidly loses angular momentum and energy due to GW emission, and eventually inspirals and merges as a GW-sources detectable by aLIGO/VIRGO.

This paper is organized as follows. In section II we present the analytic model, the basic assumptions and the calculations. In section III we present the numerical verification to the analytic model. We discuss the results and summarize in section IV. The numerical procedure, the equations we integrate and the data analysis are described at length at appendix Appendix I and Appendix II.

II Analytic model

II.1 Formation Scenario

We consider a wide BBH with semi-major axis (SMA) $a>10^{3}{\rm AU}$. The binary resides in the field of the host galaxy and therefore be affected by short duration dynamical interactions with field stars. The dynamical encounters can typically be modeled through the impulse approximation, i.e. in the regime where the interaction timescale $t_{{\rm int}}\equiv b/v_{{\rm enc}}$(where $b$ is the closest approach to the binary and $v_{{\rm enc}}$ is the velocity of the perturbing mass) is much shorter than the BBH orbital period time, $P$. These perturbations can torque the system and exchange orbital energy thereby decreasing/increasing the binary semi-major axis $a$ and the binary eccentricity $e$. If these interactions drive the system to a sufficiently small pericenter passage, $q$, then the system can merge via GW emission within Hubble time.

There are four relevant timescales for this impulsive treatment: the interaction timescale $t_{{\rm int}}\equiv b/v_{{\rm enc}}$; the binary orbital period $P$; the merger time from a specific binary configuration via GW emission $T$; and the time between two consecutive encounters of the system and a fly-by perturber, $t_{{\rm enc}}=1/f=\left(n_{*}\sigma v_{{\rm enc}}\right)^{-1}$ where $n_{*}$ is the stellar number density, $\sigma$ is the geometric cross-section of the binary and the stellar fly-by.

We restrict our model to the impulsive regime, namely $t_{{\rm int}}\ll P$. This gives upper bound to the closest approach distance $b$ and hence limits the average time between encounters. For example, for a BBH with a total mass of $20M_{\odot}$, SMA of $a\sim 10^{4}{\rm AU}$ (hence $P\approx 3\times 10^{5}{\rm yr}$) and a typical velocity encounter of $v_{{\rm enc}}=50{\rm kms^{-1}}$(velocity dispersion in the field) we can restrict $b$ such that $t_{{\rm int}}=P/10$ $\left(P/100\right)$. Hence we get $b=t_{{\rm int}}\times v_{{\rm enc}}=3\times 10^{5}{\rm AU}\;\left(3\times 10^{4}{\rm AU}\right)$. Farther out flybys can also perturb the system and further excite the system, but at very large separations the interaction become adiabatic and the effects become small. We neglect the intermediate regime in which the perturbation time and the orbital times are comparable, which are likely to somewhat enhance the perturbation rates explored here.

II.2 Analytic Description

We consider the evolution of an ensemble of wide BBH binaries with initial separations $a>10^{3}{\rm AU}$ and comparable component masses $m_{1}\sim m_{2}=m_{{\rm BH}}$. For simplicity we assume all binaries to have the same SMA, and a thermal distribution of orbital eccentricities, $f(e)de=2ede$. In the following we derive the fraction of merging systems within this ensemble and find its dependence on the SMA of the binaries, $a$ and the environmental conditions, namely the stellar density $n_{*}$ and velocity dispersion $\sigma_{v}=v_{{\rm enc}}$.

The timescale for a GW-merger of an isolated binary is given by Peters (1964)

[TABLE]

with

[TABLE]

where $G$ is Newton’s constant and $c$ is the speed of light and $m_{1}=m_{2}=m_{{\rm BH}}$. Given a binary with SMA $a$ we can solve equation (1) for the critical eccentricity $e_{c}$ required for the binary to merger within some merger time ${\rm T}=t_{{\rm merge}}$

[TABLE]

All systems with eccentricities equal or greater than $e_{c}$ would therefore merge within this time-frame $T$. Hence, given a thermal distribution of eccentricities we find the fraction of system that merge within a time $T$ and lost from the ensemble to be:

[TABLE]

Following previous studies we term this “loss” region the “loss-cone” (see e.g. (Hills, 1981)); after time $T$ all binaries in the loss-cone merge via GW-emission and this phase-space region become empty. However, binaries outside the loss-cone which do not merge within this timescale, can be perturbed by a flyby encounter as to change their angular momentum, and thereby enter and replenish the loss cone. The average size of the phase-space region into which stars are perturbed during a single orbital period is termed the smear cone, defined by

[TABLE]

where $v_{k}$ is the Keplerian velocity of the binary. The value of $v_{k}$ can be calculated given that the average separation of a Keplerian orbit is $\left\langle r\right\rangle=a\left(1+1/2e^{2}\right)$ and we approximate $e\rightarrow 1$, namely

[TABLE]

$\left\langle\Delta v\right\rangle$ is the average change in the velocity over an orbital period due to perturbations (Hills, 1981). Let us consider fly-by interactions using the impulse approximation. Hills (1981) showed that on average the velocity change (for a binary with SMA, $a$), to the binary components is of the order of

[TABLE]

where $v_{{\rm enc}}$ is the velocity of the fly-by star with respect to the binary center of mass, $m_{p}$ is the perturber mass and $b$ is the closest approach distance of a fly-by. Therefore, the square of the angular size of the smear cone cause by the impulse of the fly-by on the binary is

[TABLE]

and for $\theta\ll 1$ we get the fractional size of the smear-cone velocity space over the $4\pi$ sphere to be after a single passage of the perturber

[TABLE]

It is evident from (6) that for a given binary the size of the smear cone depends on the perturber quantities, i.e. mass, velocity and the closest approach. The ratio of $F_{s}$ to $F_{q}$ indicates the fraction of the loss cone filled after a single fly-by.

[TABLE]

In the case where the loss cone is continuously full ($F_{s}=F_{q}$) the depletion rate only depends on the loss cone size, $F_{q}$ and the merger time , $T$. Hence the loss rate for the full lose cone is given by:

[TABLE]

Note that the loss rate is independent of the stellar density in the field, i.e. once the stellar density is sufficiently large as to fill the loss-cone, the loss rate is saturated, and becomes independent of the perturbation rate. Furthermore, one can see from Eq. (8) that the full loss-cone rate scales like $\dot{L}\propto F_{q}\propto a^{-8/7}$, i.e. the full loss-cone rate decreases with increasing SMA.

On the other hand, tighter binaries are less susceptible for change due to a fly-by, this is evident from equation (6). Therefore closer from a critical SMA we expect that the loss cone will not be full all the time, in this “empty loss cone” case the loss rate depends on the rate of orbits being kicked into the loss cone:

[TABLE]

where $n_{*}$ is the stellar density, $\sigma=\pi b^{2}$ is the geometric cross-section.

The condition for the loss cone to be continuously full is that the loss-cone orbits are replenished at least as fast as they are depleted due to the GW emission. This occurs when the rate of fly-by’s that enter orbits to the loss cone,$f$ is equal to the rate of which orbits are depleted from the loss cone, $1/T$:

[TABLE]

Furthermore, the condition for the loss cone to be continuously full is that the sizes of the lose-cone and the smear-cone are equal. This equilibrium occurs when

[TABLE]

A stellar fly-by is sufficiently strong as to replenish the loss cone if

[TABLE]

Plugging this to equation (10) we get an equation for the critical SMA that separates the empty and full loss-cone regimes:

[TABLE]

Using the critical SMA we can calculate the merger probability in each of these regimes, $a<a_{{\rm crit}}$ (empty) and $a>a_{{\rm crit}}$ (full). We denote $F_{q}$ as the fraction of wide binaries destroyed after time $T$, and therefore $\left(1-F_{q}\right)$ represents the fraction of binaries that survive as wide binaries at the relevant timescale. For the empty loss-cone regime ($a<a_{{\rm crit}})$ the relevant timescale is $1/f$; for the full loss-cone regime ($a>a_{{\rm crit}})$ the relevant timescale is $T$. Therefore, $\left(1-F_{q}\right)$ is a monotonically decreasing function of time, and the probability for a merger of a wide binary is

[TABLE]

where $t$ is the time since birth of the binary. As one can expect the probability only depends on the size of the loss cone and the rate of interactions. For the limit of $t\cdot f\cdot F_{q}\ll 1$ we can expand equation (14) and take the leading term, to find the loss probability to be approximated by

[TABLE]

given that

[TABLE]

together with equation (2) we get:

[TABLE]

In the full loss-cone regime the limiting factor is not the value of $f$, but the merger timescale $T$. Therefore, the full expression for the loss probability for $a>a_{{\rm crit}}$ is

[TABLE]

In the limit of $F_{q}\cdot t/T\ll 1$ we can approximate the probability by

[TABLE]

The above treatment neglects the fact that perturbations may also “ionize” a binary and destroy it, namely, the binary is disrupted by the random fly-bys. Such ionization process decreases the available number of wide binaries. To account for the ionization process we consider the finite lifetime of wide binaries due to fly-bys using the approximate relation given by (Bahcall et al., 1985) for $t_{1/2}$, the half-life time of a wide binary evolving through encounters

[TABLE]

Taking this into account we can correct for eq. (17) and eq. (19) to get

[TABLE]

and

[TABLE]

where $\tau=t_{1/2}/\ln 2$ is the mean-lifetime of the binary.

In order to estimate the number of systems observable within a year in aLIGO we first calculate the number of systems merging in a Milky-Way (MW)-like galaxy per unit time. In order to do that we need to integrate over all SMA in a given stellar density and over all stellar densities in the galaxy we model. We follow a similar calculation from Michaely and Perets (2016). We model the Galaxy in the following way, let $dN\left(r\right)=n_{*}\left(r\right)\cdot 2\pi\cdot r\cdot h\cdot dr$ be the the number of stars in a region $dr$ (and scale height $h$), located at distance $r$ from the center of the Galaxy. Following (Kaib and Raymond, 2014) and references within we model the Galactic stellar density in the Galactic disk as follows

[TABLE]

where $n_{0}=0.1{\rm pc^{-3}}$ is the stellar density near our sun, $R_{l}=2.6{\rm kpc}$ (Jurić et al., 2008) is the galactic length scale and $r_{\odot}=8{\rm kpc}$ is the distance of the sun from the galactic center. Integrating over the stellar densities throughout the Galaxy we can obtain the total number of mergers through this process. Next we account for the fraction of wide BBH systems from the entire population of starts in the Galaxy. We use the following standard values. Given a Kropa initial mass function (Kroupa, 2001), the fraction of the stars that evolve to become BHs is $f_{{\rm primary}}\approx 10^{-3}$. If we assume most BHs form without any natal-kick ( similar assumptions were taken in other works (e.g. Mandel, 2016; Belczynski et al., 2016)), we can expect all BHs to be in binary (or higher multiplicity) systems and the fraction $f_{{\rm bin}}=1$, consistent with the binary fraction inferred for the O-stars progenitors of BHs (Moe and Di Stefano, 2016; Duchêne and Kraus, 2013; Sana et al., 2014). Next we assume a uniform distribution of the mass ratios, $Q\in\left(0.1,1\right)$ (Moe and Di Stefano, 2016; Duchêne and Kraus, 2013), to get a fraction of secondaries that evolve into BHs of $f_{{\rm secondary}}\approx 0.4$. We also assume that the SMA has a log-uniform distribution (Opik law) and therefore the fraction of systems with SMA larger than $10^{3}{\rm AU}$ is $f_{{\rm wide}}\approx f_{{\rm bin}}\times 0.2$. This value is actually a lower limit for massive binaries, recently Igoshev and Perets (2019) found that a wide binary fraction for massive B-stars to be $f_{{\rm wide}}\approx 0.5$, and theoretical models suggest that the fraction of wide binary O-stars and BHs could be close to unity (Perets and Kouwenhoven, 2012), $f_{{\rm wide}}=1$, and we therefore expect a wide-binary fraction in the range $0.2-1$, in the following we use $f_{{\rm wide}}=0.5$.

[TABLE]

[TABLE]

The number of merging BBH per ${\rm Myr}$ from this channel for a MW-like Galaxy is

[TABLE]

Following Belczynski et al. (2016) we calculate the merger rate, $\mathscr{R}$ per ${\rm Gpc^{3}}$ by using the following estimate

[TABLE]

[TABLE]

while $\rho_{{\rm gal}}$ is local density of the MW-like galaxies with the value of $\rho_{{\rm gal}}=0.0116{\rm Mpc^{-3}}$(e.g. Kopparapu et al. (2008)) and $\Gamma$ is given in the units of ${\rm Myr^{-1}.}$

III Numerical calculation

In this section we describe the numerical calculation we preform. We simulate the evolution of a binary BH with masses of $m_{1}=m_{2}=10M_{\odot}$ for $10{\rm Gyr}.$ We treat the evolution by considering both the evolution of the binary between encounters, and in-particular the effects of GW-emission, as well as the change of the binary orbital elements due to the impulsive fly-by encounters with field stars.

In order to calculate the average time between encounters we use the rate $f=n_{*}\sigma\left\langle v_{{\rm enc}}\right\rangle$, where $n_{*}$ is the stellar number density, taken to be the solar neighborhood value of $n_{*}=0.1{\rm pc^{-3}}$; $\left\langle v_{{\rm enc}}\right\rangle$ is the velocity of the perturber as measured from the binary center of mass, where we set $\left\langle v_{{\rm enc}}\right\rangle=50{\rm kms^{-1}}$ similar to the velocity dispersion in the solar neighborhood; and $\sigma$ is the interaction cross-section. We focus on the impulsive regime, namely $t_{{\rm int}}\ll P$ (see subsection II.1). With these values the largest closest approach distance $b$ for which an encounter can be considered as impulsive is $b_{{\rm max}}=5\times 10^{4}{\rm AU}.$** **The average time between such impulsive encounters is given by $t_{{\rm enc}}=1/f\approx 1{\rm Myr}$. Therefore we randomly sample the time between encounters from an exponential distribution with a mean $f$ (due to the Poisson distribution of encounter times).

We initialize the wide binary with a SMA $a$ and eccentricity $e$. At each step we first find the next encounter time $t_{{\rm enc}}$, and evolve the binary for $t_{{\rm enc}}$ through the equations of motion given by Peters (1964),

[TABLE]

[TABLE]

where $G$ is Newton’s constant and $c$ is the speed of light. If the binary did not merge through GW-emission by the time of the next encounter we simulate the impulsive interaction with a perturber with velocity of $v_{{\rm enc}}$ drawn from a Maxwellian distribution with velocity dispersion $\left\langle v_{{\rm enc}}\right\rangle$ and a mass of $m_{p}=0.6M_{\odot}$, typical for stars in the field. After changing the binary orbital parameters due to the encounter we continue to evolve the binary until the next encounter and so on, until the binary merges, disrupts or the maximal simulation time of 10 Gyrs is reached.

The numerical results are presented in Figure 2; the numerical result are highly consistent with the result of the analytic model. A more detailed technical description of the numerical procedure, equations and analysis is given in the appendix Appendix I and Appendix II.

IV Discussion and summary

In this paper we explore a novel channel for the production of BBH GW-sources from wide ($>10^{3}$AU) binaries in the field. Such binaries are sensitive to perturbations by stellar-fly-bys even in the low density environment in the field. We find that a fraction of all wide-binaries attain sufficiently close-approaches (as their orbits are excited to very high eccentricities) as to inspiral and eventually merge through GW-emmision before any consecutive encounter can change the orbit.

The merger rate strongly depends on the natal-kicks given to BHs at birth, which are poorly constrained (e.g. Repetto et al., 2012, 2017). In particular, it is still unknown whether a BH receives a momentum kick at birth like a NS, or forms without any natal-kick following a failed supernova or a large amount of fallback (e.g. Ertl et al., 2015; Belczynski et al., 2004, 2008). Previous models that were able to produce rates comparable to the rate inferred from observations had typically taken similar assumptions of zero kick velocities (for all BHs, or at least for all BHs more massive than 10 ${\rm M}_{\odot}$), while models assuming higher natal kicks produced significantly lower rates (Belczynski et al., 2008; Dominik et al., 2015). In our case low natal-kicks can unbind the wide binaries, lowering their fraction. In principle, in models where wide-binaries form following the dispersal of their birth-cluster on longer time-scales (Kouwenhoven et al., 2010; Perets and Kouwenhoven, 2012), BH may acquire wide companions well after their formation. Nevertheless, even in these cases the BHs need to be retained in the cluster until its dispersal, and therefore the natal kick needs to be sufficiently low for a BH not to escape the cluster. We conclude that adapting similar no-kick assumptions for BHs (as done by other potentially successful scenarios) suggests the wide-binary channel explored here can give rise to a high production rate of GW-sources from BBH mergers of perturbed ultra-wide binaries.

Beside the rate estimate, $\sim 10\times f_{{\rm wide}}{\rm\,yr^{-1}}{\rm Gpc^{-3}}$, our proposed evolutionary channel gives rise to specific characteristics of the BBH mergers, which together can provide a distinct signature for this channel, as we discuss in the following.

Eqs. (21) and (22) describe the probability dependency for a given environment, namely the stellar density $n_{*}$ and the encounter velocity $v_{{\rm enc}}.$ We note the in both equations there is a $\tau\left(1-e^{-t/\tau}\right)$ dependency. Hence following eq. (20) the merger probability increases with the encounter velocity. For example, taking the same environment as assumed in section II.2 but with $v_{{\rm enc}}=200{\rm kms^{-1}}$ gives a factor of $\sim 1.92$ higher rate of BBH GW-sources. We therefore expect a preference for host galaxies with higher velocity dispersion.

This model is sensitive to extreme mass ratio. The equations that govern the rates depend on $\beta^{\prime}$. When the binary mass is kept constant but the mass ratio $Q$ is varied we get the following dependence on $Q$

[TABLE]

and since $Q$ is a monotonically increasing function, equal mass components have the highest probability to merge. Moreover, the merger rate also has a monotonic dependence on the total binary mass, due to the complex mass dependence in the loss-cone analysis and the effects of ionization (see Eqs. 21 and 22). Hence overall we expect a preference towards GW-sources from more massive binaries and higher mass-ratios. Furthermore, in this channel the spins of the BHs are likely uncorrelated given the origin of the BH components from very wide separations (or a random capture) and we therefore expect the spins of the merging BBH components to be randomly (mis-)alligned, in contrast with e.g. the isolated binary evolution channel; (Mandel and de Mink, 2016, e.g.). Given the long time-scale for inspiral from large separations we also expect BBHs to fully circularize by the time they reach the aLIGO band and to not produce any eccentric binaries at these frequencies, in contrast with some of the dynamical channels. Finally, unlike the isolated binary channel, which predicts a delay time dependence of $\propto t^{-1}$ (Dominik et al., 2015) our model, which have no time dependency on the merger probability, generally predicts a uniform delay time distribution.

In summary, the wide-binary origin for BBH GW mergers can give rise to a potential rate of $\sim 10\times f_{{\rm wide}}$ {\rm yr}^{-1}$${\rm Gpc^{-3}} (where $f_{{\rm wide}}$ can plausibly reside in the range $0.2-1$), comparable to (the lower range of) the observationally inferred rate of $\sim 10-110$ {\rm yr}^{-1}$${\rm Gpc^{-3}} from aLIGO/VIRGO detection, and is strongly dependent on the natal-kicks imparted to BHs at birth. It can be characterized by the following signatures: (1) A slight preference for high mass ratio BBH GW-sources. (2) A preference for more massive BBH. (3) Typically randomly misaligned spin-orbits BHs. (4) Circular orbits in the aLIGO band. (5) Preference for high velocity dispersion host galaxies/environments. (6) A uniform delay-time distribution. We acknowledge support from the ISF-ICORE grant 1829/12. The authors acknowledge the University of Maryland supercomputing resources (http://hpcc.umd.edu) made available for conducting the research reported in this paper. EM would like to thank Coleman Miller and Johan Samsing for stimulating discussions regarding this work.

Appendix I

Post-interaction

orbital elements

In order to calculate the post-interaction orbital elements of a pertubed binary, we apply the following procedure. We first randomize the perturber trajectory; we randomly sample the closest approach point, $\vec{b}=b\hat{b}$ of the perturber from an isotropic distribution by randomizing the two spherical coordinate angles, uniformly distributed $\cos\alpha\in\left(-1,1\right)$ and uniformly distributed $\beta\in\left(0,2\pi\right)$, which determine a plane perpendicular to the closest approach vector $\vec{b}$. The perturber trajectory can have any direction within this plane, hence we randomize an additional angle, uniformly distributed in the range $\xi\in\left(0,2\pi\right)$ as to choose an arbitrary axis in the plane. Next we randomize $b$ by setting the distribution function to be $f\left(b\right)\propto b_{{\rm max}}$.

We randomize the state of the binary by randomizing its mean anomaly, ${\scriptscriptstyle M}$ from a uniform distribution between $\left(0,2\pi\right)$ and calculate the post-interaction binary orbital elements, $a_{{\rm post}}$ and $e_{{\rm post}}$ using the impulse approximation. In the impulse approximation we neglect any motion of the binary during the passage of perturber. In this approximation we can calculate the velocity change of each of the components of the binary [Rickman, 1976, Collins and Sari, 2008]. The velocity vector change for $m_{1}$ is given by

[TABLE]

where $\vec{r}_{p}$ is the position of the fly-by perturber set to be at the closest approach at $t=0$. In the impulse regime we can approximate the trajectory of the fly-by during the interaction as a straight line thus

[TABLE]

The velocity change for $m_{2}$ is then

[TABLE]

where $b_{2}$ is the closest approach of the fly-by perturber. We can find the relation between $\vec{b}$ and $\vec{b}_{d}$ in the straight line trajectory approximation by [Collins and Sari, 2008]:

[TABLE]

The change in the relative velocity is simply

[TABLE]

Given the relative velocity vector change we can calculate the post-interaction eccentricity vector

[TABLE]

to find the post-interaction eccentricity which is the norm of the eccentricity vector $e_{{\rm post}}=\left|\vec{e}\right|$.

For the post-interaction SMA we calculate the change in orbital energy. The energy change is due to the velocity kick imparted by the perturber. In the impulse approximation we model the interaction only via a velocity change and not by a change in the separation itself during the encounter. The specific orbital energy is

[TABLE]

where $r$ is the separation. The velocity change $\vec{\Delta v}_{r}$ can be written as a sum of the parallel (to the instantaneous orbital velocity) and the perpendicular vectors $\vec{\Delta v}_{r}=\Delta v_{r,\parallel}+\Delta v_{r,\perp}$. Hence the specific energy change is given by

[TABLE]

which translates to change in the SMA of

[TABLE]

to give us the final post-interaction SMA $a_{{\rm post}}=a+\Delta a$.

Appendix II

Numerical Results

In order to validate our analytic calculation we numerically verify equations (21) and (22) as plotted in Figure 1. We calculate the merger probability for the same binary set up and stellar environment, specifically for a BBH with $m_{b}=20M_{\odot}$ pertubed by stars of typical mass $m_{p}=0.6M_{\odot}$ and velocity dispersion of $\left\langle v_{{\rm enc}}\right\rangle=50{\rm kms^{-1}}$ the stellar number density is taken to be $n_{*}=0.1{\rm pc^{-3}}$. The final merger probability is calculated after $t=10{\rm Gyr}$ since the BBH was initialized. We calculate the GW merger probability for several SMA values, $a\approx a_{{\rm crit}}=2\times 10^{4}{\rm AU}$ and $a=3\times 10^{4}{\rm AU},\ 1\times 10^{4}{\rm AU}.$ For each value of $a$ we consider a range of initial eccentricities $e$, and follow the evolution of each of the systems for $10{\rm Gyr}$. In Table 1 we present the number of systems for each combination of $a$ and $e$. For each of the modeled initial conditions we calculate the fraction of merged systems in the following way. We record the number of merged systems $n\left(a,e\right)$ to find the fraction of merged systems out of the total number of modeled systems

[TABLE]

We can then also find the standard deviation given by

[TABLE]

In Figure 3 we present the estimated probability as a function of initial eccentricity for the case of $a\approx a_{{\rm crit}}=2\times 10^{4}{\rm AU}.$

In order to calculate the overall probability, $P$ we weigh $f_{{\rm merged}}$ with a thermal distribution of the initial eccentricities. Analytically this corresponds to the following integral

[TABLE]

where $f\left(e\right)=2e$ is the thermal distribution of the initial eccentricities. Numerically we approximate the integral by the following sum,

[TABLE]

where $\Delta e_{i}=e_{i+1}-e_{i}$ and $\max\left\{e_{i}\right\}<1$. Another approximation is

[TABLE]

where $\Delta e_{i}=e_{i+1}-e_{i}$ and $\max\left\{e_{i}\right\}=1$. It is clear that $P_{1}$ overestimates the integral while $P_{2}$ underestimates $P$. Hence we take the arithmetical mean of $P_{1}$ and $P_{2}$ and assign it as the calculated probability, $P$.

In order to calculate the standard deviation we use the following estimate

[TABLE]