Gravitational Waves from Binary Black Hole Mergers Inside of Stars

TL;DR

This study uses numerical relativity to examine how stellar gas density affects binary black hole mergers and gravitational wave signals, ruling out certain stellar formation scenarios for GW150914.

Contribution

It provides the first detailed numerical analysis of BBH mergers inside stellar gas, demonstrating how high-density environments alter GW signals and constraining formation models.

Findings

High-density gas significantly alters GW signals during BBH mergers.

LIGO data excludes stellar gas densities above 10^7 g/cm^3 for GW150914.

Stellar core fragmentation is unlikely the formation channel for GW150914.

Abstract

We present results from a controlled numerical experiment investigating the effect of stellar density gas on the coalescence of binary black holes (BBHs) and the resulting gravitational waves (GWs). This investigation is motivated by the proposed stellar core fragmentation scenario for BBH formation and the associated possibility of an electromagnetic counterpart to a BBH GW event. We employ full numerical relativity coupled with general-relativistic hydrodynamics and set up a $30 + 30 M_\odot$ BBH (motivated by GW150914) inside gas with realistic stellar densities. Our results show that at densities $\rho \gtrsim 10^6 - 10^7 \, \mathrm{g \, cm}^{-3}$ dynamical friction between the BHs and gas changes the coalescence dynamics and the GW signal in an unmistakable way. We show that for GW150914, LIGO observations conclusively rule out BBH coalescence inside stellar gas of $\rho \gtrsim 10^7 \, \mathrm{g\,cm}^{-3}$. Typical densities in the collapsing cores of massive stars are in excess of this density. This excludes the fragmentation scenario for the formation of GW150914.

Figures

Tables

Table 1. Table 1: Model Summary. ρ 0 subscript 𝜌 0 \rho_{0} is the initial gas density, t merge subscript 𝑡 merge t_{\mathrm{merge}} the merger time, and ℳ ZDHP subscript ℳ ZDHP \mathcal{M}_{\mathrm{ZDHP}} and ℳ 150914 subscript ℳ 150914 \mathcal{M}_{150914} are the GW mismatches with the vacuum waveform for Advanced LIGO design noise and noise at the time of GW150914, respectively. For GW150914, a mismatch ℳ ≳ 0.0017 greater-than-or-equivalent-to ℳ 0.0017 \mathcal{M}\gtrsim 0.0017 becomes noticable.

Model	${\rho }_0{\left(M/60M_{\odot }\right)}^{-2}$	$t_{\mathrm{merge}}$	${ℳ}_{\mathrm{ZDHP}}$	${ℳ}_{150914}$
	$\left[{\text{g cm}}^{-3}\right]$	[ms]
G0 Vacuum	0	510	0	0
G1	$1.72\times {10}^4$	508	$8\times {10}^{-5}$	$3\times {10}^{-5}$
G2	$1.72\times {10}^5$	490	$0.0058$	$0.0016$
G3	$1.72\times {10}^6$	369	$0.1882$	$0.0665$
G4	$1.72\times {10}^7$	186	$0.3718$	$0.2386$

Equations

$$\mathcal{M}(h^1,h^2)=1-\underset{\{{\chi }_i\}}{\max }\mathcal{O}(h^1,h^2),$$

$$\mathcal{O}(h^1,h^2)=\frac{⟨h^1|h^2⟩}{\sqrt{⟨h^1|h^1⟩⟨h^2|h^2⟩}}.$$

$$⟨a|b{⟩}_f=4\mathrm{Re}{\int }_0^{\infty }\frac{\tilde{a}(f){\tilde{b}}^{\ast }(f)}{S_n(f)}df.$$

$$h(t)=F_{+}{h}_{+}+F_{\times }{h}_{\times },$$

$$⟨h^1|h^2{⟩}_{2det}=⟨h^{1,\alpha }|h^{2,\alpha }{⟩}_s+⟨h^{1,\beta }|h^{2,\beta }{⟩}_s.$$

$$⟨h^1|h^2{⟩}_{\mathrm{opt}}=⟨h_{+}^1|h_{+}^2{⟩}_s+⟨h_{\times }^1|h_{\times }^2{⟩}_s,$$

