The Population of Eccentric Binary Black Holes: Implications for mHz Gravitational-wave Experiments

To illustrate the scalings from Section 2 we assume a generic eccentric BBH population motivated by dynamically formed systems in dense stellar environments and few-body systems undergoing LK oscillations. We assume equal mass BBHs with ${m}_{1}={m}_{2}=30\,{M}_{\odot }$, an initial semimajor axis distribution that is log-uniform between 1 and 1000 au, and a thermal eccentricity distribution (Jeans 1919), i.e., e² is uniform in [0, 1]. Only systems with t_gw less than the Hubble time are included, since only these will contribute to the observed merger rate. The thermal eccentricity distribution is motivated by the properties of dynamically formed binaries produced in dense stellar systems (e.g., Samsing & D'Orazio 2018), and by the fact that it produces a uniform distribution of r_p0 in the $e\to 1$ limit, equivalent to a uniform distribution of the angular momentum squared J², which is a natural consequence of non-secular stochastic angular momentum kicks due to the tertiary in hierarchical triple systems (e.g., Katz & Dong 2012).

We exclude systems from our sample with

$$\begin{eqnarray}&&{r}_{p}\lt 6.4\times {10}^{-5}\,\mathrm{au}{\left(\displaystyle \frac{{m}_{1}{m}_{2}a}{{30}^{2}{M}_{\odot }^{2}\mathrm{au}}\right)}^{2/7}{\left(\displaystyle \frac{M}{60{M}_{\odot }}\right)}^{1/7},\end{eqnarray} \tag{ 10 }$$

in the $e\to 1$ limit, because the fractional change in the orbital energy per orbit is of order unity and the secular equations break down (i.e., t_gw ≤ P). Such systems emit GWs at a peak frequency of

$$\begin{eqnarray}&&{f}_{p}\gt 0.95\,\mathrm{Hz}\,{\left(\displaystyle \frac{{30}^{2}{M}_{\odot }^{2}\mathrm{au}}{{m}_{1}{m}_{2}a}\right)}^{3/7}{\left(\displaystyle \frac{M}{60{M}_{\odot }}\right)}^{2/7}.\end{eqnarray} \tag{ 11 }$$

These extreme systems are not of relevance for the main comparison in the LISA band from 0.1 to 10 mHz, but may arise in nature (Silsbee & Tremaine 2017; Samsing & D'Orazio 2018; Kremer et al. 2019).

We evolve any given binary system "i" in the sample using the secular equations and calculate ${\dot{f}}_{p}^{(i)}$ as a function of f_p. We assume each system represents an equilibrium population with the same initial conditions undergoing migration toward coalescence. Thus, at any peak frequency, the population i has a distribution ${{dN}}^{(i)}{/{df}}_{p}\propto 1/{\dot{f}}_{p}^{(i)}$, and the total number of systems has a frequency distribution ${dN}/{df}={\sum }_{i}{{dN}}^{(i)}{/{df}}_{p}\,=({\rm{\Gamma }}/N){\sum }_{i}[1/{\dot{f}}_{p}^{(i)}]$, where N is the sample size, and the distribution is normalized to the LIGO rate of $50\,{\mathrm{Gpc}}^{-3}\,{\mathrm{yr}}^{-1}$. To obtain the number of systems N₁₂ in frequency bin [f₁, f₂], we integrate the frequency distribution over the bin, ${N}_{12}=({\rm{\Gamma }}/N){\sum }_{i}{t}_{12}^{(i)}$, where ${t}_{12}^{(i)}$ is the time spent in the frequency bin for population i.

Figure 2 shows the histograms of systems in the frequency range [0.01, 10] mHz assuming that all the LIGO mergers come from either the eccentric channel or the circular channel. The upper panel shows the number of systems per logarithmic frequency bin, with the eccentric systems broken into sub-samples by orbital period. For the BBH population we consider here, when a > 11.5 au, the orbital period exceeds 5 yr, which is roughly the operation timescale of LISA. Like single-transit planet detections in transit surveys (e.g., Villanueva et al. 2018), only a single pulse may be seen during the mission. However, systems with orbital periods of days or months will provide many repeated pulses during the entire mission (Section 4.1).

Zoom In Zoom Out Reset image size

The bottom panel of Figure 2 shows the ratio of the number of eccentric systems to the number of circular systems in bins of frequency. Note that this ratio does not depend on the overall normalization of the LIGO rate. We find ≃6–10 times more BBHs from the eccentric channel than that predicted by the circular channel at around 0.1 mHz, decreasing to ≃2 times more at ≃3 mHz. In absolute numbers, we find that ≃45, 90, and 116 eccentric BBHs in our Galaxy are currently emitting in the 0.1–1 mHz range with orbital periods P ≤ 1, 10, and 100 days, respectively. These numbers are ≃2.5, 5, and 6.5 times more than that predicted from the circular case in the same frequency band.

The distribution of BBHs with frequency, and thus the enhancement with respect to the circular channel, depend on the initial distribution of (a₀, r_p0). For more realistic estimates, in Sections 3.2 and 3.2 we recompute the equilibrium distributions for eccentric BBHs arising from triple systems and dynamical interactions in GCs, respectively.

Binaries in hierarchical triple systems are subject to gravitational perturbations from their tertiaries, and can be driven to high eccentricities due to the LK mechanism (Kozai 1962; Lidov 1962). In secular calculations where both the inner and outer orbits are averaged, the time over which the angular momentum of the inner binary is changed by order unity by the tertiary (the instantaneous LK timescale), in the ${m}_{2}\to 0$ limit, is given by (e.g., Bode & Wegg 2014; Antognini 2015)

$$\begin{eqnarray}&&{t}_{\mathrm{LK}}^{(\mathrm{ins})}\sim \displaystyle \frac{8\sqrt{2}}{15\pi }\left(1+\displaystyle \frac{M}{{m}_{3}}\right)\displaystyle \frac{{P}_{\mathrm{out}}^{2}}{P}{\left(1-{e}_{\mathrm{out}}^{2}\right)}^{3/2}\sqrt{\displaystyle \frac{{r}_{p}}{a}},\end{eqnarray} \tag{ 12 }$$

where m₃ is the tertiary mass, and e_out and P_out are the eccentricity and the outer orbital period, respectively.

The equilibrium argument in Section 2 relies on the assumption of dynamically isolated binaries whose frequencies evolve monotonically as a result of GW emission. This is only true for triple systems whose inner binaries are dynamically decoupled from the outer tertiary. A criterion for decoupling is that the inner binary is driven to sufficiently high eccentricity that ${t}_{\mathrm{LK}}^{(\mathrm{ins})}$ becomes longer than t_gw or the general relativistic (GR) precession timescale ${t}_{\mathrm{prec}}\sim 2{c}^{2}{a}^{3/2}{r}_{p}/[3{({GM})}^{3/2}]$.

In order to make a first estimate of the BBH population produced by triple systems, we run a secular calculation for triple systems with masses m₁ = m₂ = m₃ = 30 M_⊙. The eccentricities of the inner and outer orbits, e and e_out, are both sampled from a thermal distribution. The semimajor axis of the inner orbit a is sampled from a log-uniform distribution in [10,1000] au, and the semimajor axis ratio of outer to inner orbit, a_out/a is sampled from a log-uniform distribution in [10, 1000]. We discard systems with ${a}_{\mathrm{out}}(1-{e}_{\mathrm{out}})\lt 10a$ to make sure the validity of the secular calculation, and discard systems with a_out > 10⁵ au because they are too wide to make up an important fraction of triple systems. The orientations of both the inner and outer orbits are sampled randomly. We turn on the quadrupole-order term in the Newtonian three-body perturbing Hamiltonian, which leads to the LK effect, the 1PN term (GR precession), and the 2.5PN GW dissipation terms for the inner orbit. We run 10⁷ systems for 10 Gyrs, and find that ∼1.14% systems experience a decrease in their semimajor axis of order unity. These systems dynamically decouple from the tertiary and will merge within a relatively short time. We take the last eccentricity maximum of each such system and its semimajor axis a, and use them to set the initial conditions $({a}_{0},{r}_{p0})$ of "isolated binaries" in our calculation of the equilibrium distribution of the population.

Following the same procedure as in Section 3.1, we obtain the peak frequency histograms as shown in Figure 3. We again find a significant enhancement of the eccentric BBH population in the frequency range 0.1–1 mHz, in which the absolute number of systems with orbital periods within P ≤ 1, 10, and 100 days is ≃20, 80, and 130. Compared to Figure 2, the systems at f_p ≲ 0.05 mHz disappear. This comes from the fact that such systems still have perturbations from their tertiaries and are thus not dynamically isolated. However, the enhancement between 0.1 and 1 mHz persists, and is a factor of ≃10 at ∼0.3 mHz.

Zoom In Zoom Out Reset image size

Note that many triple systems undergoing LK oscillations may emit GWs in the ≃0.1–10 mHz band, but may not be dynamically decoupled by our criterion. Such systems are not included here because they do not obey the equilibrium assumptions set out in Section 2. The total population of GW emitters in the 0.1–10 mHz band (whether dynamically decoupled or not) has yet to be computed for a realistic and evolving distribution of massive triple systems as the Galaxy forms over cosmic time.

Figure 3 gives just one minimal estimate for the distribution from triple systems. Different component binary masses, which lead to octupole-order terms in the three-body Hamiltonian, tertiary masses, initial orbital parameter distributions, and cuts on the resulting population can quantitatively affect the results. Several variations are presented in the Appendix, with maximum and minimum enhancements relative to the circular case of ≃2–15.

BBHs arising in GCs and other dense stellar environments provide another important channel for eccentric BH migration. Recent numerical studies show that three populations of BBHs are produced during few-body scattering in GCs. One is produced by chaotic three-body motion, leading to BBH mergers in the cluster at very high eccentricity such that f_p ≳ 1 Hz and t_gw becomes less than the orbital period P (as in Equation (10)). These systems evolve dynamically and never enter the 0.1–1 mHz band. A second physical class of mergers are those from BBHs excited to high enough eccentricity within the cluster that t_gw becomes shorter than the time between two interactions. The third class is those BBHs ejected from the cluster. Adopting the nomenclature of (Samsing & D'Orazio 2018) we refer to these three classes as "three-body mergers," "two-body mergers," and "ejected mergers," respectively. The latter two classes evolve secularly through the 0.1–1 mHz band and are dynamically isolated (Samsing & D'Orazio 2018; Kremer et al. 2019).

We set aside the dynamically merging three-body mergers and consider only two-body mergers and ejected mergers. As an illustration, for these two physical categories, we adopt the distribution of BBH parameters resulting from the semi-analytic model described in Samsing & D'Orazio (2018). The binary component masses are assumed equal, with m = 30 M_⊙. The semimajor axis and eccentricity distributions for the BBHs when they are dynamically isolated are given by Samsing & D'Orazio (2018), and are used to set the initial conditions of "isolated binaries" in our calculation of the equilibrium BBH distribution, just as in Section 3.1.

Figure 4 shows the peak frequency histogram (top) and the ratio of eccentric to circular systems (bottom) for dynamically formed eccentric BBHs, normalized to the LIGO rate, as in previous figures. The shape of the number of systems per bin encodes the formation channel. As more clearly shown in the ratio plot (bottom), the two-body BBH mergers inside GCs dominate the distribution from 0.4–10 mHz, producing a peak relative to the circular case at f_p ≳ 1 mHz. The ejected BBH mergers resulting from binaries kicked out of GCs contribute to a much wider range of frequencies, with a low-frequency cutoff at ≃0.02 mHz for the P ≳ 10 day binaries.

Zoom In Zoom Out Reset image size

To estimate the two-detector sky-averaged signal-to-noise ratio (S/N) for eccentric BBHs, we start from Equation (45) in Barack & Cutler (2004), i.e., summing over contributions from all harmonics of the orbital frequency f_orb,

$$\begin{eqnarray}&&{\rm{S}}/{{\rm{N}}}^{2}=2\displaystyle \sum _{n=1}^{{n}_{\max }}\int \displaystyle \frac{{h}_{c,n}^{2}}{{f}_{n}{S}_{h}({f}_{n})}d\mathrm{ln}{f}_{n},\end{eqnarray} \tag{ 13 }$$

where n_max is the maximum harmonic used in fitting, ${f}_{n}\equiv {{nf}}_{\mathrm{orb}}$, S_h is the full strain spectral sensitivity density including the LISA instrumental noise and the confusion noise from the unresolved galactic binaries (e.g., Robson et al. 2018).⁷ The characteristic amplitude ${h}_{c,n}$ is given by ${h}_{c,n}={(\pi D)}^{-1}\sqrt{2{\dot{E}}_{n}/{\dot{f}}_{n}}$, where the unit G = c = 1 has been applied. D is the distance of the source. ${\dot{E}}_{n}$ is the GW energy emission rate in the nth harmonic and is given by

$$\begin{eqnarray}&&{\dot{E}}_{n}=\displaystyle \frac{32}{5}{M}_{c}^{10/3}{\left(2\pi {f}_{\mathrm{orb}}\right)}^{10/3}{g}_{n}(e),\end{eqnarray} \tag{ 14 }$$

where ${M}_{c}\equiv {({m}_{1}{m}_{2})}^{3/5}{M}^{-1/5}$ is the chirp mass (for ${m}_{1}={m}_{2}=30\,{M}_{\odot }$, ${M}_{c}=26\,{M}_{\odot }$), g_n(e) is given by Equation (20) in Peters & Mathews (1963). Since the eccentric migration timescale ($\sim {t}_{\mathrm{gw}}$) is much longer than the mission lifetime τ, the orbital frequency does not change much during the mission, amd the integral in the S/N becomes

$$\begin{eqnarray}&&\int \displaystyle \frac{{h}_{c,n}^{2}}{{f}_{n}{S}_{h}({f}_{n})}d\mathrm{ln}{f}_{n}\simeq \displaystyle \frac{{h}_{c,n}^{2}}{{f}_{n}{S}_{h}({f}_{n})}\displaystyle \frac{{\dot{f}}_{n}\tau }{{f}_{n}}=\displaystyle \frac{2{\dot{E}}_{n}\tau }{{\pi }^{2}{D}^{2}{f}_{n}^{2}{S}_{h}({f}_{n})}.\end{eqnarray} \tag{ 15 }$$

Combining it with Equation (14) we obtain

$$\begin{eqnarray}&&{\rm{S}}/{{\rm{N}}}^{2}=\displaystyle \frac{512\tau }{5{D}^{2}}\displaystyle \frac{{\left({{GM}}_{c}\right)}^{10/3}}{{c}^{8}}{\left(2\pi {f}_{\mathrm{orb}}\right)}^{4/3}\displaystyle \sum _{n=1}^{{n}_{\max }}\displaystyle \frac{{g}_{n}(e)}{{n}^{2}{S}_{h}({f}_{n})}\,.\end{eqnarray} \tag{ 16 }$$

As discussed by Gould (2011) in the context of eccentric binary white dwarfs, the GW emission is dominated by pulses at periastron.

Taking τ = 4 yr, D = 10 kpc, n_max = 10⁵, f_p = 0.5 mHz, M_c = 26 M_⊙, for an eccentric BBH with P = 1, 10, and 100 days, we obtain S/N = 26, 8.5, and 2.7, respectively, implying that a number of eccentric BBHs will be detectable by LISA. In Figure 5, we show the S/Ns for systems at 10 kpc with different peak frequencies and P = 1, 10, and 100 days. Requiring S/N = 5(2) for detection and assuming a distance of 10 kpc (which gives a conservative estimate of the number of observable systems), we estimate that ∼7, 7, and 8 (15, 17, and 17) eccentric BBHs with P less than 1, 10, and 100 days will be detected in our Galaxy in 0.1–1 mHz in the case of the simplest distribution, while ∼5, 6, 7 (11, 14, 15) systems in the triple case, and ∼4, 4, 5 (8, 9, 10) systems in the GC case. While the number of detectable eccentric systems is largely limited by the S/N at f_p ≲ 0.2 mHz, where many more systems exist, the detection of the few systems at higher frequencies immediately implies the existence of the eccentric BBH population. Also note that although less systems are present at higher frequencies in our Galaxy, a larger volume is accessible due to the larger S/Ns, hence more extragalactic sources may be detectable (e.g., Rodriguez et al. 2018; Samsing & D'Orazio 2018; Kremer et al. 2019). These sources could also be interesting to experiments sensitive to slightly higher frequency bands, such as DECIGO, Taiji, and TianQin. The right axis of Figure 5 shows that eccentric systems with P < 1 day can be discovered by LISA out to ∼8 Mpc distances.

Zoom In Zoom Out Reset image size

In Table 1, we show the estimated numbers of LISA-detectable galactic eccentric and circular BBHs, assuming different LIGO rates due to the uncertainty of the measured LIGO merger rate. Although the actual number of circular systems in each frequency bin is less than that of eccentric systems, the detectable numbers may be similar, due to the better S/Ns when detecting circular binaries. Note that the numbers are obtained by assuming all the systems are at distance 10 kpc, a typical value of their average distance. While closer systems have a higher S/N, much less systems exist at closer distances if a uniform spatial distribution of BBHs in the stellar halo is assumed, hence resulting in little changes in the estimated numbers of observable systems.

A primary assumption in the results presented here is that the numbers of BBH systems are based on the equilibrium assumption, which will break down on 10 Gyr timescales. However, most of the eccentric migrating systems in triples and GCs with f_p > 0.1 mHz have merger times less than 1 Gyr, during which the variations expected as a result of the time history of star formation in the Galaxy may not be significant.

Additional uncertainties lie in the merger rate. There is a factor of a few uncertainty in the LIGO BBH merger rate, which is a function of BBH mass, and we have also assumed that the LIGO BBH merger rate applies to our Galaxy. While these uncertainties will affect the absolute numbers of systems expected, the ratios between the eccentric and circular cases presented are robust. In reality, the merger rate may be a mixture of all the possible channels. For the GC case, most of cluster simulations predict merger rates less than 50 Gpc⁻³ yr⁻¹ (e.g., Rodriguez & Loeb 2018). Thus, the combined number of eccentric BBHs and their frequency distribution may depend on the fraction of each channel, which might in turn be used to probe the relative importance of various channels by LISA.

For the case of eccentric BBHs produced in triple systems, there are a number of uncertainties and complexities. These include (1) octupole-order perturbations for a realistic population, (2) evection, and (3) the astrophysics of realistic triple system masses and their evolution. The octupole-order perturbation (1) arises when the inner binary has unequal masses (e.g., Naoz et al. 2013a) and may change the frequency distributions as indicated in the tests in the Appendix. A more comprehensive analysis with a realistic BH mass distribution is needed to fully explore its effect. Evection (2) induces eccentricity oscillations of the inner orbit on timescale of P_out, which may be important at high eccentricities where the inner orbit has small angular momentum, and is thus prone to torque from the tertiary (e.g., Ivanov et al. 2005; Katz & Dong 2012; Antognini et al. 2014; Fang et al. 2018). However, we neglect it here because evection may be considered as random kicks to the inner orbit and will produce a thermal distribution of e, which is already assumed in our initial distribution. Thus, while evection may change the absolute number of triple systems that produce merging and dynamically isolated BBHs, and while individual systems may experience non-secular changes of r_p, the overall distribution should not be modified by evection. The inclusion of evection may also introduce non-negligible secular effects when the triple systems are moderately hierarchical (e.g., Luo et al. 2016b; Grishin et al. 2018; Lei et al. 2018), which could suppress the octupole-order oscillations, and thereby affect the resulting distribution of systems as a function of GW frequency. Finally, (3) we neglect stellar evolution of the massive star progenitors, including possible mass transfer, adiabatic and dynamical mass-loss, and natal kicks due to the recoil during asymmetric supernova explosions. These effects may change the orbital parameter distributions or unbind the systems, which may suppress this formation channel (e.g., Silsbee & Tremaine 2017). Additionally, as the Appendix shows (see Run 3), a realistic distribution of tertiary masses can significantly change the relative enhancement of eccentric to circular systems.

Our major findings in this paper are as follows.

• 1.  
  Assuming the distribution of BBHs is in equilibrium, producing a steady merger rate as seen by LIGO, we show that eccentric BBH formation channels predict a larger number of systems relative to circular BBH formation channels throughout the 0.1–1 mHz GW frequency band. Equations (8) and (9) and Figure 1 show that this predicted enhancement is generic, and follows from the fact that eccentric systems spend more time radiating in the low-frequency GW band than circular systems.
• 2.  
  We estimate the absolute number of radiating systems in the circular and eccentric cases in the Galaxy. Figures 2–4 show the eccentric and circular cases for a generic eccentric population, for eccentric BBHs produced by triple systems undergoing Lidov–Kozai oscillations, and BBHs formed dynamically in GCs, respectively. Assuming that both eccentric and circular channels produce the observed LIGO rate, we find that eccentric systems outnumber circular systems by a factor of 2–10 throughout the 0.1−10 mHz GW band. Under these assumptions, there are ∼10–100 eccentric BBHs with GW peak frequencies of 0.1–1 mHz in the Galaxy. Dozens of these systems have orbital periods of the order of days to months.
• 3.  
  Eccentric BBH systems emit GW pulses at periastron. We calculate the S/N for detecting eccentric binaries with a LISA-like sensitivity curve, and estimate that ≃7 (15) eccentric systems should be seen in the Galaxy with S/N ≥ 5(2), slightly less than the number of observable circular systems (11 and 16 for S/N ≥ 5 and 2) in the range of 0.1–1 mHz. See Figure 5 and Table 1. For the rarer eccentric systems with higher peak GW frequency of 1–10 mHz, the detection volume increases to 0.8–8 Mpc for systems with orbital periods less than ≃1 day.