Probing Massive Black Hole Binary Populations with LISA

TL;DR

This study uses cosmological simulations to predict LISA's detection rates of merging massive black hole binaries, analyzing how different evolutionary models influence these rates and properties.

Contribution

It introduces a novel analysis of LISA detection prospects using Illustris simulation data and compares binary evolution models affecting detection predictions.

Findings

LISA detection rates of ~0.5-1 per year for black holes >10^5 M_sun

Detection rates are likely underestimated due to mass range limitations

Different binary evolution models significantly impact detection predictions

Abstract

ESA and NASA are moving forward with plans to launch LISA around 2034. With data from the Illustris cosmological simulation, we provide analysis of LISA detection rates accompanied by characterization of the merging massive black hole population. Massive black holes of total mass $\sim10^5-10^{10} M_\odot$ are the focus of this study. We evolve Illustris massive black hole mergers, which form at separations on the order of the simulation resolution ($\sim$kpc scales), through coalescence with two different treatments for the binary massive black hole evolutionary process. The coalescence times of the population, as well as physical properties of the black holes, form a statistical basis for each evolutionary treatment. From these bases, we Monte Carlo synthesize many realizations of the merging massive black hole population to build mock LISA detection catalogs. We analyze how our massive black hole binary evolutionary models affect detection rates and the associated parameter distributions measured by LISA. With our models, we find massive black hole binary detection rates with LISA of $\sim0.5-1$ yr$^{-1}$ for massive black holes with masses greater than $10^5M_\odot$. This should be treated as a lower limit primarily because our massive black hole sample does not include masses below $10^5M_\odot$, which may significantly add to the observed rate. We suggest reasons why we predict lower detection rates compared to much of the literature.

Tables

Table 1. Table 1: Merger and detection rate calculations per year are shown using the integral calculation (Equation 27 ). The evolution prescriptions are listed in the first column. The top row shows the “no delays” model: ND. The second row shows the ND-6 model, displayed in italics because it represents a subset of the ND model with m 1 , m 2 ≥ 10 6 ​ M ⊙ subscript 𝑚 1 subscript 𝑚 2 superscript 10 6 subscript 𝑀 direct-product m_{1},m_{2}\geq 10^{6}M_{\odot} . The final two rows show the DA17 (Dosopoulou & Antonini, 2017 ) and K17 (Kelley et al., 2017a ; Kelley et al., 2017b ) models, respectively. The merger rate gives the rate of coalescences without considering LISA detectability. The remaining columns show the detection rate ( ρ > 8 𝜌 8 \rho>8 ). Additionally, detection rates are separated into signal types: All, Ins, MR, and BHS. “All” indicates reaching detection threshold using the entire signal. “Ins” and “MR” represent detection rates of inspiral signals and signals from the merger and ringdown, respectively. These categories are not independent: a single binary can add to the rate in both categories. “BHS” is the detection rate of MBH binaries where black hole spectroscopy is possible (see Equation 24 ). We estimate our errors in these predictions to all fall below 0.01 according to equation 17 in Salcido et al. ( 2016 ) .

Prescription	Merger Rate	All	Ins	MR	BHS
ND	0.98	0.75	0.45	0.77	0.74
ND-6	0.57	0.44	0.23	0.45	0.43
DA17	0.80	0.70	0.42	0.70	0.67
K17	0.55	0.44	0.28	0.45	0.43

Table 2. Table 2: A collection of papers quoting rates for the detection of MBH binaries by LISA is shown above. We focus here on papers analyzing the large-seed ( 10 4 − 10 6 ​ M ⊙ superscript 10 4 superscript 10 6 subscript 𝑀 direct-product 10^{4}-10^{6}M_{\odot} ) formation channels since this is similar to the seeding mechanism in Illustris. The base population tells if the models are based on SAMs or hydrodynamic simulations. Short descriptions of the MBH binary evolutionary prescriptions are also given. “None” indicates that no delays between galactic and MBH binary mergers were included. “Constant” indicates that a constant delay was used for all binaries. “DF” is a dynamical friction prescription. “LC” indicates inclusion of a stellar hardening or loss-cone scattering model. A prescription involving torque from a gas disc is indicated with “VD.” “GW” indicates a gravitational wave driven inspiral. A prescription containing triple MBH systems is expressed with “Tri.” 1 For detectable ringdown signals. For the mass ranges we consider in this paper, ringdown signals will be measurable for all detectable sources. 2 A constant delay was chosen based on gas-heavy versus gas-poor host galaxies.

Reference	Base Population	MBHB Evolution Prescription	Merger Rate (yr-1)
Arun et al. (2009)	SAM	None	$\sim 22$
Sesana et al. (2011)	SAM	None	$\sim 25$
Klein et al. (2016)	SAM	DF,LC,VD,GW,Tri	$\sim$8
Berti et al. (2016)1	SAM	DF,LC,VD,GW,Tri	$\sim$8
Salcido et al. (2016)	Hydrodynamic	Constant2	$\sim$2
Bonetti et al. (2019)	SAM	DF,LC,VD,GW,Tri	$\sim$23
This Paper	Hydrodynamic	DF,LC,VD,GW	$\sim 0.5-1$

Table 3. Table 3: Monte Carlo results for the observed detection rate per year from our 10000 catalogs are shown above (see 3.3 ). We focused our Monte Carlo calculation on the overall signal detection rate as well as sources only detected in their inspiral stage. The coalescence timescale prescriptions are listed in the first column: ND is the “no delay” model; ND-6 is a subset of the ND model with all constituent masses below 10 6 ​ M ⊙ superscript 10 6 subscript 𝑀 direct-product 10^{6}M_{\odot} eliminated from consideration; DA17 is the binary inspiral model from Dosopoulou & Antonini ( 2017 ) ; and K17 is the inspiral model from Kelley et al. ( 2017a ); Kelley et al. ( 2017b ) . The full signal detection rate with PL and CLLF are shown in the second and third columns, respectivley. Similary, inspiral-only detection rates are shown for PL and CLLF in the fourth and fifth columns, respectively. The last two columns show the fraction of the total detection rate contributed by inspiral-only signals using PL and CLLF. For our Monte Carlo results, the error in the predictions is approximately 1 / N = 1 / 10000 = 1 % 1 𝑁 1 10000 percent 1 \sqrt{1/N}=\sqrt{1/10000}=1\% .

Prescription	PL	CLLF	PL	CLLF	PL	CLLF
-	All Signal	All Signal	Ins Only	Ins Only	Fraction of Ins Only	Fraction of Ins Only
ND	0.79	1.08	<0.01	0.19	<1%	18%
ND-6	0.47	0.57	<0.01	0.05	<1%	9%
DA17	0.72	0.97	<0.01	0.20	<1%	21%
K17	0.47	0.64	<0.01	0.13	<1%	20%

Table 4. Table 4: This table is the same as Table 1 , with rates added for the CLLF LISA configuration. Merger and detection rate calculations per year are shown using the integral calculation (Equation 27 ). The evolution prescriptions are listed in the first column. The top row shows the “no delays” model: ND. The second row shows the ND-6 model, displayed in italics because it represents a subset of the ND model with m 1 , m 2 ≥ 10 6 ​ M ⊙ subscript 𝑚 1 subscript 𝑚 2 superscript 10 6 subscript 𝑀 direct-product m_{1},m_{2}\geq 10^{6}M_{\odot} . The final two rows show the DA17 (Dosopoulou & Antonini, 2017 ) and K17 (Kelley et al., 2017a ; Kelley et al., 2017b ) models, respectively. The merger rate gives the rate of coalescences without considering LISA detectability. The remaining columns are labeled with the sensitivity curve used (PL or CLLF) to determine detectability. Additionally, detection rates are separated into signal types: All, Ins, MR, and BHS. “All” indicates reaching detection threshold using the entire signal. “Ins” and “MR” represent detection rates of inspiral signals and signals from the merger and ringdown, respectively. These categories are not independent. A single binary can add to the rate in both categories. “BHS” is the detection rate of MBH binaries where black hole spectroscopy is possible (see Equation 24 ). We estimate our errors in these predictions to all fall below 0.01 according to equation 17 in Salcido et al. ( 2016 ) .

-	-	PL	CLLF	PL	CLLF	PL	CLLF	PL	CLLF
Prescription	Merger Rate	All	All	Ins	Ins	MR	MR	BHS	BHS
ND	0.98	0.75	0.89	0.45	0.75	0.77	0.89	0.74	0.88
ND-6	0.57	0.44	0.51	0.23	0.43	0.45	0.52	0.43	0.51
DA17	0.80	0.70	0.77	0.42	0.68	0.70	0.77	0.67	0.77
K17	0.55	0.44	0.50	0.28	0.44	0.45	0.51	0.43	0.50

Equations

$$r_{\text{infl}}=10.8\left(\frac{M}{10^8{M}_{\odot }}\right){\left(\frac{\sigma }{200{\text{ km s}}^{-1}}\right)}^{-2}\text{pc},$$

$$T_{\star }^{\text{bare}}=17\frac{6.6}{\ln \Lambda }{\left(\frac{R_e}{10\text{ kpc}}\right)}^2\left(\frac{\sigma }{300{\text{ km s}}^{-1}}\right)\left(\frac{10^8{M}_{\odot }}{m}\right)\text{Gyr},$$

$$T_{\star ,1}^{\text{gx}}=0.06\frac{2}{\ln {\Lambda }^{\prime }}{\left(\frac{R_e}{10\text{ kpc}}\right)}^2\left(\frac{\sigma }{300{\text{ km s}}^{-1}}\right)\left(\frac{10^8{M}_{\odot }}{m}\right)\text{Gyr},$$

$$T_{\star ,2}^{\text{gx}}=0.15\frac{2}{\ln {\Lambda }^{\prime }}\left(\frac{R_e}{10\text{ kpc}}\right){\left(\frac{\sigma }{300{\text{ km s}}^{-1}}\right)}^2\left(\frac{100{\text{ km s}}^{-1}}{{\sigma }_s}\right)\text{Gyr}.$$

$$T_{\star }=\text{max}\left(T_{\star ,1}^{\text{gx}},T_{\star ,2}^{\text{gx}}\right).$$

$$T_{\bullet }^{\text{bare}}=$$

$$\times {\left(\frac{M}{3\times 10^9{M}_{\odot }}\right)}^{1/2}{\left(\frac{m}{10^{8}M\odot }\right)}^{-1}{\left(\frac{r_{\text{infl}}}{300\text{ pc}}\right)}^{3/2}\text{ yr},$$

$$T_{\bullet }^{\text{gx}}=$$

$$\times \left(\frac{M}{3\times 10^9{M}_{\odot }}\right){\left(\frac{100{\text{ km s}}^{-1}}{{\sigma }_s}\right)}^3\text{ yr}.$$

$$T_{\bullet }=\text{min}\left(T_{\bullet }^{\text{bare}},T_{\bullet }^{\text{gx}}\right).$$

$$a_h\approx 36\frac{q}{{\left(1+q\right)}^2}\frac{M+m}{3\times 10^9{M}_{\odot }}{\left(\frac{\sigma }{300{\text{ km s}}^{-1}}\right)}^{-2}\text{ pc},$$

$$T_{\text{h,GW}}\approx$$

$$\times$$

$$\frac{a_h}{a_{\text{GW}}}\approx 55{\left(\frac{r_{\text{infl}}}{30\text{ pc}}\right)}^{5/10}{\left(\frac{M+m}{10^8{M}_{\odot }}\right)}^{-5/10}{\left(\frac{4q}{{\left(1+q\right)}^2}\right)}^{4/5},$$

$$T_{\text{GW}}=2.3\times 10^6\left(\frac{a_h}{10^{-3}\text{pc}}\right){\left(\frac{10^8{M}_{\odot }}{M}\right)}^2\left(\frac{10^5{M}_{\odot }}{m}\right)\left(\frac{1}{1+q}\right)\text{ yr}.$$

$$T_{\text{final}}=\{\begin{array}{ll}{T}_{\text{h,GW}},& \text{if}q\ge 10^{-3}\\ T_{\text{GW}},& \text{if}q<10^{-3}.\end{array}$$

$$t_{\text{coal}}=T_{\star }+T_{\bullet }+T_{\text{final}}.$$

$$\frac{dv}{dt}\bigg{|}_{\textrm{DF}}=-\frac{2\pi G^{2}(M+m_{\textrm{\tiny DF}})\rho}{v^{2}}\ln\Lambda_{c},$$

$$m_{\text{DF}}=m{\left(\frac{m+m_{\text{host}}}{m}\right)}^{1-t/{\text{uptau}}_{\text{dyn}}}.$$

$$\frac{da}{dt}\bigg{|}_{\textrm{SS}}=-\frac{G\rho}{\sigma}a^{2}\,H,$$

$$h_c^2=4f^2{\left|\tilde{h}(f)\right|}^2,$$

$$f_{\text{mrg}}=\frac{c^3}{G}\frac{1}{6^{3/2}\pi M_{\text{T}}\left(1+z\right)},$$

$$f(t)=\frac{c^3}{G}\frac{1}{8\pi M_T(1+z){\tau }^{3/8}}\left(1+\left(\frac{11}{32}\eta +\frac{743}{2688}\right){\tau }^{-1/4}\right),$$

$$\tau =\frac{c^3}{G}\frac{\eta t}{5M_T(1+z)}.$$

$$f_{\text{end}}=\frac{G}{c^3}\frac{0.2}{M_T(1+z)}.$$

$$⟨{\rho }^2⟩=\frac{16}{5}{\int }_0^{\infty }\frac{h_c^2}{h_N^2}\frac{1}{f}df.$$

$${\rho }_{\text{GLRT}}^{2,3}$$

$${\rho }_{\text{GLRT}}^{2,4}$$

$$\frac{d^{2}n(z)}{dzd{V}_c}\approx \frac{N(z)}{\Delta z{V}_c},$$

$$\frac{d{N}_{\text{coal}}}{d{t}_{\text{obs}}}={\int }_0^{\infty }\frac{d^2\bar{n}(z)}{dzd{V}_c}\frac{dz}{dt}\frac{d{V}_c}{dz}\frac{dz}{1+z},$$

$$\lambda =\frac{d{N}_{\text{coal}}}{d{t}_{\text{obs}}}{t}_{\text{dur}}.$$

$$W(z)=\frac{dz}{dt}(z)\frac{d{V}_c}{dz}(z)\frac{1}{1+z}.$$

Full text

††pubyear: 2018††pagerange: Probing Massive Black Hole Binary Populations with LISA–A.1

1 Introduction

With the selection of the Laser Interferometer Space Antenna (LISA; Amaro-Seoane et al. 2017) by the European Space Agency (ESA) for its L3 mission, the massive black hole (MBH) community will gain an important tool for understanding the physics and evolutionary history of MBHs. Following a galaxy merger, the two MBHs from the galactic centers eventually form a bound pair and become a binary. This evolving MBH pair radiates gravitational waves at a range of low frequencies $\sim 1\ \text{nHz}-1$ mHz. At the higher end of this range, when the two MBHs are nearer to coalescing, their gravitational wave emission will make them a prime target for the LISA mission. LISA is sensitive to MBH binaries of $\sim 10^{3}-10^{9}M_{\odot}$ (Klein et al., 2016; Amaro-Seoane et al., 2017; Katz & Larson, 2019) at frequencies from $0.1\ \text{mHz}-10$ mHz corresponding to separations of less than $10^{3}$ Schwarzschild radii. Signals from these binaries will allow the scientific community to study the origin and evolution of MBHs over cosmic time (eLISA Consortium et al., 2013; Barausse et al., 2015; Amaro-Seoane et al., 2017). Additionally, MBH binary signals have the potential to reach very high signal-to-noise ratios (SNR, also referred to as $\rho$), unattainable for ground-based gravitational wave detectors (eLISA Consortium et al., 2013; Barausse et al., 2015). This will allow for high precision measurements of cosmological parameters as well as a greater understanding of fundamental physics (Gair et al., 2013; Barausse et al., 2015).

In this paper, we use the Illustris large-scale cosmological simulations (Vogelsberger et al., 2014a, b; Genel et al., 2014; Sijacki et al., 2015) to analyze MBH binary populations, their dynamics, and LISA detection prospects. A few papers have provided similar studies. Blecha et al. (2016) performed an analysis on MBH recoil kicks using Illustris data. Kelley et al. (2017a); Kelley et al. (2017b) and Kelley et al. (2018) performed a similar analysis to ours in relation to Pulsar Timing Array (PTA) predictions. There have been many papers predicting rates for LISA (e.g. Klein et al., 2016; Berti et al., 2016; Salcido et al., 2016; Bonetti et al., 2019), as well as new predictions for the TianQin gravitational wave observatory, which is similar in construction to LISA (Wang et al., 2019). Many rate prediction papers are built from semi-analytic models (SAM). Most rate predictions with large seeds from SAMs, like those in Klein et al. (2016) and Berti et al. (2016), use seeds on the order of $10^{4}M_{\odot}$. These masses are unresolved in the Illustris simulations; therefore, our overall rate estimates will be lower than those predicted in these papers. Similarly, as delays between galaxy mergers and their central MBH mergers have been included in SAMs, the predicted rates have lowered a bit from $\sim 20$ yr*-1* (Arun et al., 2009; Sesana et al., 2011) to $\sim 8$ yr*-1* (Klein et al., 2016; Berti et al., 2016). However, the inclusion of triple MBH interactions has increased the rate to $\sim 20$ yr*-1* (Bonetti et al., 2019). In addition to using MBH binaries from Illustris, we use improved prescriptions for the MBH binary evolutionary process (delay prescriptions) from Dosopoulou & Antonini (2017) and Kelley et al. (2017a); Kelley et al. (2017b).

Black holes with masses of $\sim 10^{6}-10^{10}M_{\odot}$ are usually considered MBHs. These MBHs are believed to exist in the centers of most galaxies of considerable size (Soltan, 1982; Kormendy & Richstone, 1995; Magorrian et al., 1998). This stems from observations of dynamics in the center of other galaxies, as well as our own Milky Way Galaxy. The MBH in the center of the Milky Way Galaxy has been constrained to a mass of $\sim 4.1\times 10^{6}M_{\odot}$ (Boehle et al., 2016). Intermediate mass black holes (IMBH) with masses of $\sim 10^{2}-10^{5}M_{\odot}$ have been theorized and observations of IMBHs have been suggested (e.g. Lin et al., 2018; Bellovary et al., 2019), but remain uncertain. However, electromagnetic observations are beginning to find MBHs in between these ranges ($10^{5}-10^{6}M_{\odot}$) in dwarf galaxies (Reines et al., 2013; Moran et al., 2014; Satyapal et al., 2014; Lemons et al., 2015; Sartori et al., 2015; Pardo et al., 2016; Nguyen et al., 2018, 2019).

Within the $\Lambda$CDM paradigm (e.g. White & Frenk, 1991), galactic halos merge. If conditions are right, the MBHs in their centers can form a binary (Begelman et al., 1980). The pair of MBHs will inspiral via interactions with surrounding gas and stars until they reach close enough separations to emit detectable gravitational waves (e.g. Sesana et al., 2004; Haiman et al., 2009; Sesana, 2010; Roedig et al., 2011; Dosopoulou & Antonini, 2017; Kelley et al., 2017a; Kelley et al., 2017b; Rasskazov & Merritt, 2017). The only way to study MBHs to date is through electromagnetic observations, which lead to a potential observational bias in the MBH population measurements (Shen et al., 2008; McConnell & Ma, 2013; Shankar et al., 2016; Rasskazov & Merritt, 2017). LISA will add a strong and independent method for studying these exotic objects. Additionally, the combined measurement of the luminosity distance to a binary with LISA and the electromagnetic measurement of a binary redshift can be used as a “standard siren” to measure the Hubble parameter (Schutz, 1986; Holz & Hughes, 2005; Abbott et al., 2017). To perform this measurement, the host galaxy of the binary will generally be needed to get the redshift value. Many groups suggest that some periodic active galactic nuclei could be MBH binaries (e.g. Graham et al., 2015; Charisi et al., 2016; Liu et al., 2016). This type of measurement can help identify the host galaxy and, therefore, the redshift of the binary. However, even without an electromagnetic counterpart, LISA measurements can help constrain various cosmological parameters (Petiteau et al., 2011). In addition to cosmological parameters, EM counterparts of MBH binaries can help illuminate a variety of astrophysical processes, including accretion physics and galaxy evolution (Burke-Spolaor, 2013; Bogdanović, 2015).

Understanding formation channels of these large MBHs and how they relate to galaxy formation models is an active area of research. Leading theoretical ideas about MBH formation channels largely differ in their considerations for the mass of MBH seeds at early times in the evolution of the Universe. One such scenario involves the direct collapse of pre-galactic halos with a seed mass on the order of $10^{4}-10^{6}M_{\odot}$ at redshifts of $10-20$ (Loeb & Rasio, 1994; Begelman et al., 2006; Latif et al., 2013; Habouzit et al., 2016; Ardaneh et al., 2018; Dunn et al., 2018). Another scenario involves seeds of $\sim 10^{3}-10^{4}M_{\odot}$ from runaway cluster collapse (Omukai et al., 2008; Devecchi & Volonteri, 2009; Davies et al., 2011; Katz et al., 2015). For smaller seeds, MBH formation channels involve seeds from the collapse of large Population III stars into black holes on the order of $10^{2}M_{\odot}$. This would occur at earlier times in cosmic history at redshifts of $20-50$ (Haiman et al., 2000; Fryer et al., 2001; Heger et al., 2003; Volonteri et al., 2003; Tanaka & Haiman, 2009; Alvarez et al., 2009). Regardless of the formation channel, observations of active galactic nuclei of $\sim 10^{9}M_{\odot}$ in the centres of galaxies at $z\sim 6-7$ indicate MBHs in the universe must have formed quickly on cosmological timescales after the Big Bang, within 1 billion years (Fan et al., 2001a, b, 2006; Mortlock et al., 2011).

It may be that reality is a combination of these three theories. Additionally, different accretion types and spin values can affect the growth rate of MBHs (Plowman et al., 2010; Plowman et al., 2011; Sesana et al., 2011). Due to the possible measurement bias of AGNs with EM observations (Lauer et al., 2007; Schulze & Wisotzki, 2011), these possibly conflicting scenarios will benefit greatly from LISA detections at higher redshifts across the entire mass spectrum.

The coalescence of two MBHs occurs after a long dynamical process forcing the two MBHs to decay from $\sim$kpc separations down to merger (Begelman et al., 1980; Yu, 2002; Merritt & Milosavljević, 2005). Once two galaxies have merged, their central MBHs sink to the center via dynamical friction from interactions with surrounding stars (Chandrasekhar, 1943; Quinlan, 1996; Quinlan & Hernquist, 1997). Once the two MBHs become gravitationally bound, the dynamical friction formalism breaks down, and individual interactions between singular stars and the binary must be considered. These interactions extract angular momentum from the binary, driving them closer to each other (e.g. Merritt, 2013; Vasiliev & Merritt, 2013). This regime is the “stellar hardening” or “loss-cone scattering” regime, which refers to the specific cone in parameter space where stars have to exist in order to extract angular momentum from the binary (Frank & Rees, 1976; Lightman & Shapiro, 1977). Following these interactions, the binary becomes close enough to interact with a circumbinary gas disk if gas is present. Generally, the torque from the gas disk is expected to bring the binary closer to coalescence (Haiman et al., 2009). However, there is growing evidence that specific binary parameters and gas disk properties can lead to the gas disk forcing the binary outwards to larger separations (Muñoz et al., 2019; Moody et al., 2019). After interaction with a gas disk, or if there is little to no gas present, the binary will enter the gravitational wave regime where it will evolve until coalescence. Once the binary enters the gravitational wave regime, its dynamics follow the formalism of Peters & Mathews (1963) at small separations of $\sim 100-1000$ Schwarzschild radii (Kelley et al., 2017a). Several studies have predicted some residual eccentricity when MBH binaries enter the LISA frequency band (Porter & Sesana, 2010; Amaro-Seoane et al., 2010; Mirza et al., 2017; Dosopoulou & Antonini, 2017). However, we will treat all binaries as circular for convenience. This is also conservative as the eccentricity will cause the binaries to merge faster.

We will examine binary lifetime models from Dosopoulou & Antonini (2017) and Kelley et al. (2017a); Kelley et al. (2017b). In section 3.1 we will discuss our models and their specific mathematical approaches to binary lifetime calculations. For this paper, in order to match with cosmological parameters in the Illustris simulation, we assume a WMAP-9 cosmology with $H_{0}=70.4$ km/s/Mpc, $\Omega_{\text{M}}=0.2726$, $\Omega_{\text{b}}=0.0456$, and $\Omega_{\text{vac}}=0.7274$ (Hinshaw et al., 2013).

2 Illustris Simulation

The Illustris cosmological simulations are a suite of simulations evolving gas cells, dark matter (DM), star, and MBH particles from $z=137$ to $z=0$ in a cube of side length $106.5$ comoving Mpc. Illustris is based on the moving, unstructured-mesh hydrodynamic code Arepo (Springel, 2010). Illustris reproduced statistics of large-scale galaxy assembly as well as internal structures of elliptical and spiral galaxies (Vogelsberger et al., 2014b). In particular, we extracted data from Illustris-1, the highest resolution simulation in the suite, with $1820^{3}$ gas cells and DM particles. DM particles have a mass resolution and typical gravitational-softening length of $\sim 6.3\times 10^{6}M_{\odot}$ and 1.5 kpc, respectively. At redshift zero, there are over $3\times 10^{8}$ star particles, which have a resolution of $\sim 1.3\times 10^{6}M_{\odot}$ and a typical softening length of 700 pc.

MBHs are implemented in Illustris as massive sink particles. When halos attain a total mass of $7.1\times 10^{10}M_{\odot}$, they are seeded with an MBH of mass $1.42\times 10^{5}M_{\odot}$ if it does not already have an MBH in it (Sijacki et al., 2015); the highest density gas cell in the halo is converted to the MBH particle. At this point, the initial dynamical MBH mass will be the same as the gas cell in its previous state. However, the MBH particle is assigned an internal mass of the seed mass, which is tracked from this point (Vogelsberger et al., 2013).

MBH particles grow by Eddington-limited, Bondi-Hoyle accretion from their parent gas cells initially, and then from its nearby gas cells after the dynamical mass of the particle becomes equal to its internal mass. From this point, the dynamical mass and internal MBH mass increase in tandem. The similarity in mass between the MBH particles, gas cells, and star/DM particles would cause the MBH particle to scatter around halos in an unphysical manner without settling down in the center of the halo. Therefore, MBH particles are repositioned to the potential minimum of the host halo at every time step (we refer to this as the “repositioning algorithm”).

For more overview on the Illustris simulations, see Vogelsberger et al. (2013) and Torrey et al. (2014). See Vogelsberger et al. (2014a), Genel et al. (2014), and Sijacki et al. (2015) for detailed simulation results and comparisons of the simulations to observations. The initial data sets used in this study were all obtained at www.illustris-project.org (Nelson et al., 2015). Using this data, we perform more post-processing to create the final datasets used in our analysis. We will further explain this in the following sections.

2.1 Massive Black Hole Merger Population

Two MBH particles are merged in the Illustris simulation when they come within a smoothing length of each other ($\sim$kpc). Since these mergers occur at larger scales, we treat this merger event as the formation of the binary, from which we evolve the binary to coalescence with sub-grid models. Following Kelley et al. (2017a), we adopt the term “merger” to indicate this simulation-only process: the combination of two MBH simulation particles into one, indicating the formation of a binary on $\sim$kpc scales. We will refer to the final combination of two realistic MBHs into one as the MBH binary “coalescence.”

Over the course of the simulation, detailed MBH and host-galaxy information is saved in a series of 135 snapshots. Higher time-resolution data was saved for each merger, including the time of the merger and the constituent MBH masses (Blecha et al. 2016, Kelley et al. 2017a). We extract detailed properties for all MBHs in the simulation at each snapshot, as well as relevant, global properties for each of their host galaxies from the Illustris “Group Catalogs” (Nelson et al., 2015). In addition to the higher time-resolution merger data set, we extract information about host galaxies related to the mergers. For each merger, we locate the host galaxies of the constituent MBHs at the snapshot immediately preceding the merger, as well as the host galaxy of the remnant MBH at the snapshot immediately following the merger. For these galaxies, we not only attained global information, we also gather information about their specific distribution of gas, stellar, and DM constituents. The last data needed for our analysis is the Sublink merger trees (Rodriguez-Gomez et al., 2015) to follow galaxies from snapshot to snapshot.

Throughout the cosmic history within Illustris, there are 23,708 MBH merger events. However, a fraction of these mergers are artificial. First, the Friends-of-Friends (FOF) halo finder will occasionally associate two halos as one. When this occurs, the aforementioned repositioning algorithm will force the two MBHs in the centers of each galactic halo to the new potential minimum determined during this misstep by the FOF finder. This causes the two MBHs to merge. After this “fly-by” encounter, the two galaxies may separate into two distinct halos as seen by the FOF finder. When this occurs, there will be one galaxy without a central MBH. At this point, a new MBH is seeded in this galaxy causing future artificial mergers to inflate the merger catalog. Similarly, the FOF finder may identify a transient matter overdensity, subsequently seeding a low-mass MBH into the overdensity. This newly seeded MBH is then quickly merged into the MBH in the nearest massive halo due to the repositioning algorithm, once again adding unphysical mergers.

Previously, this was dealt with in Kelley et al. (2017a) and Blecha et al. (2016) using a cut based on mass–only $M_{\bullet}>10^{6}M_{\odot}$ are kept–to exclude low mass MBHs that are overwhelmingly the MBHs involved in these numerical issues. As the authors state, this cut had minimal effect on their predictions for the PTA background, which is dominated by high-mass MBH binaries $M\gtrsim 10^{8}M_{\odot}$). In a study about LISA, these near-seed mass MBHs play a very important role. Therefore, we designed a post-processing method to avoid removing these small MBHs in the most robust way possible.

To handle these issues we start by requiring that all merger constituent MBHs must exist for at least one snapshot prior to the merger, which is always true for MBHs above $10^{6}M_{\odot}$. This removes the MBHs seeded when the FOF finder identifies an overdensity and seeds an unphysical MBH. This is effective in removing these MBHs because the time between the seeding and the merging of the unphysical MBH is less than the duration of one snapshot. We then focus on identifying galaxies that have had their central MBHs removed by the ambiguities related to the FOF finder and the repositioning algorithm. We track the evolution of the galaxy devoid of an MBH as it continues on after its fly-by encounter where it lost its MBH. If this galaxy seeds a new MBH before it merges with another galaxy, we remove this MBH from our catalog, as well as any of its subsequent mergers.

In addition to filtering the secondary mergers, we analyzed the effects of the premature mergers created by the FOF association of two separate halos. In this process, the removal of the MBH from its galaxy occurs earlier in cosmic time than the actual galactic mergers. This causes the MBH merger to occur at slightly earlier cosmic times. This effect would increase our rate predictions: mergers at earlier cosmic times inflate the number of mergers because the volume of the observer’s past light-cone is larger at higher redshifts. However, we believe this effect to be small for two reasons: the binaries involved in this scenario tend to have largely unequal mass ratios, making them increasingly difficult to detect with LISA; and the delay between the unphysical MBH merger and the subsequent galaxy merger is on the order of $10^{8}$ yrs, which represents a small fraction of the lifetime of the universe as well as a small fraction of the average binary inspiral times predicted from our two binary inspiral models.

We will now discuss the post-processing performed on the merger host galaxies as well as the cuts made based on those galaxies.

2.2 Host Galaxy Information

To process the host galaxy data, we follow the process of Blecha et al. (2016) and Kelley et al. (2017a). We need density profiles and velocity dispersions of remnant host galaxies for input into our evolution timescale models. We use the profiles at resolvable scales to extrapolate inward to the centers of the galaxies to infer properties at unresolved scales. Therefore, we confirm each remnant host galaxy is sufficiently resolved for these calculations by requiring the galaxy contain at least 80 DM particles, 80 gas cells, and 300 star particles.

From these remaining galaxies, we construct spherically averaged, radial density profiles for stars, gas, and DM. We calculate these profiles based on the innermost shells of particles (cells for the gas) surrounding the galactic center. We assume the profile represented by the innermost shells of particles/cells extends inwards to the core of the galaxy. We require at least four particles/cells in each radial bin. The profile is then formed from the innermost 8 bins satisfying the four particle/cell minimum requirement. Binaries were excluded from the final catalogs when fits could not be constructed under these requirements. A graphical representation of this process can be seen in Figure 1 of Kelley et al. (2017a). For our evolution prescription, we constrain the density profile index to be in between 0.5 and 2.5. The distribution function for the models tested becomes unphysical below an index of 0.5. The upper end of 2.5 is determined based on observed stellar cusps of giant elliptical galaxies. After all of our cuts, we are left with 17,535 of the original 23,708 mergers. We compare this number to 9,270, which is the amount of mergers remaining after the cuts applied in Kelley et al. (2017a). Therefore, we analyze 8,265 more mergers in this work. The mergers that remain form our final merger catalog. Figure 1 compares the main binary parameters resulting from our extraction method to the flat mass cut of $10^{6}M_{\odot}$. Figure 2 shows the main properties of mergers in this final catalog.

3 Methods

3.1 Binary Lifetime Models

For the following models, our goal is to calculate the evolutionary timescale from $\sim$kpc scales to coalescence, which is usually of order $\sim$Gyr. We use this timescale to find the coalescence time by adding the evolution time to the time of binary formation (particle merger in Illustris).

The first model we examine has no evolution of the binaries from their particle mergers in the Illustris simulation. We refer to this model as “ND” for no delays. In other words, we consider the formation time to equal the coalescence time. This model is our baseline model against which we compare our more detailed models for binary MBH coalescence timescales: it represents the exact prediction from the simulation if the merger process is not modeled below $\sim$kpc scales, which is a common assumption in rate prediction papers.

In addition to our ND model, we will examine a subset of our “no delays” model requiring masses to be greater than $10^{6}M_{\odot}$. This allows us to test the difference in our extraction process by comparing the new data set to the old extraction data set similar to the one used in Kelley et al. (2017a). We will refer to this model as “ND-6.”

In the following, we will describe detailed models that have been constructed and analyzed in previous papers. Therefore, we give a quick overview of each model. For more information on the DA17 model (section 3.1.1), please see Dosopoulou & Antonini (2017). For the K17 model (section 3.1.2), please see Kelley et al. (2017a); Kelley et al. (2017b).

3.1.1 DA17 Model

The equations shown below are taken directly from Dosopoulou & Antonini (2017). In what follows, primary, $M$, (secondary, $m$) will refer to the larger (smaller) MBH. For this model, we are assuming the mergers are gas-poor.

The initial timescale in this model is the large-scale orbital decay from $\sim$kpc scales to a separation equal to the influence radius, $r_{\text{infl}}$, of the primary MBH. $r_{\text{infl}}$ is a shorter length scale than the resolution in the Illustris simulation. Therefore, we use an approximation from Merritt et al. (2009) given by,

[TABLE]

where $\sigma$ is the three dimensional stellar velocity dispersion of the primary galaxy. Modeling the primary host galaxy as a singular isothermal sphere, Binney & Tremaine (1987) show this decay timescale is given by,

[TABLE]

where $\ln{\Lambda}$ is the Coulomb logarithm and $R_{e}$ is the effective radius of the primary galaxy. However, Dosopoulou & Antonini (2017) modify this formalism to include the effect of the secondary MBH remaining embedded in a core of stars from the secondary galaxy. This added mass causes the system to sink faster towards the primary galaxy’s center, therefore, taking less time than is predicted by Equation 2. Assuming the mass of stars bound to the secondary to be a constant proportionality of $10^{3}m$ (Merritt & Ferrarese, 2001), this large-scale decay timescale becomes,

[TABLE]

where $\Lambda^{\prime}=2^{3/2}\sigma/\sigma_{s}$ with $\sigma_{s}$ representing the stellar velocity dispersion of the secondary galaxy. This equation, however, does not include tidal stripping of stars from the secondary galaxy by the primary galaxy. Including this effect, the large-scale decay timescale is given by,

[TABLE]

Using Equations 3 and 4, we can approximate the large-scale decay timescale as

[TABLE]

During the large-scale decay, the stellar velocity distribution is treated as Maxwellian. Once the secondary MBH reaches $r_{\text{infl}}$ of the larger MBH, this assumption no longer holds because the potential is dominated by the central MBH. Therefore, the stellar velocity distribution is treated according to Equation 20 in Dosopoulou & Antonini (2017). We refer to this next regime as the “dynamical friction regime” to match the conventions of the original paper. The timescale for the binary to decay to a shorter separation $r=\chi r_{\text{infl}}$ ($\chi<<1$) is given by,

[TABLE]

where $\gamma$ is the power law exponent in the stellar density profile, $\rho(r)=\rho_{0}\left(r/r_{\text{infl}}\right)^{-\gamma}$; $\alpha$, $\beta$, and $\delta$ are calculated from Equations 17-19 in Dosopoulou & Antonini (2017).111In these equations, we assume a circular orbit setting $\xi=1$ (see the paper for more details). However, the dynamical friction decay timescale is not greatly affected by the orbital eccentricity. If we include the stars bound to the secondary, as in Equation 4, this timescale becomes,

[TABLE]

Similar to Equation 5, we find $T_{\bullet}$ with,

[TABLE]

We use Equation 8 to evolve the binary down to the hardening radius, $a_{h}$, given by (Merritt, 2013),

[TABLE]

where $q$ is the mass ratio ($q\leq 1$). Therefore, we set $\chi=a_{h}/r_{\text{infl}}$.

The final phase in the DA17 model is the “hardening phase”, as it includes the effect of gravitational radiation. With dry mergers, remnant galaxies after a merger are expected to be triaxial. In this configuration, an efficient hardening of the binary is exposed to a full and consistently refilling loss-cone (Khan et al. 2011, Vasiliev et al. 2014). From $a_{h}$ until coalescence, including the gravitational wave regime, the timescale is given by (Vasiliev et al., 2015),

[TABLE]

where $\phi=0.4$ and $\psi=0.3$ are triaxial parameters estimated from Monte Carlo simulations in Vasiliev et al. (2015). In Equation 10, we left out the eccentricity factor as it is unity because we are assuming circularity. Gravitational radiation takes over at $a_{\text{GW}}$, determined from the ratio (Vasiliev et al., 2015),

[TABLE]

where we have once again left out the eccentricity factor as it is equal to unity. For $q\approx 10^{-3}$, $a_{h}$ is less than $a_{\text{GW}}$. For these binaries, we use the inspiral time for a circular binary due to only gravitational radiation according to (Peters & Mathews, 1963),

[TABLE]

Therefore, our final timescale, $T_{\text{final}}$ is given by,

[TABLE]

The DA17 model’s final coalescence timescale, $t_{\text{coal}}$, is therefore given by

[TABLE]

3.1.2 K17 Model

This section introduces the K17 model, described in detail in Kelley et al. (2017a) and Kelley et al. (2017b). Binaries are numerically integrated from their formation at large-separations until their eventual coalescence. Dynamical friction is implemented following Chandrasekhar (1943), where the deceleration is given as,

[TABLE]

where the relative velocity is taken to be the maximum of the orbital velocity and stellar velocity dispersion, i.e. $v=\max\left(v_{\textrm{orb}},\sigma\right)$, $\rho$ is the total mass-density, and the Coulomb logarithm is set to $\ln\Lambda_{c}=15$. The effective mass of the secondary, $m_{\textrm{\tiny DF}}$, assumes that the mass of the secondary host galaxy is stripped over the course of a dynamical time, i.e.,

[TABLE]

The initial host galaxy mass is measured from Illustris in the snapshot preceding the merger event. Once the binary shrinks below the “loss-cone radius” (Begelman et al., 1980), the hardening rate is calculated following the stellar-scattering prescription from Sesana et al. (2006),

[TABLE]

where $a$ is the semi-major axis of the binary which is being integrated, and $H$ is a dimensionless coefficient calculated from numerical scattering experiments.

The accretion rate calculated in Illustris provides an estimate of the presence of circumbinary gas. To model the energy extraction from this material, we assume the gas settles into a geometrically-thin alpha-disk (Shakura & Sunyaev, 1973) with different regions corresponding to the dominant components of pressure (radiation vs. thermal) and opacity (Thomson vs. free-free), following Shapiro & Teukolsky (1986). The hardening rates in each regime are calculated in Haiman et al. (2009), as a function of disk surface-density and binary mass-ratio. Numerical simulations have found these analytic prescriptions to be quite accurate over the parameter ranges studied (Tang et al., 2017; Fontecilla et al., 2019). We assume these disks extend out to the radius at which they become Toomre unstable. GW energy extraction is implemented at all radii, following Peters & Mathews (1963).

3.2 Determining Detectability

3.2.1 Characteristic Strain

MBH binaries provide a variety of signals measurable by LISA since their chirp evolution in the frequency domain occurs near the low-frequency band edge of the LISA sensitivity curve. Binaries with $\sim 10^{5}-10^{7}M_{\odot}$ total mass will provide a measurable inspiral, merger, and ringdown leading to very loud signals even out to the cosmic horizon (Amaro-Seoane et al., 2017).

The binary inspiral is the initial stage of binary black hole coalescence when the two MBHs orbit one-another at separations greater than the innermost stable circular orbit (ISCO; $R=6GM/c^{2}$). At these separations, the orbit is usually treated with a post-Newtonian formalism.

The merger stage follows the binary inspiral with a highly non-linear relativistic process. This process continues until the MBHs have contacted each other to form a single event horizon, leading to ringdown. The dominant mode of the ringdown spectrum is expected to be the $l=m=2$ quasinormal mode. Deviations from general relativity can be measured if LISA can detect subdominant modes in the ringdown spectrum. This process is referred to as the so-called “black hole spectroscopy” (Berti et al., 2006; Berti et al., 2016; Baibhav et al., 2018; Baibhav & Berti, 2019).

We use the characteristic strain, $h_{c}$, to model the binary signal which accounts for the time the binary spends in each frequency bin (Finn & Thorne, 2000). The characteristic strain is given by (Moore et al., 2015),

[TABLE]

where $\tilde{h}(f)$ represents the Fourier transform of a time domain signal. To find $\tilde{h}(f)$, we use the phenomenological waveform PhenomD (Husa et al., 2016; Khan et al., 2016). PhenomD is based on fitting analytical templates to numerical relativity waveforms. For a detailed description of its constructions, please see Husa et al. (2016) and Khan et al. (2016). Here, we focus on how the waveform is determined based on the parameters of the MBHs in our population.

To generate the waveforms, we use the gwsnrcalc Python package from the BOWIE analysis tool (Katz & Larson, 2019). gwsnrcalc takes as inputs the masses of the MBHs, $M$, $m$; the dimensionless spin of each MBH $a_{1}$, $a_{2}$; the redshift of the binary, $z$; and the start and end times of the binary’s orbit, in relation to the merger of the binary, $t_{\text{st}}$ and $t_{\text{end}}$.

The dimensionless spin of each MBH is $a_{i}=J_{i}/m_{i}^{2}$, where $J$ is the magnitude of the spin angular momentum. $a$ ranges from -1.0 (anti-aligned to the orbital angular momentum) to 1.0 (aligned to the orbital angular momentum). For convenient use of PhenomD, we treat the spins as aligned. Measurements of MBH spins have shown spins near maximal (Miller, 2007; Reynolds, 2013). For this reason, we choose to model spins of $a_{1}=a_{2}=a=0.8$. As the spin magnitude is raised, the waveform will gain more signal. For near-equal mass systems, which represent a majority of systems in our catalog, the difference in the spin does not change the signal significantly. For systems of mass ratio farther from unity, the spin can have a significant impact on their detectability because the signal peak can increase by an order of magnitude from the spin-down ($a=-1$) to the spin-up case ($a=1$). Therefore, applying this spin configuration ($a=0.8$) represents the optimistic case for these systems. The choice to use the same spin for both MBHs is made because PhenomD was calibrated in mostly equal-spin configurations. Within its calibration range, PhenomD performs accurately matching waveforms to better than $\sim 1\%$ error. Outside of its calibration range, it produces physically reasonable results, indicating it can be useful for basic studies (Khan et al., 2016). See Figure 3 for examples of characteristic strain curves.

3.2.2 Start and End Times

The start times, $t_{\text{st}}$, and end times, $t_{\text{end}}$, both represent the time until merger for a specific binary at which LISA begins and ends its observation of the binary’s signal. These times are effectively a map to the frequency bounds of the PhenomD waveform model. These times will be on the order of years. Since the merger and ringdown timescale is on the order of minutes to hours for MBHs, we do not include this timescale in $t_{\text{st}}$ and $t_{\text{end}}$. The main reason for this construction is it allows us to refrain from assuming a specific observation time for each binary: it allows us to test, within our Monte Carlo sample, binaries that merge at early and late times in the LISA observation window, as well as binaries that merge after the LISA observation window. Early merging binaries will have less time where their inspiral signal can be observed, compared to later mergers. Similarly, LISA will not be able to detect the merger or ringdown for binaries that merge after the LISA observing window, leading to inspiral-only signals (if the inspiral is detectable over the observation time). Figure 4 displays a diagram showing how $t_{\text{st}}$, $t_{\text{end}}$, and $T_{\text{obs}}$ are related for various sources. Please see Section A in the Appendix for our analysis related to inspiral-only signals.

As a point of reference, the merger frequency separating inspiral from merger, is given by

[TABLE]

where $M_{\text{T}}=M+m$ is the total mass of the binary in its source frame. The $(1+z)$ term redshifts this mass to the detector frame.

The frequency at a time before merger to 1PN order is given by (Blanchet, 2014),

[TABLE]

where,

[TABLE]

$t$ is given in the detector frame. The start frequency of the waveform, $f_{\text{st}}$, is, therefore, $f(t_{\text{st}})$. For signals that exhibit inspiral, merger, and ringdown, $t_{\text{end}}$ is zero. In this case the end frequency, $f_{\text{end}}$, is the highest frequency used in the PhenomD model, representing the end of the ringdown, given by (Khan et al., 2016),

[TABLE]

If $t_{\text{end}}$ is not zero, indicating that the source is only detected in the inspiral stage, $f_{\text{end}}$ is $f(t_{\text{end}})$.

3.2.3 LISA Sensitivity

The LISA sensitivity configuration used is from the LISA Mission Proposal (Amaro-Seoane et al., 2017). We refer to this sensitivity as “PL” in Katz & Larson (2019). This sensitivity is based on a 3-arm triangular configuration with 2.5 million km armlengths, 30 cm diameter telescopes, and 2 W end-of-life laser power. In Section A of the Appendix, we perform an analysis comparing PL to an older iteration of the LISA configuration (Larson et al., 2000) to show how the LISA configuration changes will affect various aspects of LISA MBH analysis.

We show our sensitivity curve in Figure 3 in terms of the sky-averaged characteristic strain, $h_{N}$. The sky-averaging factor is 3/20 (Robson et al., 2019). Sensitivity curves are generally presented in terms of the power spectral density (PSD) of the noise, $S_{N}$. To convert from $S_{N}$ to $h_{N}$, we use $h_{N}=\sqrt{fS_{N}}$ (Moore et al., 2015).

We also include the effect of the Galactic background noise in addition to the instrumental noise. We use the analytical approximation of Hiscock et al. (2000) to the Galactic background noise suggested in Bender & Hils (1997). This is shown in Figure 3. Compared to recent predictions from Robson & Cornish (2017), this is a conservative estimate of this noise contribution. The contribution of this background can be decreased with proper global fitting methods and a longer observation window (Robson & Cornish, 2017).

3.2.4 Signal-to-Noise Ratio

We use the SNR to determine the detectability of the sources in our catalog. The SNR is estimated by integrating the ratio of the signal to the noise in the frequency domain. The sky, orientation, and polarization averaged SNR is given by (Robson et al., 2019),

[TABLE]

The SNR is then multiplied by a factor of $\sqrt{2}$ because we consider a 2 channel interferometer.

An additional question we analyze is how many sources exist in our models with a high enough SNR to perform black hole spectroscopy. To do this, we use the General Likelihood Ratio Test (GLRT) formalism suggested in Berti et al. (2016). Using the GLRT, the SNR of the $l=m=2$ ringdown mode is used as a proxy for determining the detectability of the $l=m=3$ or $l=m=4$ modes. Sources can be spectroscopically measured if $\rho_{\text{l=m=2}}>\rho_{\text{GLRT}}\equiv\text{min}\left(\rho_{\text{GLRT}}^{2,3},\rho_{\text{GLRT}}^{2,4}\right)$, where $\rho_{\text{GLRT}}$ for each mode is is given by (Berti et al., 2016),

[TABLE]

3.3 Monte Carlo Analysis

We use a Monte Carlo analysis technique based on Poisson statistics to characterize the range of possibilities resulting from the Illustris output. We do this for multiple reasons. The primary reason is the Illustris output is one iteration of the evolution of a fractional volume within the Universe. We want to understand how the detection rate and source characteristics will vary with the Illustris output as the statistical backdrop. Additionally, we wanted to create a catalog generator for LISA MBH binary signals, which requires a Monte Carlo draw of a new sample each time. As we will discuss in section 3.3.2, the Monte Carlo sampling allows us to refrain from assuming an observable duration of the waveform for each of the binaries. Most detection rate predictions assume an observable time for each MBH of 1 year. This does not account for binaries that will merge before one year of LISA observation. It also does not include the longer measurement of an inspiral signal if the binary signal is observable at times longer than 1 year before merger. In other words, our method provides a more realistic basis for assessing binary detectability during the LISA mission.

3.3.1 Merger Rate Prediction

The first parameter in the Monte Carlo sampling process is the coalescence rate of MBH binaries in the Illustris simulation. We calculate this parameter one time for each evolutionary prescription. At this stage of the sampling process, we are not considering detectability; we consider any coalescence that occurs prior to $z=0$. For the ND and ND-6 models, all binaries in our sample coalesce before $z=0$ because these two models assume no inspiral time for all MBH binaries. For the binary inspiral models of DA17 and K17, we find that 84% and 66% of binaries coalesced before $z=0$, respectively. See section 4.1 and Figure 6 for more information on how the inspiral models affected the population of binaries coalescing before $z=0$. We determine the number of coalescences, $N$, in a given redshift interval $z+\Delta z$. We then compute the number of coalescence events across redshift intervals per comoving volume element,

[TABLE]

where $V_{c}$, the term on the right hand side, is the comoving volume of the Illustris simulation, (106.5 Mpc)3. We can then calculate the number of coalescences per observing time interval given by,

[TABLE]

where the $1+z$ redshifts the infinitesimal time element in $dz/dt$ to the observer frame time interval. When we refer to the “integral rate calculation,” we are referring to Equation 27. For our Monte Carlo catalogs, this quantity becomes our input into our Poisson rate calculator.

3.3.2 Poisson Sampling

The two parameters needed to perform the desired Poisson sampling is rate of coalescences determined from Equation 27 and the duration for which we want to draw potential sources, $t_{\text{dur}}$. If we only wanted to draw sources for the observation window, we would set $t_{\text{dur}}=T_{\text{obs}}$. However, this would only focus on sources coalescing within the observation window. We also want to test for inspiraling sources that would coalesce some time after the LISA observing window. Therefore, we choose $t_{\text{dur}}>T_{\text{obs}}$. We tested a variety of values for $t_{\text{dur}}$. We found $t_{\text{dur}}=10^{2}$ yrs encompassed all of the observable systems in our catalogs, while maintaining computational efficiency.

Our final Poisson parameter, $\lambda$, is given by,

[TABLE]

In other words, this is the expected number of coalescence events over $10^{2}$ years. For each catalog, we draw the number of sources occurring within our 100 year window from this Poisson distribution.

3.3.3 Event Times

When we have the number of sources drawn, we assign each a random coalescence time, $t_{\text{ev}}$, between zero and $10^{2}$ yrs. At the timescales we are considering ($\sim 10^{2}$ yrs), the distribution of sources over time will not be affected by the evolution of the Universe. $t_{\text{ev}}=0$ indicates an event occurring at the moment LISA begins observations. $t_{\text{ev}}=T_{\text{obs}}$ represents an event occurring at the moment the LISA observation window ends.

When considering the waveform described in section 3.2.1, we use $t_{\text{ev}}$ and $T_{\text{obs}}$ to determine the start ($t_{\text{st}}$) and end ($t_{\text{end}}$) times related to waveform creation. $t_{\text{st}}=t_{\text{ev}}$ because the event time represents the time before merger at the start of LISA observation. Therefore, $t_{\text{end}}=t_{\text{ev}}-T_{\text{obs}}$ if $t_{\text{ev}}>T_{\text{obs}}$. If $t_{\text{ev}}\leq T_{\text{obs}}$, $t_{\text{end}}=0$ (see Figure 4).

3.3.4 Resampling Binary Parameters

After sampling the number of binaries and the event times for each event, we need to sample binary parameters of $M$, $m$, and $z$. To do this, we use kernel density estimation methods. However, there is a key distinction that needs to be made for sampling these binary parameters: we must incorporate the volume and time redshifting factors implicit in the expansion of the universe as weights, $w$, in the density estimation. If you assume an infinitesimal redshift bin width for Equation 27 ($\Delta z\rightarrow dz$), this weighting factor as a function of redshift is given by,

[TABLE]

The weight applied to the $i$th binary is then $w_{i}=W_{i}/\sum_{i=1}^{N}W_{i}$. We also include the covariance across these parameters in the KDE, so as to sample the population accurately.

4 Results

4.1 Binary Lifetimes

We first test and compare our evolutionary prescriptions to understand how the initial population of binaries will change when evolved to coalescence with different sub-grid models. In the ND and ND-6 models, all binaries are considered to be coalesced prior to $z=0$. By modeling the sub-grid physics as in the DA17 and K17 models, some binaries will no longer merge before $z=0$ and will, therefore, deflate the merger rate. Additionally, the mass distributions of the coalesced binaries may change because the prescriptions have different dependencies on the masses. Figure 5 shows the evolutionary timescales calculated for all binaries binned by total mass and mass ratio. The DA17 model generally results in a more peaked distribution, while K17 shows a flatter profile across all plots. However, at higher total masses, the two prescriptions become very similar in their predictions. Figure 6 shows the effect of these evolutionary timescales on the coalescence fractions of our population. This figure illustrates the global differences, in terms of binary parameters, between the two prescriptions. DA17 favors near-equal mass and low total mass systems, while K17 favors near-equal mass systems with total masses towards the higher end. With the K17 model, larger masses are favored because they are embedded in higher density, more centrally concentrated stellar cores. On the other hand, the specific dynamical friction prescription used in the DA17 model (see Equations 6, 7, and 8) causes high-mass systems to exist for longer times in the dynamical friction stage between the influence radius of the larger MBH (Equation 1) and the hardening radius (Equation 9) (please see the original papers for more details). It is also clear the overall coalescence fractions are higher with DA17 than with K17. The overall coalescence fraction for DA17 was 84%. With the K17 model, only 66% of all binaries coalesced before $z=0$.

4.2 Rate Predictions

We calculate merger rates with two methods: integrating the redshift distributions with Equation 27 and Monte Carlo sampling (section 3.3). The results are similar as expected; however, the Monte Carlo aspect allowed for more freedom in terms of not setting specific values for observation duration and starting times before merger. It also ensures we are examining different realizations of the merging MBH population. Table 1 shows integral merger rate results for each prescription. It also shows detection rates for each stage of binary black hole coalescence as well as the rate of sources where spectroscopic measurements (Equation 24) are possible. Similar results for the Monte Carlo method are shown in Table 3. The redshift distributions used in these calculations are shown in Figure 7. For the main integral merger rate and detection rate results, the rate is quoted as per year. For the following integral rate calculations, we estimate our standard deviation in our predictions to be less than 0.01 (based on our chosen redshift bin width), according to Equation 17 in Salcido et al. (2016).

First, comparing ND with ND-6 (“ND-6” represents the subset of binaries in “ND” which have each constituent mass above $10^{6}M_{\odot}$), we see the advanced extraction was important for LISA-related analysis since the ND merger rate was almost two times the merger rate of ND-6. This holds true for the detection rates as well. An interesting statistic here is the inspiral detection rate, 0.23 yr*-1* with ND-6 versus 0.45 yr*-1* with ND, as this demonstrates the importance of the lower-mass systems retained by our advanced extraction technique.

Comparing all four models, we see some interesting results. The hierarchy predicted in section 3.1 is apparent. The DA17 rates resemble much more strongly the ND rates compared to K17 due to the inability of K17 to coalesce the large number of low mass systems. Interestingly, the K17 rates strongly resemble the ND-6 rates, indicating the loss of the low mass systems, as well as the general loss in coalescing systems due to adding sub-grid modelling, caused the rate to decrease to similar levels as without the low mass systems entirely. These aspects can be seen clearly in Figure 7. After calculating the SNR and making our cut at $\rho=8$, we see that the DA17 curve tracks the ND curve, while the K17 curve mirrors the ND-6 curve.

When analyzing the mass distributions, we find all models produce similar results. Mass distributions from our Monte Carlo analysis can be seen in Figure 8. The black solid line shows the limit imposed by the Illustris simulation seed mass at $\sim 10^{5}M_{\odot}$. Binaries with total mass and mass ratio values above this line cannot exist. Similarly, the black dashed line shows the limit imposed by the $10^{6}M_{\odot}$ cut. This once again highlights the effect of the more advanced extraction. ND, DA17, and K17 exhibit roughly the same structure. Specifically, their mean values and higher order moments about the mean values are within a small percentage of each other; however, DA17 has a slightly smaller kurtosis in the mass ratio. This general similarity does indicate the K17 model has a relatively flat effect across the parameter space, where sources are detectable by LISA, suppressing each mass and mass ratio regime in an equivalent manner.

5 Discussion

The merger rates and LISA detection rates of MBH binaries from the Illustris simulation are low compared to the majority of the literature on this subject. Table 2 summarizes predicted rates in the literature, including whether their base population is from a semi-analytical model (SAM) or a hydrodynamic simulation. Additionally, the type of delay prescription employed is mentioned. For this summary, we focused on predictions related to high-mass seeds since this seeding prescription better matches our setup in this paper. Even with our ND model, which represents the rates directly from the simulation without any sub-grid modeling, our predicted rate is low. LISA will be sensitive to binaries $\leq 10^{5}M_{\odot}$. Therefore, any study that cannot resolve the galaxies and halos where these low masses evolve (see section D in Klein et al., 2016) will under-predict the total detection rate. Examining Figure 8, we see LISA primarily observes binaries of $1\leq q\leq 10^{2}$. Therefore, if we consider our population to be complete when we can resolve all constituent MBHs at mass ratios up to $10^{2}$, then our study is only complete above $10^{7}M_{\odot}$ leading to a further underestimation of the overall detection rate. Another way to consider this effect is by understanding how the seeding prescription will affect the base population. Without seeds below $10^{5}M_{\odot}$, the detection rate for masses above $10^{5}M_{\odot}$ will also be deflated: seeds below this mass may grow to a mass above $10^{5}M_{\odot}$ before undergoing a merger. This aspect significantly affects our merger rate prediction as it removes MBHs from our population at the most common and detectable masses in our sample.

The literature references for SAMs shown in Table 2 employ seeds of $\sim 10^{4}M_{\odot}$. Within the subset of these SAMs predictions, the higher rate predictions occur with no binary inspiral model employed (similar to our ND model). However, early SAMs predictions did analyze various delay prescriptions, but for smaller seeds ($\sim 150M_{\odot}$) (e.g. Volonteri et al., 2003; Sesana et al., 2004; Sesana et al., 2005). Historically, SAMs have predicted rates about an order of magnitude larger than our findings; however, our results are within an order of magnitude of those predicted in Salcido et al. (2016) for the EAGLE simulations (McAlpine et al., 2016), indicating similarity between two hydrodynamic-based populations. While hydrodynamic simulations strive to produce populations from simple physics, they are restricted to a much more limited parameter space than SAMs (discussed further below). This means that rate predictions from SAMs and hydrodynamic-based models fundamentally differ.

One clear difference is the ability of SAMs to explore an arbitrarily large range of masses, while hydrodynamic simulations are resolution limited, typically to above $\sim 10^{5}M_{\odot}$. If we consider the increased prevalence of MBHs as we move towards smaller masses, we would expect that models with access to binaries below $10^{5}M_{\odot}$ would greatly inflate the merger rate in a similar fashion to what we have seen when including the lower masses from Illustris through our advanced extraction. Along the same lines, Illustris, and simulations like it, do not access the scales needed to examine populations of dwarf galaxies. As evidence mounts that dwarf galaxies house MBHs in their center in both observations (Reines et al., 2013; Moran et al., 2014; Satyapal et al., 2014; Lemons et al., 2015; Sartori et al., 2015; Pardo et al., 2016; Nguyen et al., 2018, 2019) and simulations (Volonteri et al., 2008; van Wassenhove et al., 2010; Bellovary et al., 2019), we must improve models since these dwarf galaxy sources are enitrely missing from our analysis. Due to the prevalence of dwarf galaxies as well as the general understanding that dwarf galaxies consistently merge into larger “host” galaxies over time, missing the dwarf galaxy MBHs could deflate our rate calculations significantly. This is especially true because these dwarf galaxies will house smaller MBHs that will produce strong signals in the LISA frequency band. Similarly, a number of potential systems may still be missed as we approach the resolution limit of $\sim 10^{5}M_{\odot}$.

Similar to this issue with dwarf galaxies, the MBH seeding mechanism plays a large role in rate calculations. The seeding mechanism in Illustris (and EAGLE) is ad-hoc. The seeding model chosen for Illustris produces seeds at later times than those seen in other simulations like Tremmel et al. (2018). Additionally, the Illustris seeding prescriptions will produce seeds at much later times when compared to SAMs. Later seeding means a smaller volume accessible to LISA observations. Therefore, merger rates from SAMs may also be intrinsically higher at higher redshits.

The sub-grid models also affect our rates given that they cause coalescences to occur after $z=0$. The particular model chosen has the potential to deflate the rate significantly. As previously mentioned, SAM predictions decreased as more “delay” models were incorporated. In Tremmel et al. (2018), their simulation works to numerically resolve smaller scales in order to track MBH pair formation. They show that the expected number of close MBH binary pairs is lower than expected, which would lead to a depletion in the LISA rate. However, as seen in Bonetti et al. (2019), including triple interactions can remedy this issue quite significantly, producing rates similar to those found when not considering a delay prescription ($\sim 23$). Additionally, including prescriptions for gas-driven migration can also increase our rates with our DA17 model (it is included in the K17 model). Without the inclusion of gas-driven migration or triple interactions, our rate is truly a lower limit. However, it can be seen that inclusion of these prescriptions would not significantly alter our results since the rate from our ND model is only slightly larger than those found with the DA17 and K17 prescriptions.

For the EAGLE simulations, Salcido et al. (2016) test seeds of $\sim 10^{5}M_{\odot}$, similar to the seed model in Illustris. Therefore, we believe two reasons other than the seed mass are the main factors in our predicted rate difference. First, Salcido et al. (2016) perform a different extraction analysis, compared to the analysis done for this paper, to deal with numerical issues similar to those suggested in section 2.1. This may inflate the true number of mergers seen in the EAGLE simulations, especially at masses near the seed mass. Second, we use more detailed binary inspiral models than the model used in Salcido et al. (2016). In their paper, they choose an inspiral time for all binaries as flat values based on if a galaxy merger was gas-rich (0.1 Gyr) or gas-poor (5 Gyr). When we compare these inspiral times with those from our more detailed models (see Figure 5), we see that our models will generally predict longer inspiral times, as well as inspiral times that are longer than the age of the Universe. Therefore, the fraction of coalescing binaries before $z=0$ is larger with their choice of inspiral model, which would lead to a higher predicted detection rate.

If we consider our rates in a relative way, even though our predicted rates are low, the comparison of our binary inspiral models can have an impact on LISA MBH science. We have shown the difference between merger models is relatively small, but still varies by a factor of two. Since most detections are from low-mass systems (see Figure 6), much of this factor of two difference will be lost in the low-mass regime. This means the K17 binary inspiral prescription will lead to more difficulty characterizing the low-mass population, as well as seeding models which will predominantly be constrained by this lower mass regime.

6 Conclusions

We have presented new MBH binary LISA rate calculations based on the Illustris cosmological simulations. Our MBH catalog is determined from the Illustris output using a new advanced extraction method allowing us to probe masses down to the simulation seed mass of $\sim 10^{5}M_{\odot}$. Previous extraction methods made a mass cut requiring $m_{1},m_{2}\geq 10^{6}M_{\odot}$. By strictly following the interaction between MBH mergers and the MBH host galaxies, we were able to refine the analysis and retain an additional 8,265 mergers in our sample ($\sim 50\%$ of mergers analyzed). This doubled the rate of predicted detections for LISA. Binaries containing an MBH of $\sim 10^{5}M_{\odot}$, especially near equal-mass binaries, are prime targets of the LISA mission due to their time spent evolving in the LISA band, their high SNR potential, and their ability to better probe MBH seeding models.

With this MBH merger catalog from Illustris, we tested four evolutionary prescriptions to understand how binary inspiral models affect our rate predictions, as well as the binary properties of the detected population. Our base model was the ND model, which represented the exact prediction from the Illustris simulation without any delays between the mergers in the simulation (at $\sim$kpc scales) and true coalescence of the MBH pairs. In order to understand the effect of our new extraction method, we also tested a model without delays while requiring $m_{1},m_{2}\geq 10^{6}M_{\odot}$. We found the detection rate was diminished by a factor of 2 when we included this mass cut.

In addition to our base models without delays, we tested two recent sub-grid models proposed in Dosopoulou & Antonini (2017) (DA17; section 3.1.1) and Kelley et al. (2017a); Kelley et al. (2017b) (K17; section 3.1.2). When comparing the evolution timescales predicted by these two models, we find that DA17 produces a smaller spread in orders of magnitude of this timescale when looking at different masses and mass ratios. K17 produces a relatively flat profile of timescales when looking at all masses and mass ratios in our catalog. When examining coalescence fractions in Figure 6, we find that DA17 favors lower total masses while K17 favors higher total masses. They both favor more equal mass ratio binaries. Additionally, it was clear from our analysis that DA17 predicts that more binaries from our catalog will merge prior to $z=0$ (84%) compared to K17 (66%).

Due to the ability to coalesce lower mass binaries, the DA17 model resembled the ND model in terms of detection rate at $\sim 0.7$ yr*-1* (integral rate calculation). The K17 model, with its lower overall coalescence fraction and inability to retain the low-mass binaries, led to a detection rate prediction similar to the ND-6 model at $\sim 0.4$ yr*-1*. The rates predicted for black hole spectroscopy were similar in magnitude to the overall detection rates. These rates represent lower limits for similar rate predictions (see Section 5).

We also examined the probability density functions of the total masses and mass ratios of the detected binaries using Monte Carlo generated catalogs. We found that all models with low-mass support (ND, DA17, and K17) produced similar detectable populations in terms of these parameters. This indicates that the surpression of sources by K17 is effectively equivalent across mass regimes, leading to a lower detection rate, while maintaining a similar total mass and mass ratio probability density function to the DA17 and ND models. The ND-6 model cannot match this because it has an entirely inaccessible region where at least one constituent MBH of less than $10^{6}M_{\odot}$ is required. Overall, we show that these two detailed models for the evolution of MBH binaries from Dosopoulou & Antonini (2017) and Kelley et al. (2017a); Kelley et al. (2017b) lead to differences in the detection rate and the observable population by LISA.

Acknowledgements

We would like to thank Marta Volonteri and the referee for helpful suggestions to improve our paper. M.L.K. acknowledges support from the National Science Foundation under grant DGE-0948017. This research was supported in part through the computational resources and staff contributions provided for the Quest/Grail high performance computing facility at Northwestern University. F.D. acknowledges support from PCTS and Lyman Spitzer Jr fellowships. L.B. acknowledges support from National Science Foundation grant AST-1715413. Astropy, a community-developed core Python package for Astronomy, was used in this research (Astropy Collaboration et al., 2013). This paper also employed use of Scipy (Jones et al., 2001), Numpy (Walt et al., 2011), and Matplotlib (Hunter, 2007).

Appendix A Effect of LISA Configuration on Rate Predictions

Here, we analyze the difference between the classic LISA sensitivity (Larson et al., 2000) and the proposed LISA sensitivity (Amaro-Seoane et al., 2017). We do this purely to draw comparisons and understand how the measurement of MBH binaries has changed with the change in mission design. For this section, we refer to the Proposal sensitivity curve (Amaro-Seoane et al., 2017) as “PL.” The classic LISA curve has an unrealistic low-frequency behavior following an $f^{-2}$ power law to infinitely low frequencies. To correct for this, we copy the low-frequency band edge behavior of PL, move it to lower strains, and spline it together with the classic LISA curve. We refer to this curve as “CLLF.” See Katz & Larson (2019) for more information about this construction. Both PL and CLLF are shown in Figure 9. The basic difference between these two curves is better low-frequency performance exhibited by CLLF due to a longer armlength (5 million km).

The CLLF LISA configuration allows for observation of sources earlier in their inspiral as well as stronger overall measurements of this signal in the low-frequency regime. Additionally, the low-frequency difference can allow for detection of larger total mass sources (Katz & Larson, 2019). The difference between these two curves can be seen visually in Figure 9 as the space between the two curves on the low-frequency end. For this paper, we study MBHs of masses greater than $\sim 10^{5}M_{\odot}$. This mass regime radiates gravitational waves observable at frequencies below $10^{-3}$ Hz. Therefore, we want to focus our sensitivity analysis on the low-frequency band edge. Katz & Larson (2019) show that the high mass range observable by LISA is strongly dependent on the low-frequency band edge behavior and can reach masses of $\sim 10^{9}M_{\odot}$.

We also tested the ability of each detector configuration to observe partial signals as well as binaries earlier in their evolution during only their inspiral phase. Figure 10 shows the detection rate of binaries versus their time-to-merger at the start of LISA observation. The binaries shown in this Figure will not reach their coalescence by the time LISA observing terminates. Therefore, the rates shown represent detections from only the inspiral portion of coalescence. Immediately to the right of the observation time, the rate drops off significantly with the loss of the merger and ringdown signals. Additionally, the inspiral signal is not observable for larger binaries, meaning LISA will not detect any of the larger mass systems if they do not merge in the observing window. At longer times before merger, the rate is further decreased as signals from lower mass systems fall below the noise. Therefore, in order for these detections to occur, they have to be lower in mass, sufficiently close in luminosity distance, and close enough in time-to-merger to be detectable. Here, we see a stark difference between the two detector configurations. The PL configuration almost drops off entirely to the right of the observation time. CLLF displays a different behavior: it shows a more gradual decrease in the detection rate with increasing start time. Detections at start times between 4 and 7 years occur at rates between $\sim 10^{-2}$ yr*-1* and 1 yr*-1*. Therefore, while our magnitude of the overall rates is low, the detection rate of these inspiral-only sources can enhance the overall detection rate by $\sim 20\%$ if the detector’s low-frequency performance is closer to the classic LISA configuration. If lower mass binaries were included, this percentage would increase because at masses lower than those tested here, we enter a regime where the inspiral stage can be observed, but the merger and ringdown are no longer detectable.

All of our delay models have the same general behavior for both LISA configurations. Since ND-6 bottoms out at $10^{6}M_{\odot}$, the difference between ND-6 and the other models can be see as the low-mass ($<10^{6}M_{\odot}$) contribution added by the advanced extraction. Consequently, ND-6 does have a much steeper drop for CLLF because it does not have the low-mass support to boost detections of inpiral-only sources. Once again, the advanced extraction helps to establish a more complete LISA analysis. Similarly, due to the inability of K17 to maintain the low-mass systems, the K17 curve is below the DA17 and ND model curves at all start times.

A.1 Monte Carlo Rate Results

Table 3 shows our Monte Carlo results for 10000 sampled catalogs. We estimate the errors in our Monte Carlo results to be approximately $\sqrt{1/N}=\sqrt{1/10000}=1\%$. The rates for the measurement of the entire signal roughly match our integral calculations (Table 4). Testing the inspiral-only results was useful in this setting because we can get a rate without assuming specific start times. We see that PL measures a negligble rate of inspiral-only sources for all of our timescale models. CLLF, on the other hand, shows that its low-frequency performance constitutes $\sim 20\%$ of all detections for models that include the lower mass binaries (9% of detections for ND-6). These systems had the lowest masses in our catalogs.