$$\mathcal{L}(d|\vec{\vartheta })\propto \exp \left[-\frac{1}{2}⟨h^{\mathrm{M}}(\vec{\vartheta })-d|h^{\mathrm{M}}(\vec{\vartheta })-d⟩\right].$$

$$\begin{array}{r}\log \mathcal{L}=\mathrm{C}-[\frac{1}{2}⟨h^{\mathrm{M}}|h^{\mathrm{M}}⟩+\frac{1}{2}⟨h^{\mathrm{GR}}|h^{\mathrm{GR}}⟩\\ +\frac{1}{2}⟨n|n⟩-⟨n|h^{\mathrm{M}}-h^{\mathrm{GR}}⟩-⟨h^{\mathrm{M}}|h^{\mathrm{GR}}⟩],\end{array}$$

$$h^{\mathrm{M}}=(1+{ϵ}_1)h^{\mathrm{GR}}+{ϵ}_2{h}^{\perp },$$

$$\begin{array}{r}\log \mathcal{L}=\log {\mathcal{L}}_0-\frac{1}{2}{ϵ}_1^2⟨h^{\mathrm{GR}}|h^{\mathrm{GR}}⟩-\frac{1}{2}{ϵ}_2^2⟨h^{\perp }|h^{\perp }⟩\\ +{ϵ}_1⟨n|h^{\mathrm{GR}}⟩+{ϵ}_2⟨n|h^{\perp }⟩,\end{array}$$

$$\mathrm{E}[\delta \log \mathcal{L}]=\frac{1}{2}{ϵ}_1^2⟨h^{\mathrm{GR}}|h^{\mathrm{GR}}⟩+\frac{1}{2}{ϵ}_2^2⟨h^{\perp }|h^{\perp }⟩.$$

$$\mathrm{E}[\delta \log \mathcal{L}]=\frac{1}{2}{ϵ}_2^2⟨h^{\perp }|h^{\perp }⟩+\mathcal{O}({ϵ}^3).$$

$$\mathcal{M}(h^{\mathrm{GR}},h^{\mathrm{M}})$$

$$=\frac{1}{2}{ϵ}_2\frac{⟨h^{\perp }|h^{\perp }⟩}{⟨h^{\mathrm{GR}}|h^{\mathrm{GR}}⟩}+\mathcal{O}({ϵ}^3),$$

$$\mathrm{E}[\delta \log \mathcal{L}]\approx {\varrho }^2\mathcal{M}.$$

$$\mathcal{M}\gtrsim \frac{1}{{\varrho }^2}.$$

$$A(t)=\sqrt{\mathrm{Re}[{\psi }_4(t){]}^2+\mathrm{Im}[{\psi }_4(t){]}^2},$$

$$\varphi (t)={\tan }^{-1}\left(\frac{\mathrm{Im}[{\psi }_4(t)]}{\mathrm{Re}[{\psi }_4(t)]}\right),$$

$$Q_s=\frac{\Delta x_3^n-\Delta x_2^n}{\Delta x_2^n-\Delta x_1^n}=0.3505,$$

Full text

Introduction. With the recent detection of the first gravitational wave (GW) events by LIGO Abbott et al. (2016a), the era of GW Astronomy has begun. An extensive multi-wavelength network of astronomical observatories is following up each candidate GW event with the hope of catching an electromagnetic (EM) counterpart. This is very well motivated for GWs from neutron star (NS) mergers (e.g., Metzger and Berger (2012)), but for observed GWs from the merger of stellar-mass binary black holes (BBHs), no EM counterpart is expected (e.g., Abbott et al. (2016b)).

However, the first observed BBH GW event, GW150914 Abbott et al. (2016a), was possibly connected with a $\gamma$-ray event seen by the Fermi satellite Connaughton et al. (2016) (though note it was not observed by other $\gamma$-ray satellites Abbott et al. (2016b); Hurley et al. (2016)). If directly related, this would be a totally unexpected EM counterpart to what was believed to be a BBH merger in pure vacuum.

To explain such an EM counterpart, Loeb Loeb (2016) proposed that the coalescing BHs formed via dynamical fragmentation in a very massive star undergoing gravitational collapse. This scenario is tentatively supported by the simulations of Reisswig et al. (2013), who found BBH formation by dynamical fragmentation in pair-unstable supermassive primordial stars. The result of Loeb’s scenario would be a BBH system embedded in high-density stellar gas whose coalescence could drive outflows giving rise to the $\gamma$-ray transient observed by Fermi.

There are arguments from stellar evolution Woosley (2016) suggesting it may be difficult to obtain collapsing stellar cores permitting dynamical fragmentation. However, this possibility is not conclusively ruled out by theory.

In this Letter, we consider the scenario in which a BBH was formed inside a collapsing massive star and conduct the first numerical relativity simulations of BBH mergers in the presence of gas with densities comparable to those in the cores of collapsing massive stars. The results of our simulations show that the GWs observed from GW150914 are inconsistent with this event having taken place inside a collapsing massive star, ruling out the dynamical fragmentation scenario.

Methods and Initial Data. We employ the open-source Einstein Toolkit and evolve Einstein’s equations in the BSSN formalism Shibata and Nakamura (1995); Baumgarte and Shapiro (1999) with fourth-order finite differences and adaptive mesh refinement (AMR). We include general-relativistic (GR) hydrodynamics in the finite-volume approach with piecewise parabolic reconstruction at cell interfaces and the Marquina flux formula for intercell fluxes Aloy et al. (1999). Inside the BH apparent horizons, we correct unphysical states using the methods detailed in Faber et al. (2007); Farris et al. (2010). Spacetime and hydrodynamics evolution are coupled in a fourth-order Runge-Kutta integrator.

For generality, we describe our setup in $G=c=1$ units and measure quantities in terms of the ADM mass $M$. We employ BBH puncture initial data and carry out a vacuum simulation (model G0) and four simulations in which we embed the BBH system in gas of constant density $\rho_{0}=\{10^{-10},10^{-9},10^{-8},10^{-7}\}M^{-2}$ initially at rest, labeled G1–G4 in the order shown. We use TwoPunctures Ansorg et al. (2004); Löffler et al. (2006) to solve for constraint satisfying quasi-circular initial data, taking into account the gas, and placing the two equal-mass, non-spinning punctures at a coordinate separation of $11.6\,M$. In the vacuum case, this corresponds to 8 orbits to merger. We employ a $\Gamma$-law equation of state $P=(\Gamma-1)\rho\epsilon$ for the gas. We set $\Gamma=4/3$ and obtain the initial $\epsilon$ by assuming a gas dominated by relativistic degenerate electrons (e.g., Shapiro and Teukolsky (1983)). We smoothly reduce $\rho$ to an atmosphere value of $10^{-16}\,M^{-2}$ outside of $80\,M$ by applying $X(R)=0.5[1+\tanh([R-80M]/15M)]$.

We employ 7 levels of AMR with the outer boundary placed at $320\,M$. The punctures are covered with a finest grid of $\Delta x=0.0195\,M$, which corresponds to approximately 45 grid points across each apparent horizon after an initial gauge adjustment. The fine grid has a linear extent of $3\,M$ to provide high resolution for the gas dynamics near the horizons. It is embedded in 5 coarser AMR levels tracking the punctures’ orbital motion. The outermost two levels are fixed. We extract GWs at $R=120\,M$ where $\Delta x=0.625\,M$ using the Newman-Penrose $\Psi_{4}$ formalism Löffler et al. (2012); Newman and Penrose (1962). We obtain the GW strain via fixed-frequency integration Reisswig and Pollney (2011).

Rescaled to a BBH mass of $M=60\,M_{\odot}$ for comparison with GW150914, each puncture has an approximate initial mass of $30\,M_{\odot}$, the initial separation is $1030\,\mathrm{km}$, with gas densities $\rho_{0}=1.72\times\{10^{4},10^{5},10^{6},10^{7}\}\,\mathrm{g\,cm}^{-3}$. The typical central density in a presupernova star is $10^{9}-10^{10}\,\mathrm{g\,cm}^{-3}$. At a radius of $1000\,\mathrm{km}$ it is in the range $10^{7}-10^{9}\,\mathrm{g\,cm}^{-3}$, depending on stellar mass (e.g., Woosley et al. (2002)). We choose $10^{7}\,\mathrm{g\,cm}^{-3}$ as the highest simulated density since it is a reasonable and arguably low value for the density of outer stellar core material left surrounding the BBH formed in dynamical fragmentation. With the above choices, the total gas mass on the computational grid is \sim$$13.8\,[M/(60\,M_{\odot})][\rho_{0}/(1.72\times 10^{7}\,\mathrm{g\,cm^{-3}})]\,M_{\odot}.

We provide a convergence study and analysis details in the Supplemental Material to this Letter.

Dynamics. In Fig. 1, we show orbital-plane snapshots of the rest-mass density at various times in model G3’s coalescence. In cgs units and for the $M=60\,M_{\odot}$ case, its merger time is $\sim 390\,\mathrm{ms}$ (we define merger time based on the peak amplitude of the $(2,2)$ GW mode). That is $\sim 142\,\mathrm{ms}$ faster than the pure-vacuum case G0.

The density colormaps in Fig. 1 reveal that soon after the start of the simulation, an ellipsoidal high-density structure surrounds the BHs. The central high-density band visually connecting the BHs is due to the gravitational focusing of gas into this region, where acceleration toward one BH is partially cancelled by the other. This feature was also observed in BBH mergers in very low-density gas (e.g., Bode et al. (2012); Farris et al. (2010) and references therein).

The ellipsoidal stucture surrounding the BBH in Fig. 1 forms because each BH accelerates the surrounding gas, dragging it along in its gravitational wake. The associated drag force, closely related to dynamical friction (e.g., Escala et al. (2004); Chandrasekhar (1943); Ostriker (1999); Dai et al. (2017); Barausse (2007); Barausse et al. (2014)), converts orbital energy into kinetic energy and internal energy of the gas (through compression and shocks). This process is what rapidly robs the BBH of its orbital energy and angular momentum. It leads to an accelerated decline of the orbital separation and an earlier merger compared to the vacuum case G0.

The BHs accrete gas during coalescence, but even in the high-density G4 case, the total mass accreted by each BH is only $\sim$4% of its initial mass. The effect of the gradually changing mass on the coalescence is much smaller than that of dynamical friction.

In Tab. 1, we summarize key properties, including the merger times, for all simulated models. The top panel of Fig. 2 shows the BBH coordinate separation $a$ as a function of simulation time. With increasing gas density, the merger is driven to earlier times. Model G1 ($\rho_{0}\sim 10^{4}\,\mathrm{g\,cm^{-3}})$ is only very mildly affected, merging some $\Delta t\sim 2$  ms earlier than the vacuum case G0. $\Delta t$ is $\sim 20\,\mathrm{ms}$, $\sim 142\,\mathrm{ms}$, and $\sim 324\,\mathrm{ms}$, for models G2, G3, and G4, respectively, which have $10$, $100$, and $1000$ times higher density than model G1. The effect of the gas on the coalescence time is thus roughly linear in density for the lower-density cases. This is qualitatively reproduced by a simple Newtonian point-particle model including GW ($\partial a/\partial t\propto a^{-3}$) and dynamical friction ($\partial a/\partial t\propto a^{5/2}$) Escala et al. (2004) terms for orbital evolution. At high density, i.e. going from G3 to G4, dynamical friction is so strong that it is no longer a linear perturbation to the GW-dominated inspiral. The point-particle model shows that in G1–G3, the dynamical friction term is always subdominant. In G4, it initially dominates over the GW term, but quickly decreases in relevance as the orbit shrinks and GW-driven evolution takes over.

Gravitational Waves. In Fig. 3, we present $h^{22}_{+}$, the real part of the $l=2,m=2$ GW mode. The low density in model G1 has a negligible effect and its waveform is essentially identical to vacuum GWs. As the density increases from G1 to G4, merger occurs progressively earlier. This leads to dramatic changes in the emitted GW train and creates an unmistakable GW signature: (1) Since all models start at the same separation, the initial GW frequency is f_{0}\sim$$26\,\mathrm{Hz} for all models. (2) Dynamical friction shortens the inspiral, leading to a faster sweep (“chirp”) of the GWs through frequency space. (3) The faster chirp is not due to a substantial increase of the BBH mass. Hence, merger and ringdown GW emission is at frequencies that change only mildly with $\rho_{0}$. Model G0 has a final BH mass of \sim$$57.6\,M_{\odot}, a dimensionless spin $a^{\star}=0.69$, and we find a ringdown GW frequency of \sim$$296\,\mathrm{Hz}. The highest-density model G4 produces a merged BH of \sim$$64.8\,M_{\odot} and $a^{\star}=0.65$, consistent (see, e.g., Berti et al. (2009)) with its ringdown GW frequency of \sim$$265\,\mathrm{Hz}.

Analysis and Observational Consequences. We compute the GW mismatch $\mathcal{M}(h_{i},h_{j})$ (see Damour et al. (1998); Reisswig and Pollney (2011) and Supplemental Material) for each waveform G1–G4 with the G0 vacuum case. GW mismatch takes into account the detector noise spectrum and we consider frequencies in the interval $[26,3000]\,\mathrm{Hz}$. We employ Advanced LIGO design noise Shoemaker (2010) ($\mathcal{M}_{\mathrm{ZDHP}}$) and the noise around GW150914 Abbott et al. (2016a) ($\mathcal{M}_{150914}$). $\mathcal{M}$ is in $[0,1]$ and $\mathcal{M}=0$ means $h_{i}$ and $h_{j}$ are identical. For an observation with signal-to-noise ratio $\varrho$, an $\mathcal{M}\gtrsim 1/\varrho^{2}$ leads to observational inconsistencies (see Lindblom et al. (2008) and Supplemental Material). GW150914 was observed with $\varrho\sim 24$, so $\mathcal{M}\gtrsim 0.0017$ will become noticable.

We summarize $\mathcal{M}_{\mathrm{ZDHP}}$ and $\mathcal{M}_{\mathrm{150914}}$ for all models in Tab. 1. The results for $\mathcal{M}_{\mathrm{150914}}$ show that for GW150914, densities $\rho_{0}\lesssim 10^{4}\,\mathrm{g\,cm}^{-3}$ (G1) are indistiguishable from vacuum. Model G2 ($\rho_{0}\sim 10^{5}\,\mathrm{g\,cm}^{-3}$) is marginally distinguishable. The situation is very different for models G3 and G4 with $\mathcal{M}_{\mathrm{150914}}$ \sim$$0.07 and \sim$$0.24, respectively. These results show that stellar densities $\rho_{0}\gtrsim 10^{6}\,\mathrm{g\,cm}^{-3}$ lead to highly significant inconsistencies with vacuum.

An additional possibility is that the G1–G4 waveforms could have lower mismatch with vacuum waveforms of BBHs with different parameters. We explore this with a 7-dimensional numerical relativity surrogate GW model Blackman et al. (2017a, b), covering BBH mass ratio $q$ (up to $q=2$) and six spin components (up to effective spin $\chi_{\mathrm{eff}}=(M_{1}a^{*}_{1}+M_{2}a^{*}_{2})/M=0.8$), assuming zero eccentricity. We vary parameters to minimize $\mathcal{M}_{\mathrm{150914}}$ and find $0.011$ and $0.061$, for model G3 and G4, respectively. For G3, the minimum $\mathcal{M}_{\mathrm{150914}}$ is at $M=70.6\,M_{\odot}$, $q\simeq 1.0$, and effective spin $\chi_{\mathrm{eff}}\simeq 0.17$. For G4, we find $M=75.3\,M_{\odot}$, $q\simeq 1.6$, and $\chi_{\mathrm{eff}}\simeq-0.47$. Even with the reduction in $\mathcal{M}_{\mathrm{150914}}$, $\rho_{0}\gtrsim 10^{7}\,\mathrm{g\,cm}^{-3}$ leads to observable differences with any waveform covered by the surrogate model.

Having established that an equal mass, non-spinning BBH merger in stellar-density gas with $\rho_{0}\gtrsim 10^{7}\,\mathrm{g\,cm}^{-3}$ is inconsistent with GW150914 and all BBH waveforms from our surrogate model, there remains the following crucial question: Are there BBH parameter choices that could make a merger in gas appear just like GW150914?

We argue that the answer is ’No’: The observational BBH parameter space encompasses total mass, mass ratio, eccentricity, and six spin components. (i) BBHs of lower $M$ inspiral to higher frequencies and have more cycles from $\sim$26 Hz to merger. Using the surrogate, we find that $M=43.7\,M_{\odot}$ extends the G0 case by $\sim$324 ms, the difference in merger times between G0 and G4. However, its time-frequency evolution and ring-down frequency (\sim$$400\,\mathrm{Hz}) are substantially different from GW150914, leading to large mismatch. (ii) Dynamical fragmentation in our scenario leads to near-equal-mass fragments (e.g., Reisswig et al. (2013)). We consider $q=2$ as an extreme limit. In the vacuum case, it extends the inspiral by \sim$$38\,\mathrm{ms} Blackman et al. (2017b), insufficient to compensate for the gas effect. (iii) High BH spin causes “orbital hangup.” The effect is largest for equal spins aligned with the orbital angular momentum. Using our surrogate and the SpEC waveforms SXS , we find that for $a^{*}=0.99$ ($a^{*}=0.4$), inspiral is prolonged by $177\,\mathrm{ms}$ ($71\,\mathrm{ms}$). The effect is linear in $a^{*}$. To explore the effect of spin in the stellar-density G4 case, we carry out a simulation with $a^{*}=0.4$ for both BHs. We find that merger is delayed by \sim$$17.2\,\mathrm{ms}. Extrapolating to $a^{*}=0.99$ from the vacuum case, spin could extend the G4 inspiral by at most \sim$$39\,\mathrm{ms}. This is insufficient to mimic GW150914.

Discussion and Conclusions. Fragmentation of a massive star’s core into clumps that collapse further to NSs or BHs is an interesting scenario for the formation of NS binaries and BBHs (e.g., Bonnell and Pringle (1995); Davies et al. (2002); Postnov et al. (2016); Reisswig et al. (2013)). While perhaps unlikely (e.g., Fryer and Warren (2004); Woosley (2016)), this scenario has not previously been ruled out observationally. As proposed by Loeb Loeb (2016), it would endow a BBH merger with the gas necessary to produce an EM counterpart. Dai et al. Dai et al. (2017) suggested, but did not show, that the gas surrounding the BBH could have observable consequences in the emitted GWs.

We employed numerical relativity coupled with GR hydrodynamics for a controlled experiment into the effects of stellar-density gas on BBH mergers. Scaled to a total system mass of $60\,M_{\odot}$ (consistent with GW150914), our results show that dynamical friction between the BHs and gas at stellar densities \rho_{0}\gtrsim$$10^{6}-10^{7}\,\mathrm{g\,cm}^{-3} profoundly affects the coalescence dynamics, drastically shortening the time to merger. This modifies the resulting GW signal in an unmistakable way, leading to differences with vacuum waveforms that can be observed by LIGO.

Our analysis furthermore suggests that it is not possible to choose BBH parameters that would yield a waveform in stellar-density gas resembling GW150914Abbott et al. (2017). Thus we conclude that it is highly unlikely that GW150914 was formed through dynamical fragmentation in a massive star and Loebs scenario Loeb (2016) is ruled out by the GW observation alone.

Future work should address the limitations of our work: We assumed the gas to be non-magnetized and initially at rest, but angular momentum and magnetic fields can have dynamical impact. We employed a constant density, but real stars have radially varying density. Finally, we used a $\Gamma$-law equation of state, ignoring microphysics such as electron capture, neutrinos, and nuclear dissociation, which all may have effects on the gas dynamics. In our analysis, we did not consider GW detector calibration uncertainties of \sim$$10\% Abbott et al. (2016a). This should affect all waveforms equally and is unlikely to alter our conclusions.

We provide waveforms and additional visualizations of our simulations at https://stellarcollapse.org/bbhgas. We thank M. Sasaki, G. Domènech, K. Kiuchi, M. Shibata, K. Ioka, T. Tanaka, E. Schnetter, E. Firing, T. Bogdanovic, and N. Deruelle for discussions. This research is partially supported by MEXT, IRU-AFS, NSF grants ACI-1550514, CAREER PHY-1151197, and PHY-1404569, and ERC-2014-CoG 646597, MSCA-RISE-2015 690904, and STFC ST/L000636/1. We used the matplotlib Python package Hunter (2007) for the figures. The simulations were performed on the cluster Wheeler, supported by the Sherman Fairchild Foundation and Caltech, and on supercomputers of the NSF XSEDE network under allocation TG-PHY100033 and TG-PHY090003. This paper has Yukawa Institute report number YITP-17-40.

I Supplemental Material

II Analysis Details

Gravitational Wave Mismatch Calculation. The mismatch between two observed waveforms $h^{1}(t)$ and $h^{2}(t)$ is defined as one minus the maximum overlap $\mathcal{O}(h^{1},h^{2})$,

[TABLE]

where the overlap is given by

[TABLE]

Here, $\langle\cdot|\cdot\rangle$ is a detector-noise weighted inner product and optimization is carried out over a set $\{\chi_{i}\}$ of parameters impacting the overlap (e.g., shifts in waveform phases, polarization angles, arrival times) Damour et al. (1998).

In the simplest case, we can choose $\langle\cdot|\cdot\rangle$ as the frequency-domain noise weighted inner product Cutler and Flanagan (1994),

[TABLE]

Here, $S_{n}(f)$ is the detector noise power spectral density and $\tilde{a}(f)$ is the Fourier transform of $a(t)$.

The real gravitational wave signal $h(t)$ observed by a single detector is given by

[TABLE]

where $F_{+}$ and $F_{\times}$ are the detector antenna pattern functions that depend on the sky location of the source and polarization basis (see, e.g., Harry and Fairhurst (2011)).

We now consider two scenarios: (1) A best case in which both $h_{+}$ and $h_{\times}$ are measured by two optimally oriented GW detectors at Advanced LIGO design sensitivity (“ZDHP” for zero-detuning, high-power Shoemaker (2010)). (2) The realistic scenario of the two Advanced LIGO interferometers with the sensitivity at the time of GW150914.

For both cases, we need the two-detector inner product for two detectors $\alpha$ and $\beta$, which is defined Harry and Fairhurst (2011) as the sum of the single-detector contributions,

[TABLE]

Here, $h^{1,\alpha}$ is waveform 1 as seen by detector $\alpha$ through Eq. 4 and so forth. The single-detector inner product $\langle\cdot,\cdot\rangle_{s}$ used is that given by Eq. 3 with the exception that we integrate over some frequency interval defined by $[f_{\mathrm{min}},f_{\mathrm{max}}]$. In practice, we obtain the necessary Fourier transforms by using the Fast Fourier Transform algorithm after tapering the ends of the time domain signal and padding with zeros for all waveforms to have the same length in the time domain.

For scenario (1), we follow Blackman et al. (2017a) and define an optimal two-detector $\mathcal{O}_{\mathrm{opt}}$ overlap by choosing detectors oriented so that one detector is maximally sensitive to $h_{+}$ (and insensitive to $h_{\times}$) while the opposite is true for the other detector. We then have

[TABLE]

with $S_{n}(|f|)$ in Eq. 3 chosen as the Advanced LIGO ZDHP noise power spectral density. $\mathcal{O}_{\mathrm{opt}}$ is then given by Eq. 2 with $\langle\cdot|\cdot\rangle_{\mathrm{opt}}$ and the mismatch is obtained as $\mathcal{M}_{\mathrm{ZDHP}}=1-\max\mathcal{O}_{\mathrm{opt}}$. We optimize over time shifts and polarization angle shifts of the waveforms. Since we consider only the $(2,2)$ GW mode, we simply assume a face-on direction of GW propagation, and orbital phase shifts are identical to polarization phase shifts. See Blackman et al. (2017a) for further details.

For scenario (2), we use the inner product of Eq. 5 with the Advanced LIGO Hanford and Livingston antenna patterns Althous et al. (1997) for GW150914 and the parameters given in Abbott et al. (2017). We employ the actual Hanford and Livingston noise power spectral densities at the time of GW150914 provided at https://losc.ligo.org/events/GW150914/. We obtain $\mathcal{M}_{\mathrm{GW150914}}=1-\max\mathcal{O}_{\mathrm{GW150914}}$ for the $(2,2)$ GW mode by optimizing over time shifts, polarization angle shifts, and orbital phase shifts. We neglect contributions from other GW modes.

Reduction in Log-Likelihood due to Mismatch. In GW parameter estimation, the posterior probability of a BBH parameter vector $\vec{\vartheta}$ is determined from the prior and likelihood. The GW likelihood function (e.g., Abbott et al. (2016)) is given by

[TABLE]

Here, $d=h^{\mathrm{GR}}+n$ is the data observed in the detectors consisting of the GR signal (we use “GR” as a synonym for “true”) and detector noise $n$. $h^{\mathrm{M}}$ is the template waveform generated by some waveform model.

The log-likelihood is then

[TABLE]

where $C$ is a constant of proportionality.

Suppose that $h^{\mathrm{M}}$ is different from the true signal,

[TABLE]

where $\langle h^{\bot}|h^{\mathrm{GR}}\rangle=0$. Here $\epsilon_{1}$ and $\epsilon_{2}$ are numbers and we consider the limit $\epsilon_{1,2}\ll 1$. Any $h^{\mathrm{M}}$ can be decomposed in this way. The log-likelihood becomes

[TABLE]

where $\log\mathcal{L}_{0}$ is the log-likelihood when $h^{M}=h^{\mathrm{GR}}$. The expected reduction in the log-likelihood is then

[TABLE]

We now allow a small bias in the distance to the source by rescaling $h^{\mathrm{M}}$ by $(1+\epsilon_{1})^{-1}$ with which we obtain the convenient expression

[TABLE]

The mismatch between $h^{\mathrm{GR}}$ and $h^{\mathrm{M}}$ is

[TABLE]

where optimization over phase shifts, time shifts, etc. is implicit.

The signal-to-noise ratio $\varrho$ is given by $\varrho^{2}=\langle h^{\mathrm{GR}}|h^{\mathrm{GR}}\rangle$. With this, we find

[TABLE]

The posterior probability will be affected by a factor of Euler’s number $e$ when $\delta\log\mathcal{L}=1$, which can be considered a mild observational inconsistency. Hence, the mismatch $\mathcal{M}$ will begin to have an effect on GW data analysis when

[TABLE]

III Numerical Convergence

We carry out additional simulations at coarse-grid resolutions $\Delta x_{1}=1.00\,M$ and $\Delta x_{3}=1.60\,M$, in addition to our standard-resolution simulations of $\Delta x_{2}=1.25\,M$. For our convergence analysis, we choose the vacuum (G0) and the highest density (G4) as two extremes of the simulations we carry out. We focus our analysis on the gravitational waveforms since these are the most important output of our simulations.

In Fig. 4, we show numerical convergence in the Newman-Penrose scalar $\psi_{4}$ between the different resolutions for the G0 vacuum simulation. We consider phase and amplitude differences separately. The amplitude is defined as

[TABLE]

while the phase is defined as

[TABLE]

where $\mathrm{Re}[\psi_{4}]$ and $\mathrm{Im}[\psi_{4}]$ are the real and imaginary parts of $\psi_{4}$, respectively. Our numerical scheme is fourth-order, hence, we expect fourth-order convergence and a self-convergence factor of

[TABLE]

where $n$ is the order of convergence. In Fig. 4, we rescale the differences between highest resolution and second-highest (i.e. standard) resolution by $1/Q_{s}$. These rescaled curves lie essentially on top of the curves for the differences between second-highest and lowest resolution, demonstrating approximate fourth-order convergence.

In Fig. 5 we perform the same analysis for the highest-density simulation G4. In this case, the hydrodynamics plays an important role in driving the coalescence. If our finite-volume implementation dominates the numerical error, we expect second-order convergence when the flow is smooth. However, soon after the start of the simulation, steep density gradients and shocks develop for which our numerical scheme (as any high-resolution shock capturing scheme) is only first-order convergent. Hence, we can only expect first-order convergence. We compute a first-order self-convergence factor $Q_{s}=0.7143$, with $1/Q_{s}=2.85$. Figure 5 shows that we obtain roughly first-order convergence in GW amplitude and phase.

In order to clarify how numerical resolution effects the main results of our paper, we have calculate mismatches between various resolutions for the G0 and G4 cases. For the G0 case we find the mismatches to be $1.6\times 10^{-3}$ between high and medium resolution, $2.9\times 10^{-3}$ between medium and low resolution, and $3.5\times 10^{-3}$ between high and low resolution. For the G4 case the mismatches are $3.5\times 10^{-5}$ between high and medium resolution, $1.0\times 10^{-4}$ between medium and low resolution, and $2.3\times 10^{-4}$ between high and low resolution. Comparing these results with the mismatches listed in Table 1 (in the main paper), we conclude that our main conclusions are independent of numerical resolution.