DISTINGUISHING BETWEEN FORMATION CHANNELS FOR BINARY BLACK HOLES WITH LISA

The era of gravitational wave (GW) astrophysics began with the discovery of the binary black hole (BBH) merger, GW150914, by Advanced LIGO (aLIGO; Abbott et al. 2016b). The subsequent detection of a BBH merger, GW151226, with smaller progenitor masses, suggests diversity in the potential formation channels of the BBHs (Abbott et al. 2016c). Parameter estimation analyses of GW150914 were done assuming a circular orbit, but eccentricities of $e\lesssim 0.1$ would not have been measurable for the event (Abbott et al. 2016a). Multiple scenarios have been proposed to produce GW150914-like BBHs (Belczynski et al. 2016; Bird et al. 2016; Mandel & de Mink 2016; Marchant et al. 2016; Rodriguez et al. 2016b). Despite the radically different environments of the proposed formation channels, aLIGO is not able to discriminate between them without a full population of measured spins (Vitale et al. 2015), mass ratios, redshift distributions, or eccentricities.

BBHs formed dynamically in dense stellar environments such as globular clusters (GCs) generally exhibit higher eccentricities than BBHs evolved in isolated environments in galactic fields. Binary evolutionary processes like tidal circularization and mass transfer affect the orbital parameters of BBHs formed in galactic fields (e.g., Postnov & Yungelson 2014, and references therein). Once the BBHs are formed, GW emission shrinks and circularizes the initial orbit over long timescales, leading to low eccentricities in both the Laser Interferometer Space Antenna (LISA) and aLIGO bands. In GCs, BBHs are formed dynamically and subsequently ejected from the cluster, with eccentricities distributed thermally after initially sinking to the cluster core (Morscher et al. 2013). The distribution of BBH separations is a function of GC mass and leads to BBHs that evolve through GW emission to the LISA frequency band with non-negligible eccentricities (Rodriguez et al. 2016a).

It has recently been suggested that hundreds of BBHs similar to the progenitor of GW150914 will be observable by a space-based GW detector like LISA and subsequently evolve into the aLIGO band within 10 yr (Sesana 2016). It has also been shown that BBHs detected by LISA may also have measurable eccentricity if $e\gtrsim 0.001$ (Nishizawa et al. 2016; Seto 2016). Thus, using eccentricity, LISA may be able to discriminate between BBH formation channels.

In this letter we show that BBHs formed in galactic fields and BBHs formed in GCs have different combinations of mass and eccentricity in the LISA frequency band, thus providing a method to distinguish between the two formation channels. We also show two distinct populations within the BBHs evolved in isolation that are separated by higher and lower eccentricities that will be potentially distinguishable by LISA.³ In Section 2 we give a brief overview of our simulated BBH populations. In Section 3 we describe the eccentricity distribution of the simulated BBHs, show how each population of BBHs evolves due to GW emission and suggest formation channels potentially distinguishable by LISA. We provide concluding remarks in Section 4.

We simulated two populations of BBHs, one evolved in isolation in galactic fields and one formed dynamically in GCs, to explore the evolution of each population from formation frequencies of $f\sim {10}^{-9}\mbox{--}{10}^{-7}$ Hz all the way to the upper end of the aLIGO band at $f\sim {10}^{3}\,\mathrm{Hz}$. The same binary evolution physics was used to create both population simulations with stellar dynamics effects turned on and off.

We used a modified version of the binary stellar evolution code BSE (Hurley et al. 2002) with updated stellar evolution prescriptions including BH formation and natal kicks (Fryer & Kalogera 2001; Belczynski et al. 2002), updated metallicity-dependent stellar-wind prescriptions (Vink et al. 2001; Vink & de Koter 2005; Belczynski et al. 2010), and additional mechanisms to account for fallback in neutrino-driven supernovae (Fryer et al. 2012). For this study, we adopt the fiducial models of Rodriguez et al. (2016a), which are based on the most recent stellar evolution prescriptions for galactic fields (Dominik et al. 2013).

Our population of dynamically formed binaries is taken from a collection of 48 GC models developed in Rodriguez et al. (2016a) with the Cluster Monte Carlo (CMC) code (Pattabiraman et al. 2013, and the  references therein). CMC employs a statistical approach to stellar dynamics, first developed by Hénon (1975), which enables the modeling of star clusters with a significantly greater number of particles than a direct N-body simulation, while still employing the necessary physics—single and binary stellar evolution with BSE, three-body and four-body strong gravitational encounters (Fregeau & Rasio 2007), and dynamical three-body binary formation (Morscher et al. 2013)—to fully characterize the dynamical BBH merger problem. We neglect long-term secular effects (e.g., Antonini et al. 2014, 2016) and relativistic dynamical scattering (e.g., Samsing et al. 2014), which are expected to contribute to the overall BBH population at the $\sim 1 \%$ or lower level, according to these studies.

These 48 GC models span a range of initial particle numbers ($N=2\times {10}^{5},\,5\times {10}^{5},\,1\times {10}^{6},\,\mathrm{and}\,2\times {10}^{6}$), initial virial radii (${R}_{v}=1,2\,\mathrm{pc}$), and stellar metallicities ($Z=0.25{Z}_{\odot },0.05{Z}_{\odot },\,\mathrm{and}\,0.01{Z}_{\odot }$), with two statistically independent models generated for each set of initial conditions. As was done in Rodriguez et al. (2016a), we assume that the GC population of the local universe is comprised of $\sim 44 \%$ high-metallicity GCs ($0.25{Z}_{\odot }$), and $\sim 56 \%$ low-metallicity GCs ($0.05{Z}_{\odot }$ and $0.01{Z}_{\odot }$). Furthermore, we assign to each GC a tidal radius based on its galactocentric distance, which we assume to be correlated to its stellar metallicity based on observations of the Milky Way and other galaxies (Harris 2010). We finally assume a log-normal GC mass function, based on recent observations of the GC luminosity function in brightest-cluster galaxies (Harris et al. 2014) and a mass-to-light ratio of 2 for old stellar systems (e.g., Bell et al. 2003).

We draw a sample of BBHs ejected from our 48 GC models by randomly selecting binaries from each model. The number of binaries selected from a given model is determined by its weight, which we assign by dividing the GC mass function into bins, the midpoints of which are set by the average mass of GC models with the same initial particle number. The integral of the GC mass function over that bin then determines the weight assigned to that model. In other words, GC models with larger masses ($3\times {10}^{5}\,{M}_{\odot }$ to $6\times {10}^{5}\,{M}_{\odot }$, corresponding to models with initial particle numbers of 1 × 10⁶ and 2 × 10⁶) contribute more binaries than clusters with smaller masses (see Rodriguez et al. 2016a for details). Once we select this population of binaries from all clusters, we generate a five-dimensional Gaussian kernel density estimate (KDE) from the formation masses, separation, eccentricity, and ejection times of the BBH population. We sample 1000 BBHs from the KDE that evolve to the LISA band at $f={10}^{-5}$ Hz in the last gigayear before the present time. We assume that all GCs were formed 12 ± 1 Gyr ago, and we assign to each binary a cluster birth time drawn from a similar Gaussian distribution. The modulated cluster ages allow us to account for the observed spread in ages of GCs (e.g., Correnti et al. 2016). We require the dynamically formed BBHs to enter the LISA band in the last gigayear before the present time.

We generate galactic field populations from four sets of 10⁶ initial binaries sampled from standard probability distribution functions to assign each binary with an initial metallicity (Z), primary mass (m), mass ratio (q), orbital separation (a), and eccentricity (e). The population models compare four different metallicities $(Z=1{Z}_{\odot },\ 0.25{Z}_{\odot },\ 0.05{Z}_{\odot },$ and $0.01{Z}_{\odot })$, with the sub-solar metallicities being consistent with the metallicities used in the GC models. For the initial primary mass we adopt the stellar IMF $\xi (m)\propto {m}^{-2.3},m\geqslant 1{M}_{\odot }$ (Kroupa 2001) with a primary mass limit of 150 M_⊙. We assume a uniform initial mass ratio distribution consistent with current observational constraints (Mazeh et al. 1992; Goldberg & Mazeh 1994; Kobulnicky et al. 2014). We assume that initial orbital separations are distributed uniformly in $\mathrm{log}(a)$ at wide separations ($10\,{R}_{\odot }\leqslant a\leqslant 5.75\times {10}^{6}\,{R}_{\odot }$) and fall off linearly at small separations as $\zeta {(a)\propto (a/{a}_{0})}^{1.2}$, $a\lt 10\,{R}_{\odot }$ (Han 1998). The initial eccentricities are distributed thermally (Heggie 1975) as $\eta (e)=2e$.

We evolve the galactic field population for 13.87 Gyr using the same binary evolution models as the GC population, creating an equivalent population to the GC population but without dynamics. We log the birth parameters of each BBH, including important formation processes like the number of common-envelope episodes and the natal kicks imparted to the binary from the birth of each black hole. Just as with the dynamically formed BBHs, we require the low-metallicity ($Z\lt {Z}_{\odot }$) BBHs to enter the LISA band in the last gigayear before the present time. We retain any solar metallicity BBHs that evolve to the LISA band over the last 10 Gyr.

Due to the eccentric nature of BBHs formed in GCs, it is useful to consider the frequency of GWs emitted at higher harmonics of a binary's orbital frequency. The frequency of maximum GW power emission from an eccentric binary is estimated as (Wen 2003)

$$\begin{eqnarray}&&{f}_{\mathrm{GW}}=\displaystyle \frac{\sqrt{G({m}_{1}+{m}_{2})}}{\pi }\displaystyle \frac{{(1+e)}^{1.1954}}{{[a(1-{e}^{2})]}^{1.5}}.\end{eqnarray} \tag{ 1 }$$

If we consider the peak GW frequency of Equation (1) for BBHs formed in GCs, we find that the GW frequencies of the population are substantially higher than the circular GW frequencies. Since the peak sensitivity of LISA falls near $f\sim {10}^{-3}\mbox{--}{10}^{-2}$ Hz, shifts to higher frequency aid in the detectability of these sources.

3.1. BBH Orbital Evolution

Since the dynamically formed BBHs are ejected from their host GC and the galactic field BBHs evolve in isolation, only GW emission will affect the evolution of each binary. Using the simulated populations from Section 2, we evolve each BBH evolved in isolation from its birth frequency to 10³ Hz and each dynamically formed BBH from the time of ejection from its host cluster to 10³ Hz using the quadrupole approximated GW orbital evolution equations (Peters 1964).

BBHs with frequencies larger than 10⁻³ Hz are expected to have measurable frequency evolution due to a GW emission known as the GW chirp. The chirp mass, ${{ \mathcal M }}_{c}={({m}_{1}{m}_{2})}^{3/5}/{({m}_{1}+{m}_{2})}^{1/5}$, and eccentricity (or lack thereof) can be measured for every BBH with a measured chirp. Figure 1 illustrates the evolution of eccentricity and GW frequency for each binary in the models from each formation channel. The minimal eccentricities measurable by a LISA-like detector are shown in red (Nishizawa et al. 2016). The LISA frequency range is highlighted in gray and the frequency range where BBHs with a chirp mass of ${{ \mathcal M }}_{c}\gtrsim 6{M}_{\odot }$ have measurable chirps is highlighted in light blue.

Zoom In Zoom Out Reset image size

In addition to dynamically formed BBHs, we consider two populations of BBHs formed in isolation: those including either a single or no common-envelope episodes (hereafter 1CE and 0CE). BBHs formed both dynamically (black) and in isolation (green, blue) have measurable eccentricities and frequencies above 10⁻³ Hz, with dynamically formed binaries having larger eccentricities than those in galactic fields. Above 10⁻² Hz, only BBHs formed in GCs fall above the e = 0.01 line. This suggests that any BBH detected above frequency ${10}^{-2}$ Hz with eccentricity $e\geqslant 0.01$ formed dynamically in a dense stellar environment.

3.2. Eccentricity Distributions

Dynamically formed BBHs are ejected from the cluster with a thermal eccentricity distribution. Once ejected, BBHs evolve only through the emission of GWs that circularize and shrink the binary orbit. This leads to BBHs with eccentricities of $e\gtrsim 0.1$ on the low end of the LISA frequency band (${f}_{\mathrm{GW}}\sim {10}^{-5}$ Hz) and eccentricities of $e\gtrsim 0.001$ on the high end of the LISA frequency band (${f}_{\mathrm{GW}}\sim {10}^{-2}$ Hz). BBHs formed in isolation generally form at lower eccentricities, but also with lower chirp masses than BBHs formed in GCs.

Figure 2 shows the cumulative distribution of the eccentricities of the simulated BBHs formed dynamically (black) and in isolation (blue, green) at different points in their orbital evolution. Highly eccentric dynamically formed binaries were recently ejected from the GC, while the less eccentric binaries have had more time to circularize through GW emission. At 10⁻³ Hz, 92% of dynamically formed BBHs have $e\geqslant 0.001$, while 32% of them have $e\geqslant 0.01$. At 10⁻² Hz, the eccentricity of the dynamically formed BBHs has decreased such that 30% have $e\geqslant 0.001$ and 7% have $e\geqslant 0.01$. For the 1CE BBHs, 91% have $e\geqslant 0.001$ and 23% have $e\geqslant 0.01$ at 10⁻³ Hz. Again, at 10⁻² Hz the eccentricity of the 1CE BBHs has decreased, with 19% having $e\geqslant 0.001$ and 0% with $e\geqslant 0.01$. There are no 0CE BBHs with measurable eccentricities at GW frequencies above 10⁻³ Hz.

Zoom In Zoom Out Reset image size

3.3. Chirp Mass and Eccentricity Correlations

A particularly useful way to separate BBHs formed dynamically or in isolation is through the correlations between the chirp mass and eccentricity of each population. BBHs with chirp masses  ${{ \mathcal M }}_{c}\lt 10{M}_{\odot }$ exclusively form in isolation in galactic fields, though we note that young, high-metallicity clusters (not included in our models) are capable of producing BBHs with lower chirp masses (Chatterjee et al. 2016).

If only the chirp mass of a BBH with ${{ \mathcal M }}_{c}\gt 10{M}_{\odot }$ is observed, it is impossible to discern which population it originated in since the chirp masses of the GC, 1CE, and 0CE populations overlap in this region. However, if the eccentricity is also measured, the three populations can be resolved. Figure 3 shows the eccentricity versus chirp mass plots of each population at ${f}_{\mathrm{GW}}={10}^{-3}$ Hz. The shape of the distributions stays constant but the eccentricity decreases as the BBHs evolve to higher frequencies through GW emission. For each population, the estimation error on the chirp mass is ${\rm{\Delta }}{{ \mathcal M }}_{c}/{{ \mathcal M }}_{c}\simeq 3({\rm{\Delta }}\dot{f}/\dot{f})/5$ where ${\rm{\Delta }}\dot{f}=0.43{(\mathrm{SNR}/10)}^{-1}{T}_{\mathrm{obs}}^{-2}$ (Takahashi & Seto 2002). In all cases, the chirp mass estimation error is smaller than the width of the data points.

Zoom In Zoom Out Reset image size

The chirp masses of the GW150914 and GW151226 progenitors are plotted in Figure 3. The GW150914 and LVT151012 progenitors are consistent with both dynamical and isolated formation channels, and the GW151226 progenitor is consistent with the 1CE isolated formation channels. If the GW150914, LVT151012, or GW151226 progenitors had been observed by LISA, an eccentricity measurement (or lack thereof) could have aided in identifying their formation histories.

3.4. Multiple Field Populations

We find that the population of isolated binaries that evolve into the LISA frequency band is split into higher- and lower-eccentricity populations. The higher-eccentricity population is comprised entirely of BBHs that have undergone a single common-envelope episode, while the lower-eccentricity population is comprised of a combination of 0CE and 1CE BBHs. The higher-eccentricity 1CE population, while still distinct from the GC BBHs, has potentially measurable eccentricities.

The natal kick imparted to a BH at birth is connected to kicks thought to be imparted to neutron stars (NSs). NS natal kicks are assumed to follow a Maxwellian distribution with a dispersion of 265 km s⁻¹ (Hobbs et al. 2005). The BH natal kick is modified to be dependent on the mass of the pre-collapse stellar core through fallback processes, with higher kicks imparted to lower-mass remnants (Fryer et al. 2012). A binary can be disrupted if the overall energy imparted through the natal kicks is higher than the binding energy of the binary orbit.

In the case of the higher-eccentricity population, the 1CE binaries are driven to small separations through the common-envelope mechanism. Since the masses of the 1CE BBH components are generally low due to low fallback, the natal kick speeds are high enough to produce higher eccentricities, but not so high as to disrupt the BBHs. This suggests that the eccentricity of the 1CE population is dependent on the natal kick physics of the formation of the second BH.

We simulated three additional populations with varying natal kick prescriptions to explore the effect of natal kicks on the eccentricity of BBHs formed in isolation. The first prescription sets the BH natal kick equal to zero, thus any orbital changes imparted to the binary are due to momentum conservation from instantaneous mass-loss in the BH formation. We also include two variants of the BH natal kick prescription: one using standard NS kicks for BHs and one that modifies the NS kick by the mass fraction $({M}_{\mathrm{NS}}/{M}_{\mathrm{BH}})$, where ${M}_{\mathrm{BH}}$ is the mass of the newly formed BH and we assume ${M}_{\mathrm{NS}}=1.4\ {M}_{\odot }$.

We compare our results for four BH natal kick prescriptions: no natal kicks, fallback modulated kicks, fractional NS kicks weighted by mass, and full NS kicks. We did not generate new GC simulations for each kick prescription, since BBHs formed dynamically in dense stellar environments form with thermally distributed eccentricities, losing the memory of kick effects on the eccentricity distribution.

Figure 4 plots the eccentricity versus chirp mass for the 0CE and 1CE formed in isolation for the four natal kick prescriptions mentioned above at their birth orbital frequency and at 10⁻³ Hz. Again, in all cases, the chirp mass estimation error bars are smaller than the width of the data points. The lower-mass BBHs in all cases have preferentially higher kicks and thus higher birth eccentricities. The 0CE BBHs generally have larger birth separations and masses, which results in low eccentricities by the time the BBH evolves through GW emission to the LISA frequency band.

Zoom In Zoom Out Reset image size

As expected, BBHs formed with full NS natal kicks retain the largest eccentricities through their GW driven orbital evolution, followed by fractional NS kicks, then fallback modulated kicks and no natal kicks. More massive BBHs evolve to lower eccentricities at a given orbital frequency than the lower-mass systems in the same population because of the mass dependence of orbital frequency. This is responsible for diminishing the high eccentricities found at birth in the 0CE population for both the full and fractional NS kicks once they reach the LISA band.

We have shown that BBHs formed both dynamically and in isolation may have measurable eccentricity in the LISA band. If BBHs are detected by LISA with eccentricities of $e\gtrsim 0.01$ at frequencies above 10⁻² Hz, they will have almost certainly originated from dynamical processes in old, dense stellar environments. If BBHs with eccentricities of $e\gtrsim 0.01$ and chirp masses of ${{ \mathcal M }}_{c}\lesssim 10{M}_{\odot }$ are detected by LISA at low frequencies, they likely originated from a common-envelope formation scenario.

In the future, we plan to extend this study by implementing a full treatment of the formation-redshift distribution of the BBHs observable by LISA originating from both dynamical processes in dense stellar environments and isolated binary evolution in galactic fields. We also plan to properly account for the detectability of each BBH using eccentricity-dependent signal-to-noise ratios (K. Breivik et al. 2016, in preparation).

K.B. and S.L.L. acknowledge support from NASA grant NNX13AM10G. C.R. and F.A.R. acknowledge support from NSF grant AST-1312945 and from NASA grant NNX14AP92G. V.K. acknowledges support from NSF grant PHY-1307020/002 and from Northwestern University. V.K. and F.A.R. also acknowledge support from NSF grant PHY-1066293 at the Aspen Center for Physics.