The Effect of Distant Encounters on Black Hole Binaries in Globular Clusters: Systematic Increase of In-cluster Mergers in the LISA band

TL;DR

This study systematically examines how weak, distant encounters influence the evolution and merger rates of binary black holes in dense stellar clusters, revealing that such interactions significantly increase in-cluster mergers and eccentric LISA sources.

Contribution

It introduces the first Monte Carlo simulations that incorporate both weak and strong interactions to model BBH evolution in clusters, highlighting the importance of weak interactions.

Findings

Weak interactions increase in-cluster BBH mergers.

More eccentric LISA sources are produced with weak interactions.

Including weak interactions alters GW inspiral times and merger probabilities.

Abstract

We present the first systematic study on how distant weak interactions impact the dynamical evolution of merging binary black holes (BBHs) in dense stellar clusters. Recent studies indicate that dense clusters are likely to significantly contribute to the rate of merging BBHs observable through gravitational waves (GWs), and that many of these mergers will appear with notable eccentricities measurable in the LISA and LIGO sensitivity bands. This is highly interesting, as eccentricity can be used to distinguish between different astrophysical merger channels. However, all of these recent studies are based on various Monte Carlo (MC) techniques that only include strong interactions for the dynamical evolution of BBHs, whereas any binary generally undergoes orders-of-magnitude more weak interactions than strong. It is well known that weak interactions primarily lead to a change in the binary's eccentricity, which for BBHs implies that weak interactions can change their GW inspiral time and thereby their merger probability. With this motivation, we perform MC simulations of BBHs evolving in dense clusters under the influence of both weak and strong interactions. We find that including weak interactions leads to a notable increase in the number of BBHs that merge inside their cluster, which correspondingly leads to a higher number of eccentric LISA sources. These preliminary results illustrate the importance of including weak interactions for accurately modeling how BBHs merge in clusters, and how to link their emitted GW signals to their astrophysical environment.

Equations

$$ϵ={\left[\frac{M^2}{m_b(m_b+M)}{\left(\frac{a_0}{r_{\mathrm{p}}}\right)}^3(1+E{)}^{-3}\right]}^{1/2},$$

$${\Delta e}_{\mathrm{FO}}=-ϵ\frac{15\pi }{4}{e}_0\sqrt{1-e_0^2}\sin 2\Omega {\sin }^{2}i,$$

$$\tau \approx \frac{768}{425}\frac{5c^5}{512G^3}\frac{a^4}{m^3}{\left(1-e^2\right)}^{7/2},$$

$$t_{\mathrm{SE}}\approx \frac{1}{12\pi G}\frac{v_{\mathrm{dis}}}{n_{\mathrm{s}}ma},$$

$$e_{\mathrm{SE}}\approx \sqrt{1-{\left(t_{\mathrm{SE}}/{\tau }_{\mathrm{c}}\right)}^{2/7}},$$

$${\Delta e}_{\mathrm{SO}}=$$

$$\displaystyle\biggl{\{}\left(5\cos i+3\cos 3i\right)\cos 2\omega+6\sin i\sin 2i\biggl{\}}$$

$$\displaystyle+4\cos 2i\biggl{\{}3\pi\left(81e_{0}^{2}-56\right)-200\left(1-e_{0}^{2}\right)$$

$$\displaystyle\cos 2\Omega\sin 2\omega\biggl{\}}+3\pi\biggl{\{}200e_{0}^{2}\sin^{4}i\cos 4\Omega$$

$$+8\left(16e_0^2+9\right){\sin }^{2}2i\cos 2\Omega$$

$$\displaystyle+\left(39e_{0}^{2}+36\right)\cos 4i-299e_{0}^{2}+124\biggl{\}}\biggl{]},$$

$$\frac{d{\Delta e}_{\mathrm{tot}}}{dt}\propto {\int }_{r_{\mathrm{p},\min }}^{r_{\mathrm{p},\max }}\frac{d{\Gamma }_{\mathrm{s}}}{d{r}_{\mathrm{p}}}{\Delta e}_{+}d{r}_{\mathrm{p}}\propto \frac{1}{\sqrt{r_{\mathrm{p},\min }}}-\frac{1}{\sqrt{r_{\mathrm{p},\max }}}.$$

$$(r_{\mathrm{p}}/a_0{)}_{\mathrm{FO}}^{\mathrm{BR}}\approx {\left[\frac{e_0\sqrt{1-e_0^2}}{e_0-1}\frac{15\pi \sin 2\Omega {\sin }^{2}i}{16\sqrt{3}}\right]}^{2/3}.$$

$$\frac{(r_{\mathrm{p}}/a_0{)}_{\mathrm{plateau}}}{(r_{\mathrm{p}}/a_0{)}_{\Delta e=0}}=2^{2/3},$$

$$T_{1PN}\approx \frac{T_0}{3}\frac{a_0}{{\mathcal{R}}_{\mathrm{m}}}(1-e_0^2),$$

$$T_3\approx \sqrt{\frac{r_{\mathrm{p}}^3}{Gm}},$$

$$\frac{r_{\mathrm{p}}}{a_0}$$

$$\gtrsim 10^3\times {\left(\frac{a_0}{0.5\text{AU}}\right)}^{2/3}{\left(\frac{m}{20M_{\odot }}\right)}^{-2/3}{\left(1-{\left(e_0/0.99\right)}^2\right)}^{2/3},$$

$$f_{\mathrm{GW}}\approx \frac{1}{\pi }\sqrt{\frac{2Gm}{r_{\mathrm{BH}}^3}}.$$

$$\frac{f_{\mathrm{GW}}}{\text{Hz}}\approx 2\cdot 10^{-5}{\left(\frac{\mathcal{T}}{10^{10}\text{yrs}}\right)}^{-3/7}{\left(\frac{a}{0.5\text{AU}}\right)}^{3/14}{\left(\frac{m}{30M_{\odot }}\right)}^{-11/14}$$

$$a_{\mathrm{HB}}\approx \frac{3}{2}\frac{Gm}{v_{\mathrm{dis}}^2},$$

$$a_{\mathrm{enc}}=a\times \delta ,$$

$$a_{\mathrm{ej}}\approx \frac{1}{6}\left(\frac{1}{\delta }-1\right)\frac{Gm}{v_{\mathrm{esc}}^2},$$

$$\Delta E_{\mathrm{GW}}\approx (85\pi /12)G^{7/2}{c}^{-5}{m}^{9/2}{r}_{\mathrm{BH}}^{-7/2},$$

$$E_{\mathrm{B}}\approx G{m}^2/(2a_0).$$

$$r_{3cap}\approx {\mathcal{R}}_{\mathrm{m}}\times {\left(a_0/{\mathcal{R}}_{\mathrm{m}}\right)}^{2/7}.$$

$$P_{3m}\approx \frac{2r_{3cap}}{a_0}\times N_{\mathrm{IMS}}.$$

$$t_{\mathrm{int}}\approx \frac{1}{6\pi G}\frac{v_{\mathrm{dis}}}{n_{\mathrm{s}}m{R}_{\mathrm{p}}}.$$

Full text

I Introduction

The origin of binary black hole (BBH) mergers is still unknown, despite several now have been directly observed through their emission of gravitational waves (GWs) (2016PhRvL.116f1102A, ; 2016PhRvL.116x1103A, ; 2016PhRvX…6d1015A, ; 2017PhRvL.118v1101A, ; 2017PhRvL.119n1101A, ; 2019arXiv190210331Z, ; 2019arXiv190407214V, ). The currently observed set shows great variety in both the BH mass spectrum, and BH spins (e.g. 2018arXiv181112907T, ). At least one neutron star (NS) merger has also been observed (2017PhRvL.119p1101A, ) in addition to several peculiar low signal-to-noise events, and this variety therefore indicates that several formation channels might be at play. Some of the recently proposed channels include: 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, ), 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, ), active galactic nuclei (AGN) discs (2017ApJ…835..165B, ; 2017MNRAS.464..946S, ; 2017arXiv170207818M, ), galactic nuclei (GN) (2009MNRAS.395.2127O, ; 2015MNRAS.448..754H, ; 2016ApJ…828…77V, ; 2016ApJ…831..187A, ; 2016MNRAS.460.3494S, ; 2017arXiv170609896H, ; 2018ApJ…865….2H, ), very massive stellar mergers (Loeb:2016, ; Woosley:2016, ; Janiuk+2017, ; DOrazioLoeb:2017, ), and single-single GW captures of primordial black holes (2016PhRvL.116t1301B, ; 2016PhRvD..94h4013C, ; 2016PhRvL.117f1101S, ; 2016PhRvD..94h3504C, ).

The key exercise is how to observationally tell these channels apart; several observables have been suggested. For example, an electromagnetic (EM) signal might be associated with both BBH mergers in AGN discs and those coming from rapidly rotating stars, in contrast to gas-free isolated BBH mergers found in, e.g., stellar clusters and the field. Tidal disruption events (TDE) associated with stars disrupted by BBHs in dense stellar clusters might also be used to probe the merger history of BBHs (2019ApJ…877…56L, ; 2019arXiv190102889S, ; 2019arXiv190406353K, ). Other ways of ‘lightning up’ single or binary BHs through TDEs include, e.g., BHs kicked into disruptive orbits from either dynamics, GW recoils or direct collapse following supernova (e.g. 2017PhRvL.118o1101S, ; 2019arXiv190509471F, ). Such multi-messenger (EM+GW) observables are interesting, but unfortunately also highly uncertain to both accurately predict and model. More promising are properties that can be read off directly from the incoming GW waveform, which include the BBH masses and spins, the orbital eccentricity, and even doppler effects related to a possible movement of the BBH’s center of mass (COM) (e.g. 2017ApJ…834..200M, ; 2018arXiv180505335R, ). For example, the spins are expected to be isotropically distributed for BBH mergers forming in clusters, in contrast to field BBH mergers (e.g. 2016ApJ…832L…2R, ; 2017arXiv170601385F, ). BBHs assembled in clusters through chaotic few-body interactions have also recently been shown to lead to a non-negligible unique population of mergers with measurable eccentricities in all bands from LISA (the Laser Interferometer Space Antenna) (2018MNRAS.tmp.2223S, ; 2018MNRAS.481.4775D, ; 2019PhRvD..99f3003K, ) to LIGO (the Laser Interferometer Gravitational Wave Observatory) (2006ApJ…640..156G, ; 2014ApJ…784…71S, ; 2017ApJ…840L..14S, ; 2017ApJ…840L..14S, ; 2018MNRAS.476.1548S, ; 2018ApJ…853..140S, ; 2019MNRAS.482…30S, ; 2018PhRvD..97j3014S, ; 2018ApJ…855..124S, ; 2019ApJ…871…91Z, ; 2018PhRvD..98l3005R, ). Eccentric sources can form in other ways, e.g. through Lidov-Kozai oscillations (e.g. 2011ApJ…741…82T, ; 2018ApJ…864..134R, ; 2018arXiv181110627F, ; 2019arXiv190208604R, ; 2019arXiv190309160F, ; 2019MNRAS.486.4443F, ; 2019arXiv190309659F, ), quadruple systems (e.g. 2019MNRAS.486.4781F, ), and single-single GW captures (e.g. 2009MNRAS.395.2127O, ; Kocsis:2012ja, ; 2016PhRvD..94h4013C, ; 2018ApJ…860….5G, ), but eccentric cluster mergers are unavoidably produced if BHs are present in clusters, in particular globular clusters (GCs), and therefore constitute a highly useful diagnostic for the cluster channel.

The topic of this paper is on how BBHs are dynamically driven to merger in GCs under the influence of both weak and strong interactions, as well as what unique GW signals to expect from this population. BBHs in GCs generally form through single-single-single interactions (e.g. Heggie:1975uy, ; 1976A&A….53..259A, ; Hut:1983js, ) in the cluster center, where the mass segregated BH population constitute a highly dense BH subsystem (e.g. 2018MNRAS.478.1844A, ; 2018MNRAS.479.4652A, ). Each BBH formed in this way is initially very wide and will not merge within a Hubble time, but subsequent interactions with incoming single BHs lead to shrinkage of the orbit, referred to as hardening, and large changes in the eccentricity (e.g. Heggie:1975uy, ). This eventually drives the BBH to merger either in-side or out-side of its cluster. Current state-of-the-art Monte Carlo (MC) cluster codes aimed at modeling the evolution of such BBHs, including the CMC code (2018PhRvD..98l3005R, ) and the MOCCA code (2013MNRAS.431.2184G, ), incorporate only the effects from strong encounters, and include therefore not the potential importance of the many more weak encounters a given BBH undergoes during its life inside its cluster.

In this work we present the first study on how the numerous number of unbound weak encounters a BBH undergoes in-between each strong encounter affects its dynamical evolution and GW merger properties. Each of these weak encounters changes the eccentricity of the BBH, which leads to a non-linear and asymmetric random walk with varying step sizes, also known as a quasi-Levy flight, in-between each strong encounter. In classical cluster dynamics where post-Newtonian (PN) corrections (e.g. 2014LRR….17….2B, ) are not included and all objects are assumed to be point-like, weak encounters are not important as the associated energy exchange is exponentially suppressed; however, when dissipative effects, such as PN corrections, tides, and finite sizes are included, objects from ordinary stars to BHs can be driven to merger through high eccentricity orbits (e.g. 2017ApJ…846…36S, ). Such eccentricity driven mergers not only give rise to interesting observables, including GWs, TDEs, and jetted EM sources, but will also change the energy balance of the cluster and thereby its large scale dynamics.

We here explore the combined effects from strong and weak three-body encounters using analytical and numerical methods, and investigate how their interplay drives BBHs to merger in dense clusters. For this we particularly make use of the strong scattering model presented in (2018PhRvD..97j3014S, ) and (2018MNRAS.tmp.2223S, ), and our new second-order weak interaction solution presented in (2019arXiv190409624H, ), referred to as HS19, which expands upon the first-order solution presented in (1996MNRAS.282.1064H, ), referred to as HR96. We find that the inclusion of weak interactions leads to notable changes in both the number of BBHs merging inside their cluster and their distribution across GW peak frequency and orbital eccentricity. This leads us to conclude that weak encounters are important to include for an accurate modeling of how BBHs merge in clusters, which is naturally an increasingly relevant inquiry now that there is a growing population of observed GW events. This poses some interesting challenges, as the leading numerical schemes which include the Hénon method in the field of PN dynamics do not easily allow for the inclusion of weak encounters. We also highlight the importance of using the second-order results from HS19, which solve some crucial problems occurring when using the expressions of HR96 for evolving BBHs in the high eccentricity regime.

The paper is organized as follows. We start by describing the outcome of a single weak encounter in Section II. This includes a discussion of why distant encounters must be important for an accurate modeling of BBH mergers, and several new dynamical effects that have not yet been properly described in the literature related to this subject. After this, we study in Section III a population of BBHs and how they evolve under the influence of both strong and weak encounters. For this we use numerical and analytical techniques combined with MC methods. We discuss and conclude our findings in Section IV.

II A Single Weak Encounter

We start by exploring the range of outcomes associated with a single weak encounter between a BBH and a single incoming BH. For this, we consider a BBH with initial orbital eccentricity $e_{0}$ and semi-major axis (SMA) $a_{0}$ interacting with an incoming BH on a parabolic orbit. Here, and in the rest of the paper, we assume the interacting BHs all have the same mass $m$. For describing the ICs we make use of the standard orbital angles shown in Figure 1, which is further described in HR96. Note here that these are different from the set used in HS19 (see Section 2.1 in the latter paper).

II.1 Background and Motivation

For binary-single interactions to be in the secular regime, the mean motion of the binary must be much greater than the orbital speed of the single at the time of peri-center passage (HR96). In this limit, the binary-single equations-of-motion (EOM) can be written out using a technique known as ‘orbital averaging’, which is widely used for e.g. studying the hierarchical bound-triple Lidov-Kozai problem (e.g. 2016ARA&A..54..441N, ). In the case where the third object is on an unbound orbit with an initial eccentricity $E$, the orbital averaged EOM naturally gives rise to the following parameter $\epsilon$, which appears as a common factor in the EOM (HS19),

[TABLE]

where $m_{b}$ is the mass of the BBH ($m_{b}=2m$), $M$ is the mass of the perturber ($M=m$), and $r_{\rm p}$ is the distance between the BBH center-of-mass (COM) and the incoming BH at the time of peri-center passage. As seen, weak interactions for which $r_{\rm p}/a{{}_{0}}\gg 1$ are characterized by having $\epsilon\ll 1$.

The main effect from a BBH undergoing a weak interaction with a single incoming BH is a change in the orbital eccentricity of the BBH, $\Delta{e}$, such that the final eccentricity after the encounter is $e_{\rm fin}=e_{0}+{\Delta{e}}$. The change $\Delta{e}$ can to first-order (FO) in $\epsilon$ in the parabolic limit be analytically expressed as (HR96),

[TABLE]

where $\Omega$ is the longitude of the ascending node measured in the plane of motion of the binary, and $i$ is the inclination between the two orbital planes of the binary and the perturber, respectively, as further described in Fig. 1. Note here that this term has been derived by expanding the EOM to lowest non-trivial order in the (assumed to be) small ratio of the binary separation to the distance between the binary COM and the perturber, i.e., it is only accurate to quadrupole order. As shown in HR96, this expression accurately fits numerical simulations for $r_{\rm p}/a_{0}\gg 1$ and intermediate values of $e_{0}$, and have therefore been widely used as a replacement for otherwise computationally demanding three-body simulations (e.g. 2019ApJ…872..165G, ). HR96 also presented solutions to the $e_{0}=0$ limit, and in regions where the secular approach breaks down. However, one consideration not addressed in HR96 is the regime where ${\Delta{e}}\approx 1-e_{0}$, which is possible if $e_{0}$ is near unity. For describing the formation of BBH mergers in dense clusters this is a highly relevant limit, as basically all mergers form through high eccentricity orbits. It is known that BBHs can be driven to very high eccentricities entirely through strong few-body encounters, i.e. without any assistance from weak interactions (see e.g. (2018PhRvD..97j3014S, )); however, binaries undergoing strong interactions will also undergo numerous weak interactions in-between each strong interaction each of which will perturb the binary (e.g. 2019ApJ…872..165G, ). The question is if a model that also includes weak interactions will lead to more or less mergers from high eccentricity orbits.

To gain insight into this and motivate our present study, we now consider the following simple cluster model: we consider a cluster core described by a BH subsystem with a constant density and velocity dispersion, and assume that all BH interactions take place in this system and only involve objects of similar mass. This model provides in fact a good description of real cluster cores, as BHs both are expected to segregate, separate in mass, and clear out the center, as a result of mass segregation. In this picture, if only strong binary-single encounters are included in the cluster dynamics, a BBH will merge in-between its interactions if its GW inspiral life time (Peters:1964bc, ),

[TABLE]

is less than its strong encounter (SE) time, $t_{\rm SE}$, defined as the time it takes for the BBH to undergo its next strong interaction (e.g. 2018MNRAS.tmp.2223S, ),

[TABLE]

where $n_{\rm s}$ is the number density of single BHs in the cluster, $v_{\rm dis}$ is the cluster velocity dispersion, $a$ is the SMA of the BBH, and $e$ is its eccentricity. This correspondingly implies that a BBH with an eccentricity $e>e_{\rm SE}$,

[TABLE]

where $\tau_{\rm c}$ denotes the time for which $e=0$ in Eq. (3), will undergo a merger before its next strong binary-single interaction. For a BBH with $a=1$ AU, $m=20\ M_{\odot}$, $n_{\rm s}=10^{5}\ \text{pc}^{-3}$, and $v_{\rm dis}=10\ \text{kms}^{-1}$, one finds $e_{\rm SE}\approx 0.99$. This is a relatively high eccentricity, which means that even small changes of order ${\Delta{e}}\sim 1-0.99\sim 10^{-2}$ can change the state and outcome of the BBH, from e.g. surviving to merging. Now, changes of such magnitude can easily be provided by weak interactions; the question is then, will weak interactions on average prevent or help a BBH to merge if its eccentricity is already near the critical value $e_{\rm SE}$? As pointed out above, the FO result ${\Delta{e}}_{\rm FO}$ presented in HR96 breaks down in this high eccentricity limit, and will therefore often lead to unphysical results such that $e_{\mathrm{fin}}\geq 1$.

To resolve this problem with the FO solution in the high eccentricity limit, we used in our companion paper HS19 second-order (SO) perturbation theory to derive higher order correction terms to ${\Delta{e}}$. In short, in the FO approximation of HR96, the binary orbital components are assumed to be fixed while integrating over the passage of the perturber (this is analogous to ‘double averaging’ in hierarchical triple systems). Instead, in HS19 we used Fourier expansions to take into account the instantaneous response of the binary to the perturber while integrating over the latter’s passage. This resulted in an expression for the eccentricity change to SO in $\epsilon$, which in the limit of small eccentricity changes and parabolic perturbers reads,

[TABLE]

where the introduced set of angles $\{i,\Omega,\omega\}$ are shown and defined in Fig. 1. We reiterate that the orbital angles used here are defined differently than in HS19. As for ${\Delta{e}}_{\rm FO}$, this term is only accurate to quadrupole order. The second term that is SO in $\epsilon$, resolves many problems associated with only using the FO term ${\Delta{e}}_{\rm FO}$, as further explained in HS19. Note here that all dependencies on mass and closest approach only enter via the $\epsilon$ term. A few examples of ${\Delta{e}}_{\rm SO}$, as a function of the peri-center distance of the perturber w.r.t to the COM of the BBH scaled by the initial SMA, $r_{\rm p}/a_{0}$, are shown and discussed in Fig. 3. As seen, the SO correction clearly breaks the simple FO power-law solution. This will be described in greater detail later. Finally, one should note that the SO solution is still only accurate to quadrupole order and will therefore break down when $r_{\rm p}/a_{0}$ starts to approach unity, and/or when the binary’s component masses are unequal. In contrast, the FO solution can break down at any $r_{\rm p}/a_{0}$, i.e. also in highly secular regions, depending on the initial eccentricity. The improvements from including octupole terms (the next non-trivial-order terms in the expansion of the small ratio of the binary separation to the distance between the binary COM and the perturber), as well has higher-order terms in $\epsilon$ (i.e,. third order), will be explored in upcoming work.

In the subsections below we systematically describe the outcomes of a BBH undergoing a single weak encounter, and the associated novel features that are not captured by the classical FO result ${\Delta{e}}_{\rm FO}$ from HR96. All our results are based on the initial set of orbital angles $i=\pi/2$, $\Omega=-\pi/4$, $\omega=0$, a set we loosely refer to as our ‘fiducial parameters’, to match previous results from HR96 for easy comparison. We therefore here focus on the dependence on $r_{\rm p}/a$, but we will in an upcoming paper present further discussions on the dependence on the initial set of orbital angles. An illustration outlining the different outcomes associated with our choice of fiducial parameters is shown in Fig. 2.

II.2 Distant weak encounters

We start by considering the top plot in Fig. 4, which shows $|{\Delta{e}}|=|{e_{\rm fin}}-e_{0}|$ as a function of $r_{\rm p}/a_{0}$. The dashed-dotted lines show the FO result ${\Delta{e}}_{\rm FO}$ from HR96 (Eq. (2)), the dotted lines show the SO result ${\Delta{e}}_{\rm SO}$ from HS19 (Eq. (6)), where the symbols show results from a full three-body integration using the $N$-body code used in (2017ApJ…846…36S, ; 2018ApJ…855..124S, ; 2018MNRAS.481.5436S, ). Note here that only one simulation is performed per point, which naturally gives rise to Poisson scatter.

The numerical and the analytical results clearly approach the simple power-law solution ${\Delta{e}}\propto(r_{\rm p}/a_{0})^{-3/2}$ given by Eq. (2) in the asymptotic limit $r_{\rm p}\gg a$. The slope value $-3/2$ is universal in the sense that it does not depend on the mass ratio or the angular configurations of the three-body system. This universality provides some insight into the role of how important a constant flux of distant encounters might be for the eccentricity evolution. Especially, we might ask the question: does a uniform and infinite sea of weak encounters lead to a divergent result? To gain some insight into this we first notice that the number of single encounters $dN_{\rm s}$ per time interval $dt$, $\Gamma_{\rm s}$, with $r_{\rm p}<R_{\rm p}$ under the assumption of constant number density and velocity dispersion is $\propto R_{\rm p}$ (2018ApJ…853..140S, ). This implies that the rate per peri-center interval $d\Gamma_{\rm s}/dr_{\rm p}\propto k$, where $k$ is a constant. Let us now consider a ‘worst case scenario’ where all single encounters lead to the same positive change in eccentricity ${\Delta{e}}_{+}\propto(r_{\rm p}/a_{0})^{-3/2}$ (note that $\langle{\Delta{e}}\rangle=0$ to FO, but not the variance). Assuming all the eccentricity changes are uncorrelated, one finds in this case that the total change in eccentricity, $d{\Delta{e}}_{\rm tot}$, per unit time, $dt$, from single encounters in the range $r_{\rm p,min}$ to $r_{\rm p,max}$ is given by

[TABLE]

This term does not diverge when setting the upper limit $r_{\rm p,max}=\infty$, which illustrates that the rate of eccentricity change from even an infinite number of distant encounters all leading to a positive change ${\Delta{e}}_{+}\propto(r_{\rm p}/a_{0})^{-3/2}$, does not lead to a diverging result. We therefore expect our results to converge as the upper value of $r_{\rm p}$ in our analysis increases; a value we generally will refer to as $R_{\rm p}$. We note here that an accurate analytical modeling of how the eccentricity changes over time from close and distant encounters can be done using a Boltzmann-approach, which requires the derivation of diffusion coefficients. We will do that in an upcoming paper. In this paper we will be using numerical MC techniques, as described later in Section III.

Finally, the discussions from this section are based on a Newtonian treatment of the problem; however, when PN terms are included several changes are expected from especially precession, which will make the ICs ill-defined and also modify the simple asymptotic solution ${\Delta{e}}\propto(r_{\rm p}/a_{0})^{-3/2}$. This will be discussed further in Section II.5.

II.3 Orbital flip and prompt GW Inspiral

We now consider lower values of $r_{\rm p}/a_{0}$ where the FO power law solution ${\Delta{e}}\propto(r_{\rm p}/a_{0})^{-3/2}$ from Eq. (2) starts to break down. Considering our examples, we see in the top plot of Fig. 4 that for BBHs with an initial eccentricity $e_{0}=0.95$, $e_{0}=0.99$, and $e_{0}=0.999$ that the FO estimates start to deviate from the numerically estimated (correct) solutions for $r_{\rm p}/a_{0}\lesssim 5$, $r_{\rm p}/a_{0}\lesssim 7$, and $r_{\rm p}/a_{0}\lesssim 20$, respectively. Generally, the value for $r_{\rm p}/a_{0}$ below which the FO result starts to break down, a value we denote $(r_{\rm p}/a_{0})_{\rm FO}^{\rm BR}$, is where the change in eccentricity ${\Delta{e}}$ is comparable to $1-e_{0}$, i.e. where the final eccentricity ${e_{\rm fin}}\approx 1$. The physical boundary $e=1$ for bound orbits is not built into the the FO solution, therefore the decrease seen in the figure for the numerically estimated value of ${\Delta{e}}$ compared to the FO power law prediction for $(r_{\rm p}/a_{0})\lesssim(r_{\rm p}/a_{0})_{\rm FO}^{\rm BR}$ originates to simply ‘prevent’ the BBH for reaching values of ${e_{\rm fin}}>1$. This statement is further supported by considering the middle plot of Fig. 4, which shows ${e_{\rm fin}}$ as a function of $r_{\rm p}/a_{0}$.

An approximate value for $(r_{\rm p}/a_{0})_{\rm FO}^{\rm BR}$ can be found by considering the relation ${\Delta{e}}=1-e_{0}={\Delta{e}}_{\rm FO}$ and then solve for $r_{\rm p}/a_{0}$. In the parabolic equal mass limit one finds,

[TABLE]

For our ‘fiducial orbital parameters’ listed in Section II.1 and for $e_{0}=0.99$, one finds using the above relation $(r_{\rm p}/a_{0})_{\rm FO}^{\rm BR}\approx 8$, which indeed provides a reasonable estimate for when FO theory starts to break down. In the limit of very high initial eccentricity one finds that $(r_{\rm p}/a_{0})_{\rm FO}^{\rm BR}\propto(1-e_{0})^{-1/3}$, which shows that for $e_{0}=1-10^{-x}$ the break value $(r_{\rm p}/a_{0})_{\rm FO}^{\rm BR}\propto 10^{x/3}$. One decade closer to $e=1$ therefore moves the break value up by about a factor of $2$. Finally, one should note that for decreasing values of $e_{0}$ the value of $(r_{\rm p}/a_{0})_{\rm FO}^{\rm BR}$ correspondingly approaches $1$ in which case the secular formalism breaks down. The above relation is therefore only useful for BBHs with an initial high eccentricity.

We now focus on our SO expression ${\Delta{e}}_{\rm SO}$ given by Eq. (6), which is shown with dotted lines in the top plot of Fig. 4. As seen here and further discussed in HS19, this relation accurately captures the correct physical form for ${\Delta{e}}$ and thereby resolves the problems associated with the FO result near and below $(r_{\rm p}/a_{0})_{\rm FO}^{\rm BR}$. The decrease in ${\Delta{e}}$ given by the SO relation that ensures that ${e_{\rm fin}}<1$ originates from the second term $\propto\epsilon^{2}$ in Eq. (6), which in this case has a sign opposite to the FO term $\propto\epsilon$. As shown in HS19, the location of $r_{\rm p}/a_{0}$ for which ${\Delta{e}}$ is maximized, i.e. the location of the plateau seen in the plot of ${\Delta{e}}$ slightly below $(r_{\rm p}/a_{0})_{\rm FO}^{\rm BR}$, can be analytically calculated. As in HS19, we refer to the $(r_{\rm p}/a_{0})$ location of the plateau by $(r_{\rm p}/a_{0})_{\rm plateau}$. The solution for $(r_{\rm p}/a_{0})_{\rm plateau}$ follows from solving the equation $\partial{\Delta{e}}_{\rm SO}/\partial(r_{\rm p}/a_{0})=0$, where ${\Delta{e}}_{\rm SO}$ is given by Eq. (6). The solution is long, but can easily be written out in closed form, which provides a clear expression for when SO terms are important (See HS19 for further details).

The dynamics near the plateau are highly interesting, and we now explore them in greater detail. First, we find that as one passes the plateau from high to low values of $r_{\rm p}/a_{0}$, the component of the BBH orbital angular momentum vector projected along the initial ‘z-axis’ (perpendicular to the orbital plane of the BBH), denoted here by $J_{\rm BHz}$, changes sign; a behavior we denote orbital flipping. Loosely speaking, the vector component $J_{\rm BHz}$ flips around, from $\text{sign}(J_{\rm BHz})=+1$ to $\text{sign}(J_{\rm BHz})=-1$, as the BBH ‘reflects off’ the $e=1$ boundary near $(r_{\rm p}/a_{0})_{\rm plateau}$. The final BBH therefore spins in the opposite direction compared to its orbital spin vector prior to the interaction. The sign change is also shown in the upper plot of Fig. 4, where the triangles pointing upwards indicate no sign change ($\text{sign}(J_{\rm BHz})=+1$ ) and the triangles pointing downwards indicate a change ($\text{sign}(J_{\rm BHz})=-1$). As seen, the change indeed happens right near $(r_{\rm p}/a_{0})_{\rm plateau}$ where ${e_{\rm fin}}\approx 1$. We note here that a similar orbital-flip effect has been observed in the Lidov-Kozai case (e.g. 2014ApJ…785..116L, ; 2018CeMDA.130….4S, ). To get a better view of the dynamics, we now consider the top plot in Fig. 5, which shows how the flip occurs using a full three-body integration of the problem. The solid lines show the orbits of the two BHs in their COM, where blue lines highlight their evolution before the orbital flip ($\text{sign}(J_{\rm BHz})=+1$) and the orange lines after the flip ($\text{sign}(J_{\rm BHz})=-1$). As seen, the flip occurs on a timescale that is less than the initial orbital time of the BBH, denoted by $T_{0}$, and therefore represents a ‘non–adiabatic maneuver’. This leads us to another important point, namely that even though the BBH changes orbital spin direction, it does not generally lead to a prompt merger during the flip, although the BBH angular momentum is able to vanish right in the flip (depending on the orbital configurations). In other words, the so-called osculating elements of the BBH $a,e,r_{\rm BH}=a(1-e)$ might predict that $r_{\rm BH}\approx 0$ during the flip which would lead to a prompt merger; however, as the flip does not happen adiabatically, but instead almost promptly in less than one orbit, the BBH is most likely never to pass through the peri-center passage $r_{\rm BH}\approx 0$ as otherwise suggested by the osculating elements. This point is illustrated in the bottom plot of Fig. 5, which shows the actual distance between the two BHs in the binary, $|\bf{r}_{1}-\bf{r}_{2}|$, and the osculating peri-center distance, $r_{\rm BH}=a(1-e)$, as a function of time. As seen, right at the flip near $t/T_{0}\approx 0$, the osculating value $r_{\rm BH}=a(1-e)$ suggests the two BHs (or stars) will undergo an extremely close passage that might lead to a merger; however, the real distance $|\bf{r}_{1}-\bf{r}_{2}|$ shows that they (in this example) never get any closer than their initial peri-center distance. Therefore, one has to be very careful when using orbital average quantities, such as osculating elements, for estimating possible mergers during an interaction.

II.4 Sign Reversal

The last classical effect we discuss is how a BBH through a weak secular encounter can end up with an eccentricity change ${\Delta{e}}$ that has an opposite sign to what the FO estimate predicts; a feature we denote as ‘sign reversal’. An example of this is shown in Fig. 5. As seen in the middle plot showing the BBH eccentricity as a function of $t/T_{0}$, the BBH’s $e$ first increases from $e_{0}=0.99$ to about $e\approx 0.995$, after which the orbit flips (indicated by blue to orange color) and the eccentricity rapidly decreases and settles at an ${e_{\rm fin}}\approx 0.9775<e_{0}$. The FO estimate completely breaks down in this limit and predicts a positive ${\Delta{e}}$ of $\approx 0.046$, where the true change is in fact negative, ${\Delta{e}}\approx-0.012$. In this particular example, where the BBH both undergoes an orbital flip and ‘exceeds’ the $e=1$ limit according to FO theory, the BBH actually ends up with an ${e_{\rm fin}}<e_{0}$ implying its GW inspiral life time $\tau$ increases as a result of the encounter; something that is far from clear when considering FO theory only. In comparison, our SO estimate from Eq. (6) predicts a change of ${\Delta{e}}\approx-0.005$, which correctly is negative, but slightly off as the initial value of $r_{\rm p}/a_{0}=3.0$ is close to where the secular approximation breaks down. Indeed, we find in this case that the exact ${\Delta{e}}$ depends on the initial angular phase of the BBH, which is a clear indication of a secular breakdown.

How and for what $r_{\rm p}/a_{0}$ a BBH undergoes a ${\Delta{e}}$ sign reversal can e.g. be seen in Fig. 4 in the top plot, where the full (empty) symbols indicate $\text{sign}({\Delta{e}})=+1$ ($\text{sign}({\Delta{e}})=-1$)). The locations of sign reversal are also clearly visible in the middle plot, where the crossing between a solid line (${e_{\rm fin}}$) and its corresponding dashed line ($e_{0}$) indicates the location where ${\Delta{e}}=0$. For example, a BBH with an $e_{0}=0.99$ that undergoes a weak encounter with $r_{\rm p}/a_{0}\lesssim 3.5$ will end up with a final eccentricity ${e_{\rm fin}}<e_{0}$, and vice versa. As explained, FO theory is not able to resolve any features related to neither the observed ‘plateau’ nor the ‘${\Delta{e}}$ sign reversal’; however, our SO theory actually predicts an elegant and simple relation between the location of the plateau, $(r_{\rm p}/a_{0})_{\rm plateau}$, and the location of sign reversal, $(r_{\rm p}/a_{0})_{{\Delta{e}}=0}$. As shown in HS19 the two are related by,

[TABLE]

that is, they are universally connected by the simple factor $2^{2/3}\approx 1.59$. Finally, if we consider the functional form of $|{\Delta{e}}|$ for $(r_{\rm p}/a_{0})<(r_{\rm p}/a_{0})_{{\Delta{e}}=0}$ we see that it has a different power law slope than the classical FO prediction $-3/2$. As described in HS19, the slope is different as this is the region where the SO term from Eq. (6) dominates the change in eccentricity. This explains why the slope is observed to be $-3$ as this is exactly the slope associated with $\epsilon^{2}\propto(r_{\rm p}/a_{0})^{-3}$.

II.5 Post-Newtonian Effects

We end this section by discussing the effects from including PN corrects in the weak interaction itself. PN corrections have been shown to be very important to include for modeling mergers in the strong interaction few-body problem (e.g. 2014ApJ…784…71S, ; 2018MNRAS.tmp.2223S, ), and have likewise been shown to play a crucial role in the hierarchical Lidov-Kozai problem (e.g. 2013ApJ…773..187N, ). The question is how important PN corrections are for our weak interaction scenario. In general, General Relativistic (GR) effects in the PN formalism can be divided into two different groups: conservative and dissipative corrections. We discuss both of these in the following subsections.

II.5.1 Conservative Terms: Precession and Shielding

The conservative corrections, which are represented by the 1,2 PN terms, lead to orbital precession. In our case, including precession in the formalism gives rise to at least two notable effects. The first effect is that the initial set of orbital angles (see e.g. Fig. 1), and especially $\Omega$, are no longer staying constant as the perturber moves in from infinity due to precession, and a unique outcome is therefore no longer given by only the classical set of ICs. For a well-defined and unique outcome when PN precession is included, one has to specify exactly where the incoming perturber is for a given set of orbital angles. This also implies that integrating from time $t\in\pm\infty$ for analytically deriving ${\Delta{e}}$ is not possible as one otherwise does in the Newtonian case (HS19). For evolving BBHs in clusters this issue of ‘ill-defined’ classical ICs is not a real problem, as the change in angles from PN precession will be randomized anyway for isotropic random interactions. The second effect appears in configurations where the time it takes for the BBH to precess one full $2\pi$ orbit due to PN apsidal precession (e.g. 1972gcpa.book…..W, ),

[TABLE]

where ${\mathscr{R}_{\rm m}}$ denotes the Schwarzschild radius of a BH with mass $m$, is shorter than the time scale it takes for the encounter to pass the BBH, a time we here refer to as the peri-center passage time $T_{3}$,

[TABLE]

where we have omitted terms of order unity. In such cases, the interaction is no longer simply described by an incoming single interacting with a BBH on a closed Keplarian orbit, as the one illustrated in Fig. 1, but instead by a single interacting with a ‘disc’ build up by ‘rosetta’ orbits with inner and outer annuli equal to the BBH’s peri-center and apo-center, respectively. This symmetry change from a non-axisymmetric closed orbit to an axisymmetric disc is expected to ‘shield’ the BBH from perturbations, in the same way as a circular BBH is much harder to perturb than an eccentric (e.g. 1996MNRAS.282.1064H, ). Using the above equations (10) and (11), one finds that such shielding will start to appear for encounters satisfying the following inequality,

[TABLE]

where we again have omitted factors of order unity. As seen, even for BBHs with initial high eccentricity, corrections from PN precession are likely not to play a significant role. However, one might still find resonant couplings between the precession and encounter frequencies for lower values of $r_{\rm p}/a$, which could lead to PN precession corrections. We do not have the necessary analytical tools to quantify these highly interesting effects yet, but we will follow up on this in upcoming work. Also, for this work we updated our PN code to also include the 1,2PN terms, but as also noticed in (2018PhRvD..98l3005R, ), it is numerically challenging to evolve highly relativistic and eccentric BBHs with these corrections.

II.5.2 Dissipative Terms: Prompt GW Inspiral

The dissipative corrections, which first appear at the 2.5PN order, lead to orbital decay for a BBH, which will result in a merger if its eccentricity is high enough for its GW inspiral life time to be shorter than its dynamical disruption time. As described in Section II.3, a BBH can be driven to a very high eccentricity $e$ by weak encounters near the plateau, and can in principle reach $e\approx 1$ during orbital flips. The question is what happens to such interactions when dissipative corrections, here included by the 2.5 PN term, are included. An example is shown in Fig. 6, which shows the orbital evolution of a BBH over the duration of a weak encounter for an interaction near the plateau. As seen, the weak encounter drives the BBH to such a high eccentricity that the 2.5PN correction term starts to take over the evolution, which leads to a clear orbital decay during the interaction. The inspiral time of such binaries are of order the peri-center passage time of the encounter w.r.t. the COM of the BBH, which interestingly fills in the ‘gap’ in GW peak frequency space between the 3-body mergers and 2-body mergers studied in (2018MNRAS.tmp.2223S, ) and later discussed in Section III.2.3. Note here that Fig. 6 shows results for a BBH that is similar to the one shown in Fig. 5, the only difference is that the one considered in Fig. 6 undergoes a weak encounter from further away; GR effects are therefore not generally associated with only close strong and weak interactions, but can play a significant role for even very distant perturbers. Finally, in Fig. 4 results from including the 2.5PN term in weak interactions are shown with black symbols, where no symbol for the upper two figures, and a cross symbol in the bottom plot, indicates the BBH merged during the weak interaction. As seen, for interactions where the BBH do not merge during the interaction the final outcomes are similar to the ones derived using Newtonian dynamics alone. For interactions where the BBH does merge during the interaction, we find notable differences when including the 2.5PN term in especially the final inspiral time, as seen in the bottom plot. As pointed out in the above paragraph, such mergers can lead to distinguishable peaks in GW peak frequency space, and are therefore in general important to keep track of.

II.5.3 Summary

To conclude, including the 2.5PN dissipative correction is easy and gives rise to potentially important and interesting prompt GW inspiral mergers forming during or right after the weak interaction. The conservative 1,2PN corrections are more difficult to include in the problem, as they also often lead to numerical challenges in the relevant regimes where the 2.5PN correction is also important. The 1,2PN terms can lead to shielding effects that protect the binary from changing eccentricity for very distant encounters; however, using simple scaling arguments, we have proven that this correction is likely not to play a significant role for astrophysical cluster BBHs, but we will study this in further detail in an upcoming paper. For these reasons, all our studies later on in this paper are based on including the 2.5PN term alone.

III Black Hole Binaries in Clusters

Having addressed in the above Section II the outcome of a single weak interaction, we now turn to the question of what the outcome is of a BBH undergoing a series of weak and strong binary-single interactions inside a dense cluster.

BBHs in dense clusters can be driven to merger during their dynamical hardening through few-body interactions such as binary-single (2014ApJ…784…71S, ) and binary-binary (2019ApJ…871…91Z, ) interactions. Following (2018MNRAS.tmp.2223S, ), and later described in Section III.1.2, we refer to a BBH that merges during a strong binary-single interaction by a 3-body merger, and a BBH that merges in-between its strong interactions by a 2-body merger. As further described in (2018MNRAS.tmp.2223S, ) and (2019arXiv190102889S, ), such mergers can only form if their GW inspiral time is less or comparable to the characteristic dynamical time of their associated interaction channel, which for 3-body mergers is $\sim 1$ year (characteristic orbital time of a BBH in a typical GC), and for 2-body mergers is $10^{7}$ years (characteristic time for a BBH to undergo its next strong encounter in a typical GC). These characteristic time scales directly maps to how the merging BBHs distribute in GW peak frequency $f_{\rm GW}$ and orbital eccentricity $e$. Following (2018PhRvD..97j3014S, ), the value of $f_{\rm GW}$ for a BBH with peri-center distance $r_{\rm BH}$ can be approximated by the GW frequency found from assuming the BBH is circular with a SMA equal to $r_{\rm BH}$. From using that the emitted GW frequency is two times the Keplerian orbital frequency one now finds,

[TABLE]

In (2018MNRAS.481.4775D, ) it was illustrated that a BBH that has to inspiral on a time scale $\mathcal{T}$, will have an $f_{\rm GW}$ equal to

[TABLE]

From this relation one concludes that 2-body mergers with $\mathcal{T}\sim 10^{7}$ years will form near the LISA band, where 3-body mergers with $\mathcal{T}\sim 1$ year form near the LIGO band. As pointed out in (2018MNRAS.tmp.2223S, ) and (2018MNRAS.481.4775D, ), these dynamically formed BBH populations will therefore leave unique observational imprints across the LISA, DECIGO and LIGO bands, which can be used to probe their dynamical history and astrophysical origins. This is highly interesting and encouraging; however, what has not been discussed in these recent studies is the effect from weak interactions, and their possible importance for an accurate modeling of the observable $f_{\rm GW}$ and $e$ distributions. We address this question in this section.

III.1 Post-Newtonian Cluster Model

In this section we start by introducing our cluster model, and how we keep track of different GW merger outcomes. We then describe in simple words the routines of our MC cluster algorithm, and especially how we include the effects from both GW emission, and strong and weak interactions111Our developed code is free to use, and can be obtained by contacting the lead author of the paper.. In the last part we present our main results on the role of weak interactions using our introduced MC cluster algorithm.

III.1.1 Strong Interactions and Dynamical Ejection

For modeling BBH mergers forming in clusters including weak and strong interactions, we make use of the following simple model inspired by (2018MNRAS.tmp.2223S, ). We consider a BBH that starts out with an initial SMA $a$ that is slightly below the hard-binary (HB) value given by

[TABLE]

We imagine this BBH lives in a population of single BHs, all with the same mass $m$, characterized by a constant number density $n_{\rm s}$, velocity dispersion $v_{\rm dis}$, and escape velocity $v_{\rm esc}$. The BBH will undergo both strong and weak interactions defined in the following way: We denote an interaction a strong interaction if the incoming single BH passes the BBH on an orbit with peri-center distance $r_{\rm p}<\mathscr{C}a$, where our fiducial value of $\mathscr{C}$ is here $\mathscr{C}=2$, and $a$ is the SMA of the BBH at the time the single passes the BBH. An interaction with $r_{\rm p}>\mathscr{C}a$ will instead be labeled a weak interaction. We assume that each strong interaction leads to a constant decrease in the SMA of the BBH from $a$ to $a_{\rm enc}$ as,

[TABLE]

where we use $\delta=7/9$, which equals the average value assuming the chaotic binary-single interaction can be modeled by an ergodic distribution as discussed in (2018PhRvD..97j3014S, ). At the same time we assign the BBH a new eccentricity $e$ that is drawn from a so-called thermal distribution $P(e)=2e$ (Heggie:1975uy, ). We choose this simple model for strong interactions due to the following three reasons. $(i)$ It provides full analytical solutions to both the stationary BBH population (e.g. 2019PhRvD..99f3006S, ) and the highly eccentric evolving population across all GW bands (e.g. 2018PhRvD..97j3014S, ; 2019arXiv190102889S, ). This allows us to not only cross check all our algorithms, but also to extract small changes coming from additional terms such as weak interactions. $(ii)$ The simple $\delta$-model has been used in all our previous studies, which allows for an easy comparison. Despite its simplicity, our previous results (e.g. 2018MNRAS.tmp.2223S, ) are in surprisingly good agreement with more sophisticated MC techniques (e.g. 2019PhRvD..99f3003K, ). $(iii)$ Adding weak interactions into the evolution of BBHs requires in principle large computational resources, as weak interactions need to be performed $\gg 1$ times in-between each strong interaction. Therefore, partly due to computational limitations, we restrict ourself to our simple $\delta$-model combined with additional hybrid techniques, as described later.

This process of a gradual decrease of the BBH’s SMA is loosely referred to as hardening. In the point-particle Newtonian limit this will always end with a dynamical ejection of the BBH when its SMA reaches the following critical value,

[TABLE]

which is where the recoil velocity from the 3-body interaction equals the escape velocity of the cluster (2018PhRvD..97j3014S, ). This classical hardening process therefore leads to an ejection of BBHs with a SMA $\sim a_{\rm ej}$ and an eccentricity distribution that is thermally distributed. As shown by (2000ApJ…528L..17P, ; 2016PhRvD..93h4029R, ), this population of ejected BBHs can explain at least part of the observed LIGO merger rate, but as pointed out by e.g. (2018PhRvD..97j3014S, ; 2018MNRAS.tmp.2223S, ), this picture is highly incomplete as it only includes (energy conserving) Newtonian dynamics. This is described in greater detail below.

III.1.2 Energy Dissipation and Few-body Mergers

When GR dissipative terms – where the leading term is the 2.5PN term in the PN expansion formalism – are included in the dynamics, a BBH can also merge inside its cluster either during or in-between its strong interactions. We describe our methods for modeling these outcomes in the following paragraphs.

3-body Mergers: For determining if a BBH merges during a strong 3-body interaction, we model this often highly chaotic and resonating state using the approach put forward in (2018PhRvD..97j3014S, ). The interaction here is divided up into $N_{\rm IMS}$ intermediate states (IMSs), each of which is described by a BBH with a bound single BH. For each IMS we assign the corresponding BBH an eccentricity sampled from the thermal distribution $P(e)=2e$, but keep the SMA fixed to its initial value $a_{0}$. To determine if the BBH undergoes a GW inspiral merger during this IMS, i.e. merge while the third objects is still bound to it, we compare the energy radiated over one orbit of the BBH,

[TABLE]

to the total orbital energy of the bound 3-body state,

[TABLE]

If $\Delta{E}_{\rm GW}>E_{\rm B}$ we label the IMS assembled BBH as a 3-body merger. This energy threshold is equivalent of saying that the BBH will undergo a GW inspiral merger during the interaction if its peri-center distance $r_{\rm BH}<r_{\rm 3cap}$, where $r_{\rm 3cap}$ is a ‘3-body capture distance’ given by (2018PhRvD..97j3014S, ),

[TABLE]

If instead $r_{\rm BH}>r_{\rm 3cap}$ the BBH does not merge during the considered IMS. We repeat this process, i.e. first assign the IMS assembled BBH an eccentricity from $P(e)=2e$ and then compare its $r_{\rm BH}$ to $r_{\rm 3cap}$, up to $N_{\rm IMS}=20$ times per interaction, which is the average number of IMSs assembled in a chaotic equal mass 3-body interaction (2018PhRvD..97j3014S, ). The total probability for a BBH to merge during a 3-body interaction, $P_{\rm 3m}$, i.e. to be formed with orbital parameters such that $r_{\rm BH}<r_{\rm 3cap}$, can in this model be expressed analytically as follows (2018PhRvD..97j3014S, ),

[TABLE]

As illustrated in (2018PhRvD..97j3014S, ), by integrating this probability up over a BBHs hardening binary-single interactions from $a_{\rm HB}$ to $a_{\rm ej}$, one finds that the total probability for a BBH to merge during a 3-body interaction relative to merging outside the cluster within a Hubble time is $\approx 10\%$.

2-body Mergers: For keeping track of the formation of 2-body mergers, we evolve the orbital parameters $a,e$ of the BBH in question in-between each interaction according to the equations given by (Peters:1964bc, ). If at any given time the GW inspiral life time of the BBH right after it has finished an interaction, $\tau$, is less than the time it takes before it undergoes its next interaction, $t_{\rm int}$, we stop the code and label the outcome a 2-body merger. To calculate the time before the next interaction, we use the same approach as for Eq. (4). Assuming that we keep track of single encounters up to a maximum peri-center distance of $R_{\rm p}$, the time in-between interactions can be approximated by,

[TABLE]

As argued in Section II.2, one does expect the results to converge for $R_{\rm p}$ going to infinity. We will address this issue later in Section III.2.3.

III.1.3 BBH Evolution with Weak Interactions

We now describe our numerical scheme used in this paper for modeling the evolution of BBHs inside a dense stellar cluster with the inclusion of weak interactions. It expands on a similar model put forward by (2018PhRvD..97j3014S, ; 2018MNRAS.tmp.2223S, ; 2019arXiv190102889S, ).

We start by first choosing the maximum relative peri-center distance within which we model the effect from interactions, $R_{\rm p}/a$, where $a$ is the SMA of the BBH in question. This rescaled value stays constant, but the physical value of $R_{\rm p}$ will decrease as the BBH hardens inside its cluster. If only strong interactions are included, then $R_{\rm p}/a\sim 1$. Therefore, a value of $R_{\rm p}/a\gg 1$ will allow the orbital parameters of the BBH to be affected by a much wider and numerous range of interactions compared to the previously studied strong interaction case. With this choice of $R_{\rm p}/a$, we follow the evolution of a BBH from its initial SMA, $a_{\rm in}\sim a_{\rm HB}$, towards its final value, $a_{\rm fin}\sim a_{\rm ej}$, assuming all properties of the cluster stay constant. As $R_{\rm p}/a$ takes a finite value we are able to evolve the BBH in finite steps. In the following we describe how we evolve the BBH from step $i$ to its next step $j$. See also Fig. 7 which shows a simplified flowchart of our algorithm.

(1) Check for 2-body merger: Based on the BBH’s orbital values at step $i$, referred to as $\{a_{i},e_{i}\}$, and the value of $R_{\rm p}/a$, we first derive $\tau$ and $t_{\rm int}$ using Eq. (3) and (22), respectively. If $\tau<t_{\rm int}$ we stop the code and record the BBH as a 2-body merger. If instead $\tau>t_{\rm int}$, we proceed as follows.

(2) GW orbital evolution: If no 2-body merger has occurred, we evolve the BBH’s orbital parameters using the equations from (Peters:1964bc, ) over the time interval it takes for it to undergo its next encounter $t_{\rm int}$. During this evolution the orbital parameters of the BBH evolve from $\{a_{i},e_{i}\}$ to $\{a_{\rm int},e_{\rm int}\}$.

(3) Pick encounter distance: We now prepare for the encounter. For this we start by randomly picking a peri-center distance for the encounter w.r.t. the COM of BBH, $r_{\rm p}/a_{\rm int}$, in the range $[0:R_{\rm p}/a_{\rm int}]$ assuming the encounters are isotropically distributed at infinity (see e.g. (2014ApJ…784…71S, )). At this step we just record the value $r_{\rm p}/a_{\rm int}$.

(4) Classify and perform encounter: At this step we start by classifying the encounter type from the above step. We do this as follows.

- Strong Interaction: If $r_{\rm p}/a_{\rm int}<\mathscr{C}$, we model the interaction as a strong interaction. To evaluate the corresponding outcome we first estimate if the interaction results in a $3$-body merger using the IMS-decomposition approach described in Section III.1.2. If a $3$-body merger does not form, we model the interaction outcome using our simple $\delta$-model described in Section III.1.1. If at this point $a_{\rm int}<a_{\rm ej}$, we stop the code after the interaction and record the BBH as being ejected from the cluster.

- Weak Interaction: If instead $r_{\rm p}/a_{\rm int}>\mathscr{C}$, we model the interaction as a weak interaction. For this we start by picking a random configuration for the angles $\{i,\omega,\Omega\}$ assuming an isotropic distribution. We then derive the outcome using one of the following ways. $(i)$: Use our few-body code described in (2017ApJ…846…36S, ). $(ii)$: Use our SO result from Eq. (6). $(iii)$: Use a ‘hybrid approach’ where we combine $(i)$ and $(ii)$. As discussed later, option $(i)$ is numerically challenging, but straight forward. Option $(ii)$ is in-accurate for $r_{\rm p}/a_{\rm int}\sim\mathscr{C}$, but extremely fast. Therefore, option $(iii)$ is often the best option, but it requires bridging $(i)$ and $(ii)$.

- Final parameters: If the BBH has not merged at this point its final set of orbital values will be $\{a_{j},e_{j}\}=\{a_{\rm int}+\Delta{a},e_{\rm int}+\Delta{e}\}$, with $\Delta{a}$ and $\Delta{e}$ derived as just described.

(5) Next step and termination: For a BBH that has survived all interactions and is still inside the cluster, we end this ‘one interaction’-routine from step $i$ to $j$ by simply putting $\{a_{\rm j},e_{\rm j}\}\rightarrow\{a_{\rm i},e_{\rm i}\}$. After this we go back to the start and repeat the process. We keep doing this until the BBH either merges inside the cluster or is ejected.

Using this algorithm we can quickly build up a distribution of BBH properties showing e.g. the distribution of BBH eccentricity $e$ across different $f_{\rm GW}$ bands from LISA to LIGO. Our model does represent a simplified picture of how BBHs undergo dynamical interactions in dense stellar systems; however, our extremely fast approach has shown to give surprisingly accurate predictions for how BBHs merge and distribute across both $f_{\rm GW}$ and $e$. Our model can easily be improved, but its simplicity as well as its corresponding list of available analytical solutions makes it an ideal testbed for including new effects, as the ones we consider in this paper. Below we study outcomes from this model.

III.2 Results from Including Weak Encounters

In this section we present our main results for how BBHs interact and merge in clusters under the influence of both weak and strong interactions. After defining our fiducial parameters and approach below, we study the dynamical history of a single BBH. After this we move to results derived from a large sample of evolved BBHs using MC techniques, from which we derive $f_{\rm GW}$ and $e$ distributions relevant for especially LISA. As illustrated, we find clear observational differences when weak interactions are included.

III.2.1 Fiducial Parameters

In all our models we use the algorithm described in the above Section III.1.3. If not stated otherwise, we further assume the following set of ‘fiducial model’ parameters: $m=20M_{\odot}$, $n_{\rm s}=10^{5}\ \text{pc}^{-3}$, $v_{\rm dis}=10\ \text{kms}^{-1}$, $v_{\rm esc}=50\ \text{kms}^{-1}$, $R_{\rm p}/a=50$, and $\mathscr{C}=2$. These are chosen for easy comparison to previous work. For our main results presented below, we also employ the following ‘hybrid approach’ for accurate and fast modeling of the outcome from weak encounters: for encounters with $2<r_{\rm p}/a<\max(5,(r_{\rm p}/a_{0})^{\rm BR}_{\rm FO})$ we use a full 2.5PN three-body integration to evaluate the outcome, where for $r_{\rm p}/a>\max(5,(r_{\rm p}/a_{0})^{\rm BR}_{\rm FO})$ we use HS19 given by Eq. (6). In Section III.3 we compare results from this hybrid approach to results using HS19 only.

III.2.2 A Single Black Hole Binary

We start by studying the evolution of a single BBH undergoing both weak and strong encounters inside a dense cluster. Results are presented in Fig. 8, which shows in the top, middle, and bottom plots the BBH’s eccentricity $e$, the value of $r_{\rm p}/a_{0}$, and the BBH’s peak GW frequency $f_{\rm GW}$, respectively, as a function of the number of scatterings, $N_{\rm sc}$ (scattering counter). For this modeling we have used our numerical scheme described in Section III.1.3, with the hybrid approach and fiducial cluster parameters listed in the above Section III.2.1.

As seen, the BBH occasionally undergoes notable evolution in eccentricity in-between its strong interactions (red crosses). For example, the BBH evolves from an eccentricity of $\sim 0.8$ to $\sim 0.9$ over the range $N_{\rm sc}\approx 400-500$ entirely through weak interactions (grey dots). A relevant question is here if weak interactions therefore are expected to lead to an observed eccentricity distribution that is notably different from the thermal distribution that often follows from strong scatterings. However, by studying a large number of runs of BBHs evolving through both weak and strong interactions, we did not find any significant differences when weak interactions are included. We will however study this in more detail in a upcoming paper. See also (2019ApJ…872..165G, ) for a similar analysis.

The importance of including weak encounters becomes clear when considering the end of the evolution of the BBH shown in Fig. 8. As seen, the last strong interaction leaves the BBH with an eccentricity of $e\approx 0.9$, which is not even high enough for it to merge before its next strong interaction (if the BBH has an inspiral time $\tau<t_{\rm SE}$, where $t_{\rm SE}$ is the time interval before the next strong encounter (see Eq. (4)), the state is marked with a grey diamond). However, the BBH is instead driven to merge entirely through weak interactions, each of which changes the BBH’s eccentricity until $\tau<t_{\rm int}$. This evolution also brings the BBH into the LISA band, as seen in the bottom plot, which will make it appear as a high eccentricity LISA source. Such eccentric sources are not a natural outcome of isolated BBH evolution, and therefore provide a test of the cluster channel (2018MNRAS.tmp.2223S, ; 2019PhRvD..99f3003K, ).

Finally, the top plot shows how the analytical estimates HR96 (green circles), and HS19 (blue circles), compare to our full numerical results (black dots). As seen, both HR96 and HS19 provide excellent estimates when $r_{\rm p}/a\gg 1$; however, both start to break down, as expected, when approaching the strong interaction limit. The important difference between HR96 and HS19 is seen for BBHs that scatter around the $e\approx 1$ limit. For example, HR96 clearly gives rise to several unphysical outcomes with $e>1$ when the BBH is already highly eccentric; a problem that is resolved by using HS19 due to its corrections from the SO term. The high eccentricity region is the most interesting and relevant for understanding and modeling how BBHs are driven to merge, which again illustrates the importance of using corrections similar to the one encoded in HS19.

III.2.3 A Population of Black Hole Binaries

We now turn to the study of what GW signals to expect from a population of BBHs evolving through weak and strong interactions, similar to the population that we potentially observe from GCs throughout the nearby observable universe. For this, we here model the dynamics of an ensemble of $15000$ individual BBHs, and follow their evolution in a similar way as the one described in the above Section III.2.2. We record 2-body and 3-body BBH mergers, as well as properties for the ejected BBH population. From this data we derive $f_{\rm GW}$ and $e$ distributions, as described in the following sections. In this sub-section all our results are based on the fiducial values and hybrid approach outlined in Section III.2.1.

We first consider Fig. 9, which shows the distribution of $f_{\rm GW}$ (top plot) and the eccentricity distribution at $10^{-2}$ Hz (bottom plot) of the merging BBH population. As described in the caption, $f_{\rm GW}$ is the GW peak frequency recorded at the time of ‘formation’, i.e. right when either the BBH in question is ejected from the cluster for ejected mergers, $\tau<t_{\rm int}$ for 2-body mergers, or right when $r_{\rm BH}<r_{\rm 3cap}$ for 3-body mergers. This distribution is therefore not universal, as $t_{\rm int}$ depends on $R_{\rm p}$ (a higher $R_{\rm p}$ will allow each BBH to inspiral over a longer time scale before the code halts, which will drive the 2-body merger distribution slightly to the right). However, it only weakly depends on $R_{\rm p}$ as each BBH inspirals with near constant $f_{\rm GW}$, and the plot therefore still provides some useful insight into how the different outcomes distribute and differ from each other. As seen, ejected mergers tend to form below both the LISA and the LIGO bands, and will therefore almost circularize before being observable, whereas 2-body mergers will appear eccentric in LISA, and 3-body mergers in LIGO. As described in (2018MNRAS.481.4775D, ), and illustrated by Eq. (14), this splitting of outcomes across $f_{\rm GW}$ is linked to the characteristic inspiral times of each population. The bottom plot of $e$ at $10^{-2}$ Hz is in principle not dependent on $R_{\rm p}$, as long as the problem converges with increasing $R_{\rm p}$. From this plot its clear that the 2-body mergers dominate the population of eccentric sources in the LISA band. The different colored lines in both plots illustrate different dynamical subsets, as described in the figure caption. Starting with comparing Subset A (last strong interaction has $\tau>t_{\rm SE}$, i.e. BBHs are here driven to merger through one or more weak encounters) and Subset B (last strong interaction has $\tau<t_{\rm SE}$, i.e. the BBHs will here still be labeled mergers using a strong interactions code only), we see that there is no major difference between the two. It will therefore generally be hard to tell the difference between strong and weak interaction driven BBH mergers from a few LISA observations alone (in this statement we have ignored rare exotic mergers such as a BBH promptly driven to merger through an interplay between PN corrections and a dynamical flip-orbit during a weak interaction). The only major difference across $f_{\rm GW}$ from the weak interaction channel is the appearance of Subset C, which originates from BBHs that merge during the weak encounter (see e.g. Fig. 6). As seen, these mergers appear at relatively high frequencies where they make up their own ‘bump’ slightly below the 3-body mergers.

To systematically explore the role of including weak encounters we now consider Fig. 10, which shows the distribution of $f_{\rm GW}$ (top plot) and the eccentricity $e$ distribution at $10^{-2}$ Hz (bottom plot), where the solid black and solid red lines show results from including weak interactions ($R_{\rm p}/a=50$) and not including weak interactions ($R_{\rm p}/a=2$) in the BBH evolution, respectively. The companion Fig. 11 shows the cumulative distribution of the eccentricity distribution at $10^{-2}$ Hz from Fig. 10. From considering these figures, one sees that including weak interactions leads to a notable increase in the number of GW sources that will appear with measurable eccentricity ($>10^{-2}$) in the LISA band. For example, the number of sources with $\log(e)>-1.0$ at $10^{-2}$ Hz increases by about $\approx 50\%$. The reason can be found from considering Fig. 12, which shows the number of in-cluster mergers (2-body + 3-body mergers) relative to the number of BBHs ejected from the cluster with $\tau<t_{\rm H}$, where $t_{\rm H}$ denotes the Hubble time, as a function of $R_{\rm p}/a$. As seen, as more and more weak encounters are included, i.e. as we increase $R_{\rm p}/a$, more of the BBHs will merge inside the cluster before being ejected. The increased number of eccentric LISA sources therefore primarily originates from an increase in the number of in-cluster mergers. The corresponding slight systematic shift towards higher $e$ seen in the bottom plot of Fig. 10 for the in-cluster mergers also help increasing the number.

We do not yet have a solid analytical understanding for how the relative importance of weak interactions change with the cluster parameters as we do for strong interactions (e.g. 2018PhRvD..97j3014S, ); however, the first results presented here are interesting, and weak interactions definitely add an interesting complexity to the problem. They are also highly important to include for an accurate modeling of possible secular effects on the BBHs if the cluster is not perfectly spherical (2019arXiv190201344H, ; 2019arXiv190201345H, ), and might also play a key role in accurately predicting the formation rate of 2. generation mergers (e.g. 2016ApJ…824L..12O, ). We reserve several of these questions for future work.

III.3 Fast Analytical Monte Carlo Approach

We finish this section by exploring how well one can do by using Eq. (6) for all weak interactions, i.e. for all $r_{\rm p}/a_{0}>\mathscr{C}(=2)$, instead of employing our hybrid approach that includes numerical simulations; an approach we refer to as our ‘Analytical MC Approach’. In this approach it is unavoidable to encounter situations where either $e>1$ or $e<0$, also when using the SO term, therefore, every time such a case appears we simply put ${\Delta{e}}=0$ without stopping the code. Results from this approach is shown in Fig. 10 and 12. As seen, despite some Poisson scatter, the Analytical MC Approach seems to perform very well. A great overlap between the hybrid and the analytical MC approach is especially seen in Fig. 10, where both of them clearly distinguish themselves from the runs where only strong interactions are included.

This is highly encouraging as the analytical MC approach is extremely fast and easy to code up, which especially allows one to explore a wide range of cluster parameters in a very short amount of time. This approach is therefore also well suited for performing a more systematic study on the importance of weak interactions as function of the properties of the cluster and the interacting BHs.

Finally, there are several easy ways of improving our proposed analytical MC approach, two of which are: $(1)$ For strong interactions use the full distribution of $\delta$ given e.g. in (Heggie:1975uy, ), instead of the average value $\delta=7/9$. $(2)$ For weak interactions, include the change in SMA $\Delta{a}$ using (2003CeMDA..87..411R, ), instead of simply assuming $\Delta{a}=0$. On top of these, it is possible to improve on the distribution of 3-body mergers, and also include 4-body mergers. In an upcoming paper we will quantify the role of single-single GW capture mergers. All of these studies will become central in developing the next generation of fast MC codes for modeling dense stellar systems and their formation of BBH GW mergers.

IV Conclusions

We have in this paper quantified the effects from weak encounters on the assembly of BBH mergers in GCs using a combination of novel analytical and numerical techniques. Current state-of-the-art codes (e.g. 2018PhRvD..98l3005R, ; 2013MNRAS.431.2184G, ), as well as recent analytical studies (e.g. 2018MNRAS.tmp.2223S, ; 2018MNRAS.481.4775D, ) have only included strong interactions in their modeling; however, BBHs in GCs will in-between each strong interaction undergo hundreds of weak interactions, each of which changes the eccentricity of the BBH slightly and thereby its GW inspiral life time. It is therefore possible for BBHs to both be driven towards merger and also away from merger entirely through weak interactions.

Using a simple cluster model with the inclusion of weak and strong interactions, we find here that the addition of weak interactions leads to a notable and systematic relative increase in the number of BBHs merging in-side the cluster. For our fiducial cluster parameters we find that the total number of observable BBH mergers goes from being dominated by ejected mergers to in-cluster mergers when weak interactions are included. Because of the relative short GW inspiral time of in-cluster mergers, this population will appear as eccentric LISA sources, and therefore add an important piece to the puzzle of how to disentangle BBH merger channels from each other by the use of GW observations.

Our main results are based on a hybrid approach where we combine full 2.5PN numerical scattering experiments with analytical expressions. We have in this regard especially demonstrated that very similar results are achieved if one uses our new SO perturbative solution (2019arXiv190409624H, ) for modeling the changes in BBH eccentricity for all weak interactions. This is highly encouraging, as this allows for a full analytical MC modeling of the problem including both weak and strong interactions as well as 2.5PN effects. This makes it e.g. possible to derive distributions of BBH merger properties in just a few seconds for any cluster parameters. Of course, this approach is based on a simplified cluster model, but previous work based on the same assumptions turns out to be in surprisingly good agreement with more sophisticated MC models (e.g. 2019arXiv190102889S, ).

We have in this paper only considered results from a single set of cluster parameters, but we plan in the near future to explore the relative importance of weak interactions for a range of cluster parameters from less dense to very dense nuclear star clusters. It might be that in some regimes weak interactions play an even larger role than found in this paper; however, we reserve that for future work.

Finally, we note that the weak interactions we have here considered are essential to include for studying the secular evolution of BBHs in e.g. non-spherical clusters (2019arXiv190201344H, ; 2019arXiv190201345H, ). The reason is that weak interactions not only give rise to a change in eccentricity, but also to a change in the orbital orientation. This change can interfere with the secular change induced from a non-spherical potential, and thereby affect the distribution of BBHs driven to merger in-side clusters through this mechanism. Whether or not this is a competing channel for forming eccentric LISA sources will be studied in upcoming work.

Acknowledgments. — We thank Stephen McMillan and Piet Hut for stimulating discussions. J.S. acknowledges support from the Lyman Spitzer Fellowship. A.S.H. gratefully acknowledges support from the Institute for Advanced Study, and the Martin A. and Helen Chooljian Membership.