Extracting equation of state parameters from black hole-neutron star mergers. I. Nonspinning black holes
Benjamin D. Lackey, Koutarou Kyutoku, Masaru Shibata, Patrick R. Brady, John L. Friedman
EExtracting equation of state parameters from black hole-neutron star mergers. I.Nonspinning black holes
Benjamin D. Lackey , Koutarou Kyutoku , Masaru Shibata , Patrick R. Brady , John L. Friedman Department of Physics, University of Wisconsin–Milwaukee, P.O. Box 413, Milwaukee, WI 53201, USA Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan
The late inspiral, merger, and ringdown of a black hole-neutron star (BHNS) system can provideinformation about the neutron-star equation of state (EOS). Candidate EOSs can be approximatedby a parametrized piecewise-polytropic EOS above nuclear density, matched to a fixed low-densityEOS; and we report results from a large set of BHNS inspiral simulations that systematically varytwo parameters. To within the accuracy of the simulations, we find that, apart from the neutron-star mass, a single physical parameter Λ, describing its deformability, can be extracted from thelate inspiral, merger, and ringdown waveform. This parameter is related to the radius, mass, and (cid:96) = 2 Love number, k , of the neutron star by Λ = 2 k R / M , and it is the same parameterthat determines the departure from point-particle dynamics during the early inspiral. Observationsof gravitational waves from BHNS inspiral thus restrict the EOS to a surface of constant Λ in theparameter space, thickened by the measurement error. Using various configurations of a singleAdvanced LIGO detector, we find that neutron stars are distinguishable from black holes of thesame mass and that Λ / or equivalently R can be extracted to 10–40% accuracy from single eventsfor mass ratios of Q = 2 and 3 at a distance of 100 Mpc, while with the proposed Einstein Telescope,EOS parameters can be extracted to accuracy an order of magnitude better. PACS numbers: 97.60.Jd, 26.60.Kp, 95.85.Sz
I. INTRODUCTION
Construction of the second-generation Advanced LIGO(aLIGO) detectors is underway, and will soon begin forAdvanced VIRGO and LCGT, making it likely that grav-itational waveforms from compact binaries will be ob-served in this decade. Plans are also in development forthe third generation Einstein Telescope (ET) detectorwith an order-of-magnitude increase in sensitivity overaLIGO. Population synthesis models predict that witha single aLIGO detector binary neutron star (BNS) sys-tems will be observed with a signal-to-noise ratio (SNR)of 8, at an event rate between 0.4 and 400 times per yearand with a most likely value of 40 per year [1]. Blackhole–neutron star (BHNS) systems are also expected, butwith a more uncertain rate of between 0.2 and 300 eventsper year at the same SNR and with a most likely valueof 10 events per year for a canonical 1.4 M (cid:12) –10 M (cid:12) sys-tem [1]. The expected mass ratios Q = M BH /M NS ofBHNS systems are also highly uncertain and may rangefrom just under 3 to more than 20 [2, 3].A major goal of the gravitational-wave (GW) programis to extract from observed waveforms the physical char-acteristics of their sources and, in particular, to use thewaveforms of inspraling and merging BNS and BHNSsystems to constrain the uncertain EOS of neutron-starmatter. During inspiral the tidal interaction betweenthe two stars leads to a small drift in the phase of thegravitational waveform relative to a point particle sys-tem. Specifically the tidal field E ij of one star will in-duce a quadrupole moment Q ij in the other star given by Q ij = − λ E ij where λ is an EOS dependent quan-tity that describes how easily the star is distorted. Amethod for determining λ for relativistic stars was foundby Hinderer [4]; its effect on the waveform was calculatedto Newtonian order (with the relativistic value of λ ) byFlanagan and Hinderer [5] and to first post-Newtonian(PN) order by Vines, Flanagan, and Hinderer [6, 7]. Thistidal description has also been extended to higher ordermultipoles [8, 9].The detectability of EOS effects have been examinedfor both BNS and BHNS systems using this analyticaldescription of the inspiral. For BNS systems, the de-tectability of λ with aLIGO was examined for polytropicEOS [5] as well as a range of theoretical EOS commonlyfound in the NS literature for aLIGO and ET [10]. Thesestudies considered only the waveform up to frequencies of400–500 Hz ( ∼ M (cid:12) equal-mass NSs). For this early part of inspiral, theyfind that the tidal deformability is detectable by aLIGOonly for an unusually stiff EOS and for low neutron-starmasses ( < . M (cid:12) ). ET on the other hand would have anorder of magnitude improvement in estimating λ , allow-ing ET to distinguish between different classes of EOS.For BHNS systems, using the recently calculated 1PNcorrections, Pannarale et al. [11] examined detectabilityfor a range of mass ratios, finding that aLIGO will beable to distinguish between BHNS and binary black hole The tidal deformability for the (cid:96) th multipole is often defined interms of the NS radius R and its dimensionless (cid:96) th Love number k (cid:96) by λ (cid:96) = (cid:96) − G k (cid:96) R (cid:96) +1 . Here we will discuss only the (cid:96) = 2term so we write λ := λ . a r X i v : . [ a s t r o - ph . H E ] S e p (BBH) systems only for low mass ratios and stiff EOSwhen considering the full inspiral waveform up to thepoint of tidal disruption.In sharp contrast to these analytic post-Newtonianresults, analysis of just the last few orbits of BNS in-spiral from numerical simulations has shown that theNS radius may be extracted to a higher accuracy, of O (10%) [12], and this is confirmed by a study based ona set of longer and more accurate waveforms from twodifferent codes [13]. In addition, comparisons betweenthe analytical tidal description and BNS quasiequilib-rium sequences [14] as well as long BNS numerical wave-forms [15, 16] suggest that corrections beyond the 1PNquadrupole description are significant and substantiallyincrease the tidal effect during the late inspiral.Numerical BHNS simulations have also been done toexamine the dependence of the waveform on mass ratio,BH spin, NS mass, and the neutron-star EOS [17–28].However, an analysis of the detectability of EOS infor-mation with GW detectors using these simulations hasnot yet been done, and the present paper presents thefirst results of this kind. EOS information from tidal in-teractions is present in the inspiral waveform. For BHNSsystems, however, the stronger signal is likely to arisefrom a sharp drop in the GW amplitude arising fromtidal disruption prior to merger or, when there is negligi-ble disruption, from the cutoff frequency at merger [29].We find from simulations of the last few orbits, merger,and ringdown of BHNS systems with varying EOS that,to within numerical accuracy, the EOS parameter ex-tracted from the waveform is the same tidal parameterΛ that determines the departure from point particle be-havior during inspiral; here Λ is a dimensionless versionof the tidal parameter:Λ := Gλ (cid:18) c GM NS (cid:19) = 23 k (cid:18) c RGM NS (cid:19) , (1)where k is the quadrupole Love number.The constraint on the EOS imposed by gravitational-wave observations of BHNS inspiral and merger is essen-tially a restriction of the space of EOS p = p ( ρ ) to ahypersurface of constant Λ, thickened by the uncertaintyin the measurement (that is, a restriction to the set ofEOS for which a spherical neutron star of the mass ob-served in the inspiral has tidal parameter Λ). We usea parametrized EOS based on piecewise polytropes [30],to delineate this region in the EOS space, but the resultcan be used to constrain any choice of parameters for theEOS space.In Sec. II we discuss the parametrized EOS used inthe simulations. We give in Sec. III an overview of thenumerical methods used and, in Sec. IV, a description ofthe waveforms from the simulations. We then discuss theanalytical waveforms used for the early inspiral and is-sues related to joining the analytical and numerical wave-forms to create hybrids in Sec. V, and we then estimatethe uncertainty in extracting EOS parameters in Sec. VI. Finally, we discuss future work in Sec. VII. In the appen-dices we discuss methods for numerically evaluating theFisher matrix, and we provide instructions for generatingeffective one body (EOB) waveforms. In a second paperwe will consider the detectability of EOS parameters forBHNS systems with spinning BHs. Conventions : Unless otherwise stated we set G = c =1. Base 10 and base e logarithms are denoted log and lnrespectively. We define the Fourier transform ˜ x ( f ) of afunction x ( t ) by˜ x ( f ) = (cid:90) ∞−∞ x ( t ) e − πift dt, (2)and the inverse Fourier transform by x ( t ) = (cid:90) ∞−∞ ˜ x ( f ) e πift df. (3) II. PARAMETRIZED EOS
To understand the dependence of the BHNS waveformon the EOS we systematically vary the free parameters ofa parametrized EOS and then simulate a BHNS inspiralfor each set of parameters. We choose the piecewise poly-tropic EOS introduced in Ref [30]. Within each densityinterval ρ i − < ρ < ρ i , the pressure p is given in termsof the rest mass density ρ by p ( ρ ) = K i ρ Γ i , (4)where the adiabatic index Γ i is constant in each interval,and the pressure constant K i is chosen so that the EOSis continuous at the boundaries ρ i between adjacent seg-ments of the EOS. The energy density (cid:15) is found usingthe first law of thermodynamics, d (cid:15)ρ = − pd ρ . (5)Ref. [30] uses a fixed low density EOS for the NS crust.The parametrized high density EOS is then joined ontothe low density EOS at a density ρ that depends on thevalues of the high-density EOS parameters. The high-density EOS consists of a three-piece polytrope with fixeddividing densities ρ = 10 . g/cm and ρ = 10 g/cm between the three polytropes. The resulting EOS hasfour free parameters. The first parameter, the pressure p at the first dividing density ρ , is closely related to theradius of a 1.4 M (cid:12) NS [31]. The other three parametersare the adiabatic indices { Γ , Γ , Γ } for the three den-sity intervals. This parametrization accurately fits a widerange of theoretical EOS and reproduces the correspond-ing NS properties such as radius, moment of inertia, andmaximum mass to a few percent [30].Following previous work on BNS [12] and BHNS sim-ulations [17, 28] we use a simplified two-parameter ver-sion of the piecewise-polytrope parametrization and uni-formly vary each of these parameters. For our two pa-rameters we use the pressure p as well as a single fixedadiabatic index Γ = Γ = Γ = Γ for the core. The crustEOS is given by a single polytrope with the constants K = 3 . × in cgs units and Γ = 1 . g/cm is 1 . × dyne/cm .(For most values of p , Γ , and Γ , the central density ofa 1.4 M (cid:12) star is below or just above ρ , so the param-eter Γ is irrelevant anyway for BNS before merger andBHNS for all times.)We list in Table I the 21 EOS used in the simulationsalong with some of the NS properties. In addition, weplot the EOS as points in parameter space in Fig. 1 alongwith contours of constant radius, tidal deformability Λ,and maximum NS mass. The 1.93 M (cid:12) maximum masscontour corresponds to the recently observed pulsar witha mass of 1 . ± . M (cid:12) measured using the Shapirodelay [32]. In this two-parameter cross section of the fullfour-parameter EOS space, parameters below this curveare ruled out. (cid:71) l og (cid:72) p (cid:76) (cid:64) d y ne (cid:144) c m (cid:68) M m a x (cid:60) . (cid:32) M (cid:159) . (cid:32) M (cid:159) R (cid:61) (cid:32) k m R (cid:61) (cid:32) k m R (cid:61) (cid:32) k m R (cid:61) (cid:32) k m (cid:76) (cid:61) (cid:76) (cid:61) (cid:76) (cid:61) (cid:76) (cid:61) FIG. 1: The 21 EOS used in the simulations are representedby blue points in the parameter space. For a NS of mass1.35 M (cid:12) , contours of constant radius are solid blue and con-tours of constant tidal deformability Λ are dashed red. Alsoshown are dotted contours of maximum NS mass. The shadedregion does not allow a 1.35 M (cid:12) NS.
III. NUMERICAL METHODS
We employ BHNS binaries in quasiequilibruim statesfor initial conditions of our numerical simulations. Wecompute a quasiequilibrium state of the BHNS binary asa solution of the initial value problem of general relativ-ity, employing the piecewise polytopic EOS described inthe previous section. The details of the formulation andnumerical methods are described in Refs. [17, 33]. Com-putations of the quasiequilibrium states are performedusing the spectral-method library
LORENE [34].Numerical simulations are performed using anadaptive-mesh refinement code
SACRA [35].
SACRA solves the Einstein evolution equations in the BSSN formalismwith the moving puncture gauge, and solves the hydro-dynamic equations with a high-resolution central scheme.The formulation, the gauge conditions, and the numericalscheme are the same as those described in Ref. [17]. Forthe EOS, we decompose the pressure and energy densityinto cold and thermal parts as p = p cold + p th , (cid:15) = (cid:15) cold + (cid:15) th . (6)We calculate the cold parts of both variables using thepiecewise polytropic EOS from ρ , and then the thermalpart of the energy density is defined from (cid:15) as (cid:15) th = (cid:15) − (cid:15) cold . Because (cid:15) th vanishes in the absence of shockheating, (cid:15) th is regarded as the finite temperature part.In our simulations, we adopt a Γ-law ideal gas EOS p th = (Γ th − (cid:15) th , (7)to determine the thermal part of the pressure, and chooseΓ th equal to the adiabatic index in the crust region, Γ ,for simplicity.In our numerical simulations, gravitational waves areextracted by calculating the outgoing part of the Weylscalar Ψ at finite coordinate radii ∼ M (cid:12) , and byintegrating Ψ twice in time as h + ( t ) − ih × ( t ) = (cid:90) t −∞ dt (cid:48) (cid:90) t (cid:48) −∞ dt (cid:48)(cid:48) Ψ ( t (cid:48)(cid:48) ) . (8)In this work, we perform this time integration with a“fixed frequency integration” method to eliminate un-physical drift components in the waveform [36]. In thismethod, we first perform a Fourier transformation of Ψ as ˜Ψ ( f ) = (cid:90) ∞−∞ dt Ψ ( t ) e − πift . (9)Using this, Eq. (8) is rewritten as h + ( t ) − ih × ( t ) = − π ) (cid:90) ∞−∞ df ˜Ψ ( f ) f e πift . (10)We then replace 1 /f of the integrand with 1 /f for | f | < f , where f is a free parameter in this method.By appropriately choosing f , this procedure suppressesunphysical, low-frequency components of gravitationalwaves. As proposed in Ref. [36], we choose f to be ∼ . m Ω / π , where Ω is the initial orbital angular ve-locity and m (= 2) is the azimuthal quantum number. IV. DESCRIPTION OF WAVEFORMS
Using the 21 EOS described in Table I, we have per-formed 30 BHNS inspiral and merger simulations withdifferent mass ratios Q = M BH /M NS and neutron starmasses M NS . A complete list of these simulations is TABLE I: Neutron star properties for the 21 EOS used in the simulations. The original EOS names [12, 17, 28] are also listed. p is given in units of dyne/cm , maximum mass is in M (cid:12) , and neutron star radius R is in km. R , k , and Λ are given forthe two masses used: { . , . } M (cid:12) . The values listed for log p are rounded to three digits. The exact values used in thesimulations can be found by adding log( c/ cm s − ) − . ≈ . p Γ M max R . k , . Λ . R . k , . Λ . p.3Γ2.4 Bss 34.3 2.4 1.566 10.66 0.0765 401 10.27 0.0585 142p.3Γ2.7 Bs 34.3 2.7 1.799 10.88 0.0910 528 10.74 0.0751 228p.3Γ3.0 B 34.3 3.0 2.002 10.98 0.1010 614 10.96 0.0861 288p.3Γ3.3 34.3 3.3 2.181 11.04 0.1083 677 11.09 0.0941 334p.4Γ2.4 HBss 34.4 2.4 1.701 11.74 0.0886 755 11.45 0.0723 301p.4Γ2.7 HBs 34.4 2.7 1.925 11.67 0.1004 828 11.57 0.0855 375p.4Γ3.0 HB 34.4 3.0 2.122 11.60 0.1088 872 11.61 0.0946 422p.4Γ3.3 34.4 3.3 2.294 11.55 0.1151 903 11.62 0.1013 454p.5Γ2.4 34.5 2.4 1.848 12.88 0.1000 1353 12.64 0.0850 582p.5Γ2.7 34.5 2.7 2.061 12.49 0.1096 1271 12.42 0.0954 598p.5Γ3.0 H 34.5 3.0 2.249 12.25 0.1165 1225 12.27 0.1029 607p.5Γ3.3 34.5 3.3 2.413 12.08 0.1217 1196 12.17 0.1085 613p.6Γ2.4 34.6 2.4 2.007 14.08 0.1108 2340 13.89 0.0970 1061p.6Γ2.7 34.6 2.7 2.207 13.35 0.1184 1920 13.32 0.1051 932p.6Γ3.0 34.6 3.0 2.383 12.92 0.1240 1704 12.97 0.1110 862p.6Γ3.3 34.6 3.3 2.537 12.63 0.1282 1575 12.74 0.1155 819p.7Γ2.4 34.7 2.4 2.180 15.35 0.1210 3941 15.20 0.1083 1860p.7Γ2.7 34.7 2.7 2.362 14.26 0.1269 2859 14.25 0.1144 1423p.7Γ3.0 1.5H 34.7 3.0 2.525 13.62 0.1313 2351 13.69 0.1189 1211p.7Γ3.3 34.7 3.3 2.669 13.20 0.1346 2062 13.32 0.1223 1087p.9Γ3.0 2H 34.9 3.0 2.834 15.12 0.1453 4382 15.22 0.1342 2324 given in Table II. For the mass ratio Q = 2 and NSmass M NS = 1 . M (cid:12) , we performed a simulation foreach of the 21 EOS. In addition, we performed simula-tions of a smaller NS mass ( Q = 2, M NS = 1 . M (cid:12) )and a larger mass ratio ( Q = 3, M NS = 1 . M (cid:12) ),in which we only varied the pressure p over the range34 . ≤ log( p / (dyne cm − )) ≤ . . l = m = 2 mode, while the disrupted TABLE II: Data for the 30 BHNS simulations. NS mass isin units of M (cid:12) , and Ω M is the angular velocity used in theinitial data where M = M BH + M NS . Q M NS EOS Ω M Q M NS EOS Ω M matter is roughly axisymmetric as it accretes onto theblack hole; and the accretion timescale of the spread-outmatter is long compared to the periods of the dominantmodes.The dependence of the waveform on the EOS canbe seen more clearly by decomposing each waveforminto amplitude A ( t ) and phase Φ( t ) with the relation h + ( t ) − ih × ( t ) = A ( t ) e − i Φ( t ) . In Fig. 3, the amplitudeas a function of time for each BHNS waveform is com-pared to a BBH waveform of the same value of Q and M NS . At early times, the waveform is almost identicalto the BBH waveform. However, a few ms before themaximum amplitude is reached, the amplitude begins todepart from the BBH case. For each Q and M NS , thisdeparture from the BBH waveform is monotonic in Λ.Neutron stars with large values of Λ merge earlier, andas a result the waveforms reach a smaller maximum am-plitude. The phase of each waveform is compared to thatof the EOB BBH waveform Φ EOB in Fig. 4. At earlytimes the phase oscillates about the EOB phase due toinitial eccentricity in the numerical waveform discussedin Sec. V B. At later times, closer to the merger, tidal in-teractions lead to a higher frequency orbit; this, togetherwith correspondingly stronger gravitational wave emis-sion, means the BHNS phase accumulates faster than theEOB phase. This continues for 1–2 ms after the wave-form reaches its maximum amplitude (indicated by thedot on each curve). Eventually the amplitude drops sig-nificantly, and numerical errors dominate the phase. Wetruncate the curves when the amplitude drops below 0.01.The monotonic dependence of the waveform on Λ canagain be seen in its Fourier transform ˜ h , shown in Figs. 5and 6, which is decomposed into amplitude and phaseby ˜ h ( f ) = A ( f ) e − i Φ( f ) . The predicted EOS dependentfrequency cutoff in the waveform [29] is clearly shown inthe amplitude . Neutron stars that are more easily dis-rupted (larger Λ) result in an earlier and lower frequencydrop in their waveform amplitude than NS with smallerΛ. The phase Φ( f ) relative to the corresponding BBHwaveform has a much smoother behavior than the phaseof the time domain waveform. This feature will be use-ful in evaluating the Fisher matrix in Sec. VI. The noisethat is seen at frequencies above ∼ Tidal disruption occurs after the onset of mass shedding of theneutron star. The frequency at the onset of mass shedding isusually much lower than that of tidal disruption for BHNS bi-naries [20]. In Ref. [29], mass-shedding frequency was identifiedas the cutoff frequency but this underestimates the true cutofffrequency. See also Refs. [37, 38] for a discussion of dynamicalmass transfer.
V. HYBRID WAVEFORM CONSTRUCTION
Since our numerical simulations typically begin ∼ Q = 4 [41]. Although we have not ex-plored them in this context, other approaches have alsobeen taken for constructing phenomenological inspiral-merger-ringdown waveforms [42–46].For equal-mass BNS, Read et al. [12] compared thenumerical BNS waveform during inspiral to a point par-ticle post-Newtonian waveform. Specifically, they usedthe 3.5 post-Newtonian (TaylorT4) waveforms matchedon to the numerical waveforms to study the measurabilityof EOS parameters. They found that differences betweenthe analytic and numerical waveforms become apparent4 − ∼
20 GW cycles), where they find agreement with thesimulations to ± .
15 rad over the full simulation up tomerger [16].For the BHNS systems discussed here, we havematched the numerical waveforms to EOB waveformsthat include inspiral, merger, and ringdown phases in-stead of post-Newtonian waveforms which are often notreliable during the last few cycles for higher mass ratios.This choice also allows us to use longer matching win-dows that average over numerical noise and the effects ofeccentricity as shown in Sec. V B. We have chosen to usethe EOB formalism to generate inspiral-merger-ringdownwaveforms, although we note that other phenomenologi-cal waveforms would probably work. For simplicity, andbecause it appears that an accurate description of the lateinspiral dynamics just before merger requires 2PN tidalcorrections [14–16] which are not yet known, we will usethe EOB waveforms without tidal corrections. Our re-sults will therefore be lower limits on the measurability Q (cid:61) M NS (cid:61) M (cid:159) p.3 (cid:71) T I T F S I S F (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:72) ms (cid:76) h (cid:43) D (cid:144) M Q (cid:61) M NS (cid:61) M (cid:159) p.9 (cid:71) T I T F S I S F (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:72) ms (cid:76) h (cid:43) D (cid:144) M FIG. 2: h + and | h | = | h + − ih × | for BHNS waveforms for ( Q = 2 , M NS = 1 . M (cid:12) ) with two different EOS are representedby solid red and blue curves respectively. The softest EOS p.3Γ2.4 is on top and the stiffest EOS p.9Γ3.0 is on bottom. AnEOB BBH waveform (dashed) with the same values of Q and M NS is matched to each numerical waveform within the matchingwindow T I < t < T F bounded by solid vertical lines. A hybrid EOB BBH–Numerical BHNS waveform is generated by splicingthe waveforms together within a splicing window S I < t < S F bounded by dotted vertical lines. The matching window is 12 mslong and ends at the numerical merger time t NR M (time when the numerical waveform reaches its maximum amplitude), whilethe splicing window is 4 ms long and begins at the start of the matching window ( S I = T I ). of EOS parameters since the EOS dependence is comingsolely from the numerical waveforms. A. Matching procedure
We use a method similar to that described by Readet al. [12] to join each of the numerical BHNS wave-forms to a reference EOB waveform, generating a hy-brid EOB–numerical waveform. Denote a complex nu-merical waveform by h NR ( t ) = h NR+ ( t ) − ih NR × ( t ) and anEOB waveform with the same Q and M NS by h EOB ( t ) = h EOB+ ( t ) − ih EOB × ( t ). A constant time-shift τ and phase-shift Φ can be applied to the EOB waveform to matchit to a section of the numerical waveform by rewritingit as h EOB ( t − τ ) e − i Φ . We hold the numerical wave-form fixed because we must specify a matching window T I < t < T F , and as discussed below, there is only asmall region of the numerical waveforms over which avalid match can be performed. Once the values of τ and Φ are determined, we will then choose to instead hold the EOB waveform fixed and shift the numericalwaveform in the opposite direction by rewriting it as h shiftNR ( t ) = h NR ( t + τ ) e + i Φ . This is done so that all ofthe numerical waveforms with the same Q and M NS arealigned relative to a single fixed reference EOB waveform.Over a matching window T I < t < T F (bounded bysolid vertical lines in Fig. 2), the normalized match be-tween the waveforms is defined as m ( τ, Φ) = Re (cid:2) z ( τ ) e i Φ (cid:3) σ NR σ EOB ( τ ) , (11)where z ( τ ) = (cid:90) T F T I h NR ( t ) h ∗ EOB ( t − τ ) dt (12)and the normalizations for each waveform in the denom-enator are defined as σ = (cid:90) T F T I | h NR ( t ) | dt (13) (cid:45) (cid:45) (cid:45) (cid:72) ms (cid:76) (cid:200) h (cid:200) D (cid:144) M Q (cid:61) M NS (cid:61) M (cid:159) p (cid:32) . (cid:32) (cid:71) (cid:32) . E O B (cid:45) (cid:45) (cid:45) (cid:72) ms (cid:76) (cid:200) h (cid:200) D (cid:144) M Q (cid:61) M NS (cid:61) M (cid:159) p (cid:32) . (cid:32) (cid:71) (cid:32) . E O B (cid:45) (cid:45) (cid:45) (cid:72) ms (cid:76) (cid:200) h (cid:200) D (cid:144) M Q (cid:61) M NS (cid:61) M (cid:159) p (cid:32) . (cid:32) (cid:71) (cid:32) . E O B FIG. 3: Amplitude of the complex waveform h = h + − ih × .Dashed curves are EOB waveforms with the same Q and M NS .Matching and splicing conventions are those of Fig. 2. and σ ( τ ) = (cid:90) T F T I | h EOB ( t − τ ) | dt. (14)The time-shift τ and phase Φ are chosen to maximizethe match m ( τ, Φ) for a fixed matching window. Ex-plicitly, the phase is determined analytically to be Φ = − arg[ z ( τ )]; plugging this result back into Eq. (11), thetime-shift is given by maximizing | z ( τ ) | / [ σ NR σ EOB ( τ )]over τ . As stated above, once τ and Φ are found weshift the numerical waveform in the opposite direction togenerate h shiftNR ( t ) = h NR ( t + τ ) e + i Φ .A hybrid waveform is generated by smoothly turn-ing off the EOB waveform and smoothly turning on (cid:45) (cid:45) (cid:45) (cid:45) (cid:72) ms (cid:76) (cid:70) (cid:72) t (cid:76) (cid:45) (cid:70) E O B (cid:72) t (cid:76) Q (cid:61) M NS (cid:61) M (cid:159) p (cid:32) . (cid:32) (cid:71) (cid:32) . (cid:45) (cid:45) (cid:45) (cid:45) (cid:72) ms (cid:76) (cid:70) (cid:72) t (cid:76) (cid:45) (cid:70) E O B (cid:72) t (cid:76) Q (cid:61) M NS (cid:61) M (cid:159) p (cid:32) . (cid:32) (cid:71) (cid:32) . (cid:45) (cid:45) (cid:45) (cid:45) (cid:72) ms (cid:76) (cid:70) (cid:72) t (cid:76) (cid:45) (cid:70) E O B (cid:72) t (cid:76) Q (cid:61) M NS (cid:61) M (cid:159) p (cid:32) . (cid:32) (cid:71) (cid:32) . FIG. 4: Cumulative phase difference Φ − Φ EOB betweenBHNS waveform and EOB BBH waveform with the same Q and M NS . The phase is defined by breaking up each com-plex waveform into amplitude and cumulative phase h + ( t ) − ih × ( t ) = A ( t ) e − i Φ( t ) . The black point on each curve indi-cates the BHNS merger time t NR M defined as the time of max-imum amplitude A ( t NR M ). The curve is truncated when theamplitude AD/M drops below 0.01. Matching and splicingconventions are those of Fig. 2. the shifted numerical waveform over a splicing window S I < t < S F (bounded by dotted vertical lines in Fig. 2)which can be chosen independently of the matching win-dow. As in Ref. [12], we employ Hann windows w off ( t ) = 12 (cid:20) (cid:18) π [ t − S I ] S F − S I (cid:19)(cid:21) (15) w on ( t ) = 12 (cid:20) − cos (cid:18) π [ t − S I ] S F − S I (cid:19)(cid:21) . (16) broadbandET (cid:45) D500 1000 2000 3000 50001 (cid:180) (cid:45) (cid:180) (cid:45) (cid:180) (cid:45) (cid:180) (cid:45) (cid:180) (cid:45) (cid:180) (cid:45) f gw (cid:72) Hz (cid:76) f (cid:144) (cid:200) h (cid:142) (cid:200) (cid:72) H z (cid:45) (cid:144) (cid:76) M (cid:87) f (cid:144) (cid:200) h (cid:142) (cid:200) D (cid:144) M (cid:72) H z (cid:45) (cid:144) (cid:76) Q (cid:61) M NS (cid:61) M (cid:159) p (cid:32) . (cid:32) (cid:71) (cid:32) . E O B broadbandET (cid:45) D500 1000 2000 3000 50001 (cid:180) (cid:45) (cid:180) (cid:45) (cid:180) (cid:45) (cid:180) (cid:45) (cid:180) (cid:45) (cid:180) (cid:45) f gw (cid:72) Hz (cid:76) f (cid:144) (cid:200) h (cid:142) (cid:200) (cid:72) H z (cid:45) (cid:144) (cid:76) M (cid:87) f (cid:144) (cid:200) h (cid:142) (cid:200) D (cid:144) M (cid:72) H z (cid:45) (cid:144) (cid:76) Q (cid:61) M NS (cid:61) M (cid:159) p (cid:32) . (cid:32) (cid:71) (cid:32) . E O B broadbandET (cid:45) D500 1000 2000 3000 50001 (cid:180) (cid:45) (cid:180) (cid:45) (cid:180) (cid:45) (cid:180) (cid:45) (cid:180) (cid:45) (cid:180) (cid:45) f gw (cid:72) Hz (cid:76) f (cid:144) (cid:200) h (cid:142) (cid:200) (cid:72) H z (cid:45) (cid:144) (cid:76) M (cid:87) f (cid:144) (cid:200) h (cid:142) (cid:200) D (cid:144) M (cid:72) H z (cid:45) (cid:144) (cid:76) Q (cid:61) M NS (cid:61) M (cid:159) p (cid:32) . (cid:32) (cid:71) (cid:32) . E O B FIG. 5: Weighted Fourier transform 2 f / | ˜ h ( f ) | of numer-ical waveforms where ˜ h = (˜ h + + ˜ h × ). Dot-dashed curvesare EOB waveforms with the same Q and M NS . The left axisis scaled to a distance of 100 Mpc, and the noise S / n ( f )for broadband aLIGO and ET-D are shown for compari-son. In each plot the numerical waveform monotonically ap-proaches the EOB waveform as the tidal parameter Λ de-creases. Matching and splicing conventions are those of Fig. 2. (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) f gw (cid:72) Hz (cid:76) (cid:70) (cid:72) f (cid:76) (cid:45) (cid:70) E O B (cid:72) f (cid:76) Q (cid:61) M NS (cid:61) M (cid:159) p (cid:32) . (cid:32) (cid:71) (cid:32) . p (cid:32) . (cid:32) (cid:71) (cid:32) . (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) f gw (cid:72) Hz (cid:76) (cid:70) (cid:72) f (cid:76) (cid:45) (cid:70) E O B (cid:72) f (cid:76) Q (cid:61) M NS (cid:61) M (cid:159) p (cid:32) . (cid:32) (cid:71) (cid:32) . p (cid:32) . (cid:32) (cid:71) (cid:32) . (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) f gw (cid:72) Hz (cid:76) (cid:70) (cid:72) f (cid:76) (cid:45) (cid:70) E O B (cid:72) f (cid:76) Q (cid:61) M NS (cid:61) M (cid:159) p (cid:32) . (cid:32) (cid:71) (cid:32) . p (cid:32) . (cid:32) (cid:71) (cid:32) . FIG. 6: Cumulative phase difference Φ − Φ EOB of the Fouriertransform between BHNS waveform and EOB BBH wave-form of the same mass and mass ratio. The phase is definedby breaking up the Fourier transform ˜ h = (˜ h + + ˜ h × ) ofeach waveform into amplitude and cumulative phase ˜ h ( f ) = A ( f ) e − i Φ( f ) . Matching and splicing conventions are those ofFig. 2. The hybrid waveform is then written h hybrid ( t ) = h EOB ( t ) t < S I w off ( t ) h EOB ( t ) + w on ( t ) h shiftNR ( t ) S I < t < S F h shiftNR ( t ) t > S F . (17)As shown in Fig. 2, we choose the start of the splicinginterval to be the same as the start of the matching win-dow S I = T I and choose the end of the splicing windowto be S F = T I + 4 ms. It is also necessary to use thesewindows to smoothly turn on the hybrid waveform at lowfrequency when performing a discrete Fourier transformto avoid the Gibbs phenomenon. Unlike the case for BNSwaveforms, it is not necessary to window the end of thehybrid waveform as the amplitude rapidly decays to zeroanyway during the ringdown.For concreteness we define t = 0 as the EOB BBHmerger time t EOB M when the EOB waveform reaches itsmaximum amplitude. After matching to the EOB wave-form, the time when the numerical BHNS waveformreaches its maximum amplitude is t NR M . B. Dependence on matching window
Because the numerical BHNS waveforms are close butnot identical to the EOB BBH waveform during the in-spiral and because there is some noise in the BHNS wave-forms, the time shift that maximizes the match dependson the choice of matching window. The matching windowshould exclude the first couple of cycles of the numeri-cal waveform during which time the simulation is settlingdown from the initial conditions. It should also excludethe merger/ringdown which are strongly dependent onthe presence of matter. The window must also be wideenough to average over numerical noise and, as we shallsee below, the effects of eccentricity in the simulations.The numerical merger time t NR M relative to the EOBBBH merger time t EOB M as a function of the end of thematching window T F − t NR M provides a useful diagnos-tic of the matching procedure. Results for matching two Q = 2 , M NS = 1 . M (cid:12) waveforms with different equa-tions of state to an EOB waveform are shown in Fig. 7.The horizontal axis is the end time T F of the matchingwindow relative to the numerical merger time t NR M . Fornegative values, the matching window contains the BHNSinspiral only. For positive values, the matching windowalso contains part of the BHNS ringdown. The verticalaxis is the location of the shifted numerical merger time t NR M after finding the best match. Four different windowdurations ∆ t = T F − T I are shown. The drift in the bestfit merger time t NR M most likely arises from the neglect oftidal effects in the EOB waveform which lead to an accu-mulating phase shift in the waveform, although it couldalso arise from numerical angular momentum loss fromfinite resolution of the simulations. Further work is inprogress to understand this issue [13].When the matching window duration is of order oneorbital period or shorter, the time-shift oscillates as afunction of T F − t NR M . We attribute this effect to theeccentricity in the numerical waveform that results frominitial data with no radial velocity. For larger matching-window durations, the effect of eccentricity is averagedout.To demonstrate concretely that the decaying oscilla-tions for ∆ t = 4 ms are the result of eccentricity, wematched an EOB BBH waveform with eccentricity to the (cid:68) t (cid:61) (cid:68) t (cid:61) (cid:68) t (cid:61)
12 ms (cid:68) t (cid:61)
16 ms (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) T F (cid:45) t M NR (cid:72) ms (cid:76) t M NR (cid:72) m s (cid:76) Q (cid:61) M NS (cid:61) M (cid:159) p (cid:32) .3 (cid:32)(cid:71) (cid:32) .4 (cid:68) t (cid:61) (cid:68) t (cid:61) (cid:68) t (cid:61)
12 ms (cid:68) t (cid:61)
16 ms (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) T F (cid:45) t M NR (cid:72) ms (cid:76) t M NR (cid:72) m s (cid:76) Q (cid:61) M NS (cid:61) M (cid:159) p (cid:32) .9 (cid:32)(cid:71) (cid:32) .0 FIG. 7: Dependence of time shift between numerical andEOB waveform on the end time T F − t NR M and width ∆ t ofthe matching window. Q = 2 and M NS = 1 . M (cid:12) for eachwaveform. The EOS used are p.3Γ2.4 (top panel), and p.9Γ3.0(bottom panel). The EOB waveform has zero eccentricity. equivalent zero eccentricity EOB BBH waveform. EOBwaveforms can be generated with small eccentricity bystarting the EOB equations of motion with quasicircular(zero radial velocity) initial conditions late in the inspi-ral. The result is shown in Fig. 8 for an EOB waveformwith the same quasicircular initial conditions as the sim-ulation for the EOS p.3Γ2.4 shown in Fig. 7. The oscil-lations take exactly the form of those shown in Fig. 7,except without the drift and offset.We estimate that the initial eccentricities in the sim-ulations used in this paper are e ∼ .
03. Decreasingthe initial eccentricity by about an order of magnitude,possibly using an iterative method that adjusts the ini-tial radial velocity [47], will remove this issue and allowone to determine the phase shift due to tidal interactionsduring the inspiral part of the simulation.
VI. PARAMETER ESTIMATION
The output of a gravitational-wave detector s ( t ) = n ( t ) + h ( t ) is the sum of detector noise n ( t ) and a possi-ble gravitational-wave signal h ( t ). Stationary, Gaussiannoise is characterized by its power spectral density (PSD) S n ( | f | ) defined by (cid:104) ˜ n ( f )˜ n ∗ ( f (cid:48) ) (cid:105) = 12 δ ( f − f (cid:48) ) S n ( | f | ) . (18)0 (cid:68) t (cid:61) (cid:68) t (cid:61) (cid:68) t (cid:61)
12 ms (cid:68) t (cid:61)
16 ms (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) T F (cid:45) t M NR (cid:72) ms (cid:76) t M NR (cid:72) m s (cid:76) Q (cid:61) M NS (cid:61) M (cid:159) M (cid:87) (cid:61) FIG. 8: Same as Fig. 7, but matching an eccentric EOB BBHwaveform with the quasicircular initial condition M Ω =0 .
028 to a zero eccentricity EOB BBH waveform.
The gravitational wave signal is given in terms of the twopolarizations of the gravitational wave by h ( t ) = F + h + ( t ) + F × h × ( t ) , (19)where F + , × are the detector response functions and de-pend on the location of the binary and the polarizationangle of the waves. We assume the binary is optimallylocated at the zenith of the detector and optimally ori-ented with its orbital plane parallel to that of the detec-tor. This condition is equivalent to averaging h + and h × ( F + = F × = 1 / h ( t ) in additive Gaussian noiseis ρ = ( h | s ) (cid:112) ( h | h ) (20)where the inner product between two signals h and h is given by ( h | h ) = 4Re (cid:90) ∞ ˜ h ( f )˜ h ∗ ( f ) S n ( f ) df. (21)In searches for gravitational-wave signals from compactbinary mergers, a parametrized signal h ( t ; θ A ) is knownin advance of detection, and the parameters θ A must beestimated from the measured detector output s ( t ). Theparameters θ A of an inspiral are estimated by maximiz-ing the inner product of the signal s ( t ) over the templatewaveforms h ( t ; θ A ). In the high signal-to-noise limit, thestatistical uncertainty in the estimated parameters ˆ θ A arising from the instrumental noise can be estimated us-ing the Fisher matrixΓ AB = (cid:18) ∂h∂θ A (cid:12)(cid:12)(cid:12)(cid:12) ∂h∂θ B (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) ˆ θ A . (22)Note that ˆ θ A are the parameter values that maximize thesignal-to-noise. The variance σ A = σ AA = (cid:104) (∆ θ A ) (cid:105) and covariance σ AB = (cid:104) ∆ θ A ∆ θ B (cid:105) of the parameters are thengiven in terms of the Fisher matrix by (cid:104) ∆ θ A ∆ θ B (cid:105) = (Γ − ) AB . (23)For hybrid waveforms, the partial derivatives in theFisher matrix must be approximated with finite differ-ences. It is most robust to compute the derivatives ofthe Fourier transforms used in the inner product. Werewrite the Fourier transform of each waveform in termsof the amplitude A and phase Φ as exp[ln A − i Φ] as givenin Eq. (A9). The derivatives ∂ ln A/∂θ A and ∂ Φ /∂θ A arethen evaluated with finite differencing. More details ofthis and the other methods we tested are given in Ap-pendix A.In general, errors in the parameters θ A are correlatedwith each other forming an error ellipsoid in parameterspace determined by the Fisher matrix Γ AB . The un-correlated parameters that are best extracted from thesignal are found by diagonalizing Γ AB . These new pa-rameters are linear combinations of the original parame-ters θ A . We focus attention below on the two parameterslog( p ) and Γ, and fix all other parameters as follows. Weuse the masses and spins determined from the numeri-cal simulations and fix the time and phase shifts as de-termined during the hybrid waveform construction. Wetherefore construct the error ellipses in { log( p ) , Γ } pa-rameter space and identify the parameter with the small-est statistical errors. We will leave an analysis of correla-tions due to uncertainty in masses and BH spin to futurework. A. Broadband aLIGO and ET
For the BHNS systems discussed here, the greatest de-parture from BBH behavior occurs for gravitational-wavefrequencies in the range 500–5000 Hz. As a result, de-tector configurations optimized for detection of BHNSsystems with low noise in the region below 500 Hz arenot ideal for estimating EOS parameters. We thereforepresent results for the broadband aLIGO noise curve [49]and the ET-D noise curve [50] shown in Fig. 9. Thebroadband aLIGO configuration uses zero-detuning ofthe signal recycling mirror and a high laser power, re-sulting in significantly lower noise above 500 Hz at theexpense of slightly higher noise at lower frequencies. Sev-eral noise curves have been considered for the EinsteinTelescope denoted ET-B [51], ET-C [52], and ET-D [50].We will use the most recent ET-D configuration, andnote that in the 500–5000 Hz range all of the ET config-urations have a similar sensitivity. The published noisecurves, and those used in this paper, are for a single in-terferometer of 10 km with a 90 ◦ opening angle. Thecurrent ET proposal is to have three individual interfer-ometers each with a 60 ◦ opening angle. This will shiftthe noise curve down appoximately 20% [50].In Figs. 10 and 11, we show the resulting 1- σ error el-lipses in the 2-dimensional parameter space { log( p ) , Γ } narrowbandbroadbandET (cid:45) D10 50 100 500 1000 50001 (cid:180) (cid:45) (cid:180) (cid:45) (cid:180) (cid:45) (cid:180) (cid:45) (cid:180) (cid:45) (cid:180) (cid:45) (cid:180) (cid:45) f (cid:72) Hz (cid:76) S h (cid:144) (cid:72) H z (cid:45) (cid:144) (cid:76) f R f R f R FIG. 9: Noise PSD for broadband aLIGO (dashed blue),ET-D (dotted red) and various configurations of narrowbandaLIGO (solid black). The minima of the narrowband config-uration are labeled f R for an optimally oriented BHNS with Q = 2 and M NS =1 . M (cid:12) at a distance of 100 Mpc. Surfaces of constantΛ / and NS radius, which are almost parallel to eachother, are also shown. One can see that the error ellipsesare aligned with these surfaces. This indicates that, asexpected, Λ / is the parameter that is best extractedfrom BHNS gravitational-wave observations. BecauseΛ / and R are so closely aligned we will use these twoparameters interchangeably. (cid:180)(cid:180)(cid:180)(cid:180) (cid:71) l og (cid:72) p (cid:76) (cid:64) d y ne (cid:144) c m (cid:68) (cid:76) (cid:144) (cid:61) . (cid:76) (cid:144) (cid:61) . (cid:76) (cid:144) (cid:61) . FIG. 10: Two 1– σ error ellipses for broadband aLIGO.Evenly spaced contours of constant Λ / are also shown. Eachellipse is centered on the estimated parameter ˆ θ A denoted bya × . The semimajor axes are significantly longer than thewidth of the figure, so each ellipse appears as a pair of par-allel lines. Matching and splicing conventions are those ofFig. 2. As mentioned above, there is some freedom in con-struction of the hybrid waveforms. The size and orien- (cid:180)(cid:180)(cid:180)(cid:180)(cid:180)(cid:180)(cid:180)(cid:180)(cid:180)(cid:180)(cid:180)(cid:180)(cid:180)(cid:180) (cid:71) l og (cid:72) p (cid:76) (cid:64) d y ne (cid:144) c m (cid:68) (cid:76) (cid:144) (cid:61) . (cid:76) (cid:144) (cid:61) . (cid:76) (cid:144) (cid:61) . (cid:180)(cid:180)(cid:180)(cid:180)(cid:180)(cid:180)(cid:180)(cid:180)(cid:180)(cid:180)(cid:180)(cid:180)(cid:180)(cid:180) (cid:71) l og (cid:72) p (cid:76) (cid:64) d y ne (cid:144) c m (cid:68) R (cid:61) (cid:32) k m R (cid:61) (cid:32) k m R (cid:61) (cid:32) k m R (cid:61) (cid:32) k m FIG. 11: 1– σ error ellipses for ET-D. Evenly spaced con-tours of constant Λ / ( R ) are also shown on top (bottom).Matching and splicing conventions are those of Fig. 2. tation of the error ellipses also depend on the details ofthis construction. We find that as long as the matchingwindow is longer than approximately four gravitational-wave cycles to average out the effects of eccentricity anddoes not include the first two gravitational-wave cycles,the orientation of the error ellipses does not change sig-nificantly. As expected, the size of the ellipses decreasesas more of the numerical waveform is incorporated intothe hybrid waveform. We therefore adopt the last 12 msbefore merger of each numerical waveform as the match-ing window and the first 4 ms of the matching windowfor splicing as shown in Fig. 2.We have emphasized that, to within present numer-ical accuracy, the late-inspiral waveform is determinedby the single parameter Λ / . This implies that, by us-ing countours of constant Λ in the EOS space, one couldhave obtained the constraint on the EOS, summarized inFigs. 10 and 11 by varying only a single EOS parameter.For the simulations with other mass ratios and neutronstar masses, we have used as our single parameter log( p )2and not Γ because countours of constant p more closelycoincide with contours of constant Λ and because Λ isa one to one function of log( p ) throughout the param-eter space. The one-parameter Fisher matrix can thenbe evaluated with finite differencing using the waveformsand values of Λ at two points in EOS parameter spacewith different log( p ).The uncertainties in Λ / and R are shown in Figs. 12and 13 for broadband aLIGO and for ET respectively.The uncertainty in these quantities is ∼ ∼ Q = 3 are somewhatlarger than for Q = 2, but not significantly so. It is notclear how rapidly the uncertainty in Λ / and R will in-crease as the mass ratio is increased toward more realisticvalues. On the one hand the tidal distortion is likely tobe much smaller for larger Q . On the other hand theoverall signal will be louder, and the merger and ring-down will occur at lower frequencies where the noise islower. Additional simulations for higher Q are needed toaddress this question. (cid:242) (cid:242) (cid:242)(cid:225) (cid:225) (cid:225)(cid:225)(cid:225)(cid:225)(cid:225)(cid:225) (cid:225)(cid:225)(cid:225)(cid:225) (cid:225)(cid:225)(cid:225)(cid:225)(cid:248) (cid:248) (cid:248) (cid:242) Q (cid:61) M NS (cid:61) M (cid:159) (cid:225) Q (cid:61) M NS (cid:61) M (cid:159) (cid:248) Q (cid:61) M NS (cid:61) M (cid:159) (cid:76) (cid:144) Σ (cid:76) (cid:144) (cid:242) (cid:242) (cid:242)(cid:225) (cid:225)(cid:225)(cid:225) (cid:225)(cid:225)(cid:225)(cid:225) (cid:225)(cid:225)(cid:225)(cid:225) (cid:225)(cid:225)(cid:225)(cid:225)(cid:248) (cid:248) (cid:248) (cid:242) Q (cid:61) M NS (cid:61) M (cid:159) (cid:225) Q (cid:61) M NS (cid:61) M (cid:159) (cid:248) Q (cid:61) M NS (cid:61) M (cid:159)
10 11 12 13 14 15 16012345 R (cid:72) km (cid:76) Σ R (cid:72) k m (cid:76) FIG. 12: 1- σ uncertainty σ Λ / and σ R as a function of theparameters Λ / or R for the broadband aLIGO noise PSD.Matching and splicing conventions are those of Fig. 2. (cid:242) (cid:242) (cid:242)(cid:225) (cid:225) (cid:225)(cid:225)(cid:225)(cid:225)(cid:225)(cid:225) (cid:225)(cid:225)(cid:225)(cid:225) (cid:225)(cid:225)(cid:225)(cid:225)(cid:248) (cid:248) (cid:248) (cid:242) Q (cid:61) M NS (cid:61) M (cid:159) (cid:225) Q (cid:61) M NS (cid:61) M (cid:159) (cid:248) Q (cid:61) M NS (cid:61) M (cid:159) (cid:76) (cid:144) Σ (cid:76) (cid:144) (cid:242) (cid:242) (cid:242)(cid:225) (cid:225)(cid:225)(cid:225) (cid:225)(cid:225)(cid:225)(cid:225) (cid:225)(cid:225)(cid:225)(cid:225) (cid:225)(cid:225)(cid:225)(cid:225)(cid:248) (cid:248) (cid:248) (cid:242) Q (cid:61) M NS (cid:61) M (cid:159) (cid:225) Q (cid:61) M NS (cid:61) M (cid:159) (cid:248) Q (cid:61) M NS (cid:61) M (cid:159)
10 11 12 13 14 15 160.00.10.20.30.40.5 R (cid:72) km (cid:76) Σ R (cid:72) k m (cid:76) FIG. 13: Same as Fig. 12, but with the ET-D noise PSD.Error estimates are an order of magnitude smaller than forbroadband aLIGO.
B. Narrowband aLIGO
The presence of a signal-recycling cavity in the aLIGOinstruments will allow them to be tuned to have improvednarrowband sensitivity at the expense of bandwith. Twoparameters control the narrowband capabilities of the in-struments [53–55]: the signal recycling mirror transmis-sivity effectively sets the frequency bandwidth of the in-strument, while the length of the signal recycling cavity(or equivalently the signal-recycling cavity tuning phase)controls the central frequency f R of the best sensitivity.By tuning one or more of the aLIGO detectors to oper-ate in narrowband mode, it may be possible to improveestimates of the EOS parameters.We have examined several narrowband tunings withcentral frequencies that vary between approximately f R = 500 Hz and 4000 Hz. These noise curves use a sig-nal recycling mirror transmissivity of 0.011 and a signal-recycling cavity tuning phase ranging from 10 ◦ down to1 ◦ , and were generated using the program gwinc [56].Three of these noise curves are shown in Fig. 9. In Fig. 14we plot the 1– σ uncertainty in NS radius σ R as a functionof the narrowband central frequency f R . For the wave-forms considered in this paper the optimal narrowbandfrequency is in the range 1000 Hz (cid:46) f R (cid:46) f R happens to be tuned towithin a few hundred Hz of the minimum for that BHNSevent. In Ref. [57], Hughes discussed a method for de-termining the best frequency f R to tune a narrowbandeddetector to extract an EOS dependent cutoff frequencyfrom a sequence of identical BNS inspirals. While thistechnique is not directly applicable to BHNS systems,which have different masses and spins, a similar approachcould be used to combine multiple BHNS observations. R (cid:61) (cid:72) p.35 (cid:71) (cid:76) R (cid:61) (cid:72) p.45 (cid:71) (cid:76) R (cid:61) (cid:72) p.7 (cid:71) (cid:76) f R (cid:72) Hz (cid:76) Σ R (cid:72) k m (cid:76) Q (cid:61) M NS (cid:61) M (cid:159) broad ET (cid:45) D R (cid:61) (cid:72) p.35 (cid:71) (cid:76) R (cid:61) (cid:72) p.45 (cid:71) (cid:76) R (cid:61) (cid:72) p.65 (cid:71) (cid:76) R (cid:61) (cid:72) p.8 (cid:71) (cid:76) f R (cid:72) Hz (cid:76) Σ R (cid:72) k m (cid:76) Q (cid:61) M NS (cid:61) M (cid:159) broad ET (cid:45) D R (cid:61) (cid:72) p.4 (cid:71) (cid:76) R (cid:61) (cid:72) p.6 (cid:71) (cid:76) R (cid:61) (cid:72) p.8 (cid:71) (cid:76) f R (cid:72) Hz (cid:76) Σ R (cid:72) k m (cid:76) Q (cid:61) M NS (cid:61) M (cid:159) broad ET (cid:45) D FIG. 14: 1– σ uncertainty in R for different configurations ofnarrowband aLIGO and for different EOS. f R defines the fre-quency where S n is a minimum as shown in Fig. 9. Horizontallines on the left and right indicate the corresponding 1– σ er-rors for broadband aLIGO and ET-D respectively. Matchingand splicing conventions are those of Fig. 2. VII. DISCUSSIONA. Results
Using a large set of simulations incorporating a two-parameter EOS, we have found that the tidal deforma-bility Λ / , or equivalently the NS radius R , is the pa-rameter that will be best extracted from BHNS wave-forms. These parameters can be estimated to 10–40%with broadband aLIGO for an optimally oriented BHNSbinary at 100 Mpc. The narrowband aLIGO configu-ration can do slightly better if it is tuned to within afew hundred Hz of the ideal frequency for a given BHNSevent. The proposed Einstein Telescope will have anorder-of-magnitude better sensitivity to the EOS param-eters.Although we have used a particular EOS parametriza-tion to show that Λ is the parameter that is observedduring BHNS coalescence, this result can be used to con-strain any EOS model—an EOS based on fundamentalnuclear theory in addition to a parametrized phenomeno-logical EOS. In particular, several parametrizations haverecently been developed, including a spectral represen-tation [58], a reparametrization of the piecewise poly-trope [59], and a generalization that also includes nuclearparameters [60].The results presented here can be compared with re-cent work to determine the mass and radius of individualNS in Type-1 X-ray bursts. ¨Ozel et al. [61] have obtainedmass and radius measurements from several systems bysimultaneously measuring the flux F , which is likely closeto the Eddington value, and the blackbody temperature T during X-ray bursts of systems with accurately deter-mined distances. During the burst, the emission areaof the photosphere F/ ( σT ) expands, contracts, thenreaches a constant value, and ¨Ozel et al. have arguedthat the final area corresponds to that of the NS surface.They obtain estimates of NS mass and radii with O (10%)1– σ uncertainty. Steiner et al. [60] have also consideredthese systems, but argue that the final emission area doesnot necessarily correspond to that of the NS surface, andas a result obtain slightly smaller NS radii and larger un-certainties in the mass and radius. These radius errorestimates are slightly smaller than those for the BHNSsystems we have considered at 100 Mpc. However, wenote that binary inspiral observations are subject to lesssystematic uncertainty due to questions of composition ofthe photosphere and associating it with the NS surface.The uncertainty in NS radius for the merger and ring-down of BHNS systems examined here is of roughlythe same size as that found for the last few orbits upto merger of BNS systems at the same 100 Mpc dis-tance [12, 13]. BNS inspirals, however, will likely oc-cur more frequently, and, including a tidally correctedinspiral–numerical hybrid, BNS systems are likely to haveuncertainties that are smaller than BHNS systems by afactor of a few. Considering the post-merger phase forBNS waveforms may also provide additional information.4Expected NS masses in both BNS and BHNS systems areslightly smaller than those measured for X-ray bursterswhich have accreted matter from their companion, soBNS and BHNS GW observations may complement X-ray burst observations by better constraining the lowerdensity range of the EOS which is not well constainedfrom X-ray burst observations [59, 61]. B. Future work
We have used in this paper several simplifications andconventions which can significantly effect the accuracy towhich EOS parameters can be extracted. We list thembelow and describe how changing them would effect theparameter error estimates.1.
Finite length of numerical waveforms
The BHNS waveforms used here include only thelast ∼
10 GW cycles of inspiral as well as the mergerand ringdown, of which the first few cycles of in-spiral are unreliable due to inexact initial data.Matching the numerical waveform to a tidally cor-rected inspiral waveform instead of just the point-particle waveform will increase the overall depar-ture from point-particle behavior by (i) creating aphase shift during the early inspiral, and more im-portantly (ii) adding to the phase of the late inspi-ral and merger the accumulated phase shift fromthe early inspiral – a phase shift that is not alreadyincluded in the stronger signal of the late inspiral.The tidal corrections are now known up to 1PN or-der. For simulations with the current number oforbits, however, it appears that higher order tidalcorrections will be needed to fully describe the lateinspiral where matching to numerical waveforms isdone. We leave the issue of generating tidally cor-rected inspiral-numerical hybrid waveforms to fu-ture work.2.
Event rates
Estimates of the detectability of EOS parametersin BNS systems are often given for an event at adistance of 100 Mpc, and we have used the sameconvention here to state the results above. Therelevant event rate is, therefore, the expected num-ber of detected events that will have an effectivedistance D eff ≤
100 Mpc. (The effective distance D eff depends on the location of the binary andits inclination relative to the detector. For anoptimally oriented and located binary, one finds D = D eff while typically D ≤ D eff .) The aLIGOinspiral rates for BNS systems are highly uncertainwith { low, most likely, high } estimates of { } Mpc − Myr − [1] or { } yr − witheffective distance D eff ≤
100 Mpc. Rates are evenmore uncertain for BHNS systems with rate esti-mates of { } yr − with effective dis-tance D eff ≤
100 Mpc [1]. Since the uncertainty in EOS parameters scales linearly with distance[ σ Λ / = σ Λ / , ( D/ D , the estimated detection rates ofsystems with effective distance D eff ≤
400 Mpc are { } yr − with a four-fold increase in un-certainty of Λ / . Fortunately, for N obs identicalevents and N det identical detectors, the uncertaintyalso scales as σ Λ / / √ N obs N det .3. Expected NS masses and mass ratios
The simulations we used included realistic massneutron stars of 1.2 and 1.35 M (cid:12) . On the otherhand, black hole masses are expected to be manytimes larger [3], with likely mass ratios closer to Q ∼ M (cid:12) –1.4 M (cid:12) system)than the Q = 2 and 3 systems we examined here.Additional simulations for mass ratios of 4 and 5are in progress.4. Spinning BH
In this paper we have not examined the effect of aspinning BH. The analytic results of Ref. [11] indi-cated that spin does not significantly improve thesensitivity to Λ for the inspiral up to the pointof tidal disruption. However, numerical simula-tions [24, 27, 28] have shown that spin can stronglyaffect the dynamics near tidal disruption and theamount of matter left over in an accretion disk. Wehave performed several tens of simulations of non-precessing BHNS systems with spinning BH withvarious BH spins, mass ratios, NS masses, and EOSparameters, and an analysis of how BH spin affectsthe detectability of EOS parameters will be the sub-ject of the next paper.5.
Correlations between parameters
In our Fisher analysis we have assumed that themass ratio, NS mass, and BH spin will be deter-mined to sufficient accuracy during the inspiral toseparate them from EOS effects during the mergerand ringdown. A full Fisher analysis using all of theBHNS parameters should be done to find the extentto which uncertainties in the other parameters altermeasurability estimates of the EOS parameters.Because BHNS waveforms smoothly deviate from cor-responding BBH waveforms as Λ increases, it is likelythat one can find a good analytical approximation for thefull inspiral, merger, and ringdown waveform by modify-ing analytical BBH waveforms. Accurate waveforms fornon-spinning BBH systems using the EOB approach havebeen developed [41, 63] and work to find EOB waveformsfor spinning BBH systems is in progress [64, 65]. Tidalinteractions have also been incorporated into the EOBapproach for BNS systems with good agreement withthe inspiral waveform from numerical simulations whenparametrized 2PN tidal interactions are fit to the numer-ical waveform [15, 16]. Another approach is to use phe-nomenological waveforms that fit the frequency domain5post-Newtonian inspiral waveform to a phenomenologicalmerger and ringdown for both spinning and non-spinningBBH systems [42]. Both of these approaches may workfor generating full analytic BHNS waveforms as well. Acomplete description of the BHNS waveform will likelyinclude corrections for the l = 3 tidal field and otherhigher order corrections. However, it is not clear giventhe current set of simulations that these effects wouldbe observable with either aLIGO or a third generationdetector such as ET. Acknowledgments
We thank Alessandro Nagar for significant help withunderstanding the EOB formalism, Yi Pan and Alessan- dra Buonanno for providing EOB waveforms with spin-ning black holes, Jocelyn Read for providing routinesused in the data analysis, and Jolien Creighton for gen-erating narrowband noise curves. This work was sup-ported by NSF Grants PHY-1001515 and PHY-0923409,by Grant-in-Aid for Scientific Research (21340051), byGrant-in-Aid for Scientific Research on Innovative Area(20105004) of Japanese MEXT, and by a Grant-in-Aidof JSPS. BL would also like to acknowledge support froma UWM Graduate School Fellowship and the WisconsinSpace Grant Consortium. KK is also supported by theGrant-in-Aid for the Global COE Program “The NextGeneration of Physics, Spun from Universality and Emer-gence” of Japanese MEXT.” [1] J. Abadie, B. P. Abbott, R. Abbott, M. Abernathy,T. Accadia, F. Acernese, C. Adams, R. Adhikari,P. Ajith, B. Allen, et al., Classical and Quantum Gravity , 173001 (2010), 1003.2480.[2] K. Belczynski, M. Dominik, T. Bulik, R. O’Shaughnessy,C. Fryer, and D. E. Holz, Astrophys. J. Lett. , L138(2010), 1004.0386.[3] K. Belczynski, V. Kalogera, and T. Bulik, Astrophys. J. , 407 (2002), arXiv:astro-ph/0111452.[4] T. Hinderer, Astrophys. J. , 1216 (2008), 0711.2420.[5] ´E. ´E. Flanagan and T. Hinderer, Phys. Rev. D ,021502 (2008), 0709.1915.[6] J. Vines, T. Hinderer, and ´E. ´E. Flanagan, ArXiv e-prints(2011), 1101.1673.[7] J. Vines and ´E. ´E. Flanagan, ArXiv e-prints (2010),1009.4919.[8] T. Binnington and E. Poisson, Phys. Rev. D , 084018(2009), 0906.1366.[9] T. Damour and A. Nagar, Phys. Rev. D , 084035(2009), 0906.0096.[10] T. Hinderer, B. D. Lackey, R. N. Lang, and J. S. Read,Phys. Rev. D , 123016 (2010), 0911.3535.[11] F. Pannarale, L. Rezzolla, F. Ohme, and J. S. Read,ArXiv e-prints (2011), 1103.3526.[12] J. S. Read, C. Markakis, M. Shibata, K. Ury¯u, J. D. E.Creighton, and J. L. Friedman, Phys. Rev. D , 124033(2009), 0901.3258.[13] L. Baiotti, J. Creighton, B. Giacomazzo, K. Kyu-toku, C. Markakis, J. Read, L. Rezzolla, M. Shibata,K. Taniguchi, and J. Friedman, in progress (2011).[14] T. Damour and A. Nagar, Phys. Rev. D , 084016(2010), 0911.5041.[15] L. Baiotti, T. Damour, B. Giacomazzo, A. Nagar, andL. Rezzolla, Phys. Rev. Lett. , 261101 (2010).[16] L. Baiotti, T. Damour, B. Giacomazzo, A. Nagar, andL. Rezzolla, ArXiv e-prints (2011), 1103.3874.[17] K. Kyutoku, M. Shibata, and K. Taniguchi, Phys. Rev.D , 044049 (2010), 1008.1460.[18] M. Shibata and K. Ury¯u, Phys. Rev. D , 121503 (2006).[19] M. Shibata and K. Uryu, Classical and Quantum Gravity , 125 (2007), arXiv:astro-ph/0611522.[20] M. Shibata and K. Taniguchi, Phys. Rev. D , 084015 (2008).[21] M. Shibata, K. Kyutoku, T. Yamamoto, andK. Taniguchi, Phys. Rev. D , 044030 (2009),0902.0416.[22] M. Shibata and K. Kyutoku, Progress of TheoreticalPhysics Supplement , 17 (2010).[23] Z. B. Etienne, J. A. Faber, Y. T. Liu, S. L. Shapiro,K. Taniguchi, and T. W. Baumgarte, Phys. Rev. D ,084002 (2008), 0712.2460.[24] Z. B. Etienne, Y. T. Liu, S. L. Shapiro, and T. W. Baum-garte, Phys. Rev. D , 044024 (2009), 0812.2245.[25] M. D. Duez, F. Foucart, L. E. Kidder, H. P. Pfeiffer,M. A. Scheel, and S. A. Teukolsky, Phys. Rev. D ,104015 (2008), 0809.0002.[26] M. D. Duez, F. Foucart, L. E. Kidder, C. D. Ott, andS. A. Teukolsky, Classical and Quantum Gravity ,114106 (2010), 0912.3528.[27] F. Foucart, M. D. Duez, L. E. Kidder, and S. A. Teukol-sky, Phys. Rev. D , 024005 (2011), 1007.4203.[28] K. Kyutoku, H. Okawa, M. Shibata, and K. Taniguchi,ArXiv e-prints (2011), 1108.1189.[29] M. Vallisneri, Physical Review Letters , 3519 (2000),arXiv:gr-qc/9912026.[30] J. S. Read, B. D. Lackey, B. J. Owen, and J. L. Friedman,Phys. Rev. D , 124032 (2009), 0812.2163.[31] J. M. Lattimer and M. Prakash, Astrophys. J. , 426(2001), arXiv:astro-ph/0002232.[32] P. B. Demorest, T. Pennucci, S. M. Ransom, M. S. E.Roberts, and J. W. T. Hessels, Nature (London) ,1081 (2010), 1010.5788.[33] K. Kyutoku, M. Shibata, and K. Taniguchi, Phys. Rev.D , 124018 (2009).[34] L. website, URL .[35] T. Yamamoto, M. Shibata, and K. Taniguchi, Phys. Rev.D , 064054 (2008).[36] C. Reisswig and D. Pollney, ArXiv e-prints (2010),1006.1632.[37] D. Lai, F. A. Rasio, and S. L. Shapiro, Astrophys. J. Lett. , L63 (1993).[38] D. Lai, F. A. Rasio, and S. L. Shapiro, Astrophys. J. ,811 (1994), arXiv:astro-ph/9304027.[39] M. Boyle, D. A. Brown, L. E. Kidder, A. H. Mrou´e, H. P. Pfeiffer, M. A. Scheel, G. B. Cook, and S. A. Teukolsky,Phys. Rev. D , 124038 (2007), 0710.0158.[40] T. Damour and A. Nagar, ArXiv e-prints (2009),0906.1769.[41] T. Damour and A. Nagar, Phys. Rev. D , 081503(2009), 0902.0136.[42] L. Santamar´ıa, F. Ohme, P. Ajith, B. Br¨ugmann, N. Dor-band, M. Hannam, S. Husa, P. M¨osta, D. Pollney,C. Reisswig, et al., Phys. Rev. D , 064016 (2010),1005.3306.[43] P. Ajith, S. Babak, Y. Chen, M. Hewitson, B. Krishnan,A. M. Sintes, J. T. Whelan, B. Br¨ugmann, P. Diener,N. Dorband, et al., Phys. Rev. D , 104017 (2008),0710.2335.[44] P. Ajith, M. Hannam, S. Husa, Y. Chen, B. Br¨ugmann,N. Dorband, D. M¨uller, F. Ohme, D. Pollney, C. Reiss-wig, et al., Physical Review Letters , 241101 (2011),0909.2867.[45] R. Sturani, S. Fischetti, L. Cadonati, G. M. Guidi,J. Healy, D. Shoemaker, and A. Vicer´e, Journal ofPhysics Conference Series , 012007 (2010), 1005.0551.[46] R. Sturani, S. Fischetti, L. Cadonati, G. M. Guidi,J. Healy, D. Shoemaker, and A. Vicere’, ArXiv e-prints(2010), 1012.5172.[47] H. P. Pfeiffer, D. A. Brown, L. E. Kidder, L. Lindblom,G. Lovelace, and M. A. Scheel, Classical and QuantumGravity , 59 (2007), arXiv:gr-qc/0702106.[48] L. A. Wainstein and V. D. Zubakov,Extraction of signals from noise (Prentice-Hall, En-glewood Cliffs, NJ, 1962).[49] D. Shoemaker (LSC, 2009), URL https://dcc.ligo.org/cgi-bin/DocDB/ShowDocument?docid=2974 .[50] S. Hild, M. Abernathy, F. Acernese, P. Amaro-Seoane,N. Andersson, K. Arun, F. Barone, B. Barr, M. Barsug-lia, M. Beker, et al., Classical and Quantum Gravity ,094013 (2011), 1012.0908.[51] S. Hild, S. Chelkowski, and A. Freise, ArXiv e-prints(2008), 0810.0604.[52] S. Hild, S. Chelkowski, A. Freise, J. Franc, N. Morgado,R. Flaminio, and R. DeSalvo, Classical and QuantumGravity , 015003 (2010), 0906.2655.[53] B. Meers, Phys.Rev. D38 , 2317 (1988).[54] A. Buonanno and Y.-b. Chen, Phys.Rev.
D64 , 042006(2001), gr-qc/0102012.[55] J. D. E. Creighton and W. G. Anderson,Gravitational Wave Physics and Astronomy (Wiley-VCH, Berlin, 2011).[56] L. S. Finn, Gravitational wave interferometer noise calculator(2011), URL http://gwastro.org/for%20scientists/gravitational-wave-interferometer-noise-calculator .[57] S. A. Hughes, Phys. Rev. D , 102001 (2002), arXiv:gr-qc/0209012.[58] L. Lindblom, Phys. Rev. D , 103011 (2010),1009.0738.[59] F. ¨Ozel and D. Psaltis, Phys. Rev. D , 103003 (2009).[60] A. W. Steiner, J. M. Lattimer, and E. F. Brown, Astro-phys. J. , 33 (2010), 1005.0811.[61] F. ¨Ozel, G. Baym, and T. G¨uver, Phys. Rev. D ,101301 (2010), 1002.3153.[62] C. Markakis, J. S. Read, M. Shibata, K. Uryu, J. D. E.Creighton, and J. L. Friedman, ArXiv e-prints (2010),1008.1822.[63] Y. Pan, A. Buonanno, M. Boyle, L. T. Buchman, L. E. Kidder, H. P. Pfeiffer, and M. A. Scheel, ArXiv e-prints(2011), 1106.1021.[64] Y. Pan, A. Buonanno, L. T. Buchman, T. Chu, L. E.Kidder, H. P. Pfeiffer, and M. A. Scheel, Phys. Rev. D , 084041 (2010), 0912.3466.[65] Y. Pan, A. Buonanno, R. Fujita, E. Racine, andH. Tagoshi, Phys. Rev. D , 064003 (2011), 1006.0431.[66] A. Buonanno and T. Damour, Phys. Rev. D , 084006(1999), arXiv:gr-qc/9811091.[67] T. Damour, B. R. Iyer, and A. Nagar, Phys. Rev. D ,064004 (2009), 0811.2069.[68] E. Berti, V. Cardoso, and C. M. Will, Phys. Rev. D ,064030 (2006), arXiv:gr-qc/0512160.[69] T. Damour and A. Nagar, Phys. Rev. D , 064028(2007), 0705.2519.[70] T. Damour and A. Gopakumar, Phys. Rev. D , 124006(2006), arXiv:gr-qc/0602117.[71] L. E. Kidder, Phys. Rev. D , 044016 (2008),0710.0614.[72] A. Nagar (2010), private communication.[73] T. Damour, A. Nagar, E. N. Dorband, D. Pollney, andL. Rezzolla, Phys. Rev. D , 084017 (2008), 0712.3003.[74] A. Buonanno and T. Damour, Phys. Rev. D , 064015(2000), arXiv:gr-qc/0001013. Appendix A: Numerically evaluating the Fishermatrix
When an analytical representation of a waveform isnot available, the partial derivatives in the Fisher matrixEq. (22) must be evaluated numerically. There are sev-eral possible methods one can use, and we will examinetheir accuracy below.
1. Finite differencing of h ( t ; θ ) The simplest method, and that used in Ref. [12],is straightforward finite differencing of the signal h = F + h + + F × h × . For example, for five waveforms with val-ues of an EOS parameter θ given by { θ − , θ − , θ , θ , θ } with equal spacing ∆ θ , the three and five point centraldifferences are given by dhdθ = ∆ h ∆ θ + O ((∆ θ ) ) , where∆ h ∆ θ := − h ( t ; θ − ) + h ( t ; θ )∆ θ (A1) dhdθ = ∆ h ∆ θ + O ((∆ θ ) ) , where∆ hdθ := h ( t ; θ − ) − h ( t ; θ − ) + h ( t ; θ ) − h ( t ; θ )∆ θ . (A2)This finite differencing method is useful when waveformsdiffer only slightly: at each time t , on the scale ∆ θ thefunction h ( t ; θ ) is well approximated by the low orderinterpolating polynomials used to generate the finite dif-ferencing formulas.7This assumption fails when the waveforms used in thefinite differencing are significantly out of phase with eachother . The tidal interaction leads to a monotonicallyaccumulating phase difference relative to a BBH wave-form, implying that at a fixed time t the function h ( θ ; t )is an oscillating function of θ . Now if an oscillating func-tion h [Φ( θ )] = cos[Φ( θ )] has wavenumber k = Φ (cid:48) ( θ ) thatvaries slowly compared to Φ, then h (cid:48) ( θ ) is better approxi-mated by − sin(Φ)∆Φ / ∆ θ than by ∆ cos[Φ( θ )] / ∆ θ . Theassumption that k is slowly varying is k (cid:48) (cid:28) k , k (cid:48)(cid:48) (cid:28) k ,and the error in, for example, each of the two second-order discretizations is given to order ∆ θ by dhdθ − ∆ h ∆ θ = h ( θ )[ 16 k + O ( kk (cid:48) , k (cid:48)(cid:48) )]∆ θ ,dhdθ + sin[Φ( θ )] ∆ Φ∆ θ = h ( θ ) 16 k (cid:48)(cid:48) ∆ θ , (A3)with the error in the second expression much smaller thanthat in the first. We consider two ways to take advantageof this difference in accuracy.
2. Finite differencing of amplitude and phase
The first is to decompose each complex waveform intoan amplitude A and accumulated phase Φ h + ( t ; θ ) − ih × ( t ; θ ) = A ( t ; θ ) e − i Φ( t ; θ ) , (A4)where the accumulated phase of each waveform is a con-tinuous function defined by Φ = − arg( h + − ih × ) ± nπ for some integer n , and at the starting time t i the ac-cumulated phase of each waveform is chosen to be onthe branch n = 0. The advantage of this method isthat, at a fixed time, the functions A ( t ; θ ) and Φ( t ; θ ) arenon-oscillatory functions of θ even when the accumulatedphase difference between two waveforms is significantlymore than a radian.With this decomposition the gravitiational wave signalis h ( t ; θ ) = A ( t ; θ )( F + cos Φ( t ; θ ) + F × sin Φ( t ; θ )) , (A5)and the derivative of h is approximated by dhdθ = ∆ A ∆ θ ( F + cos Φ + F × sin Φ)+ A ( − F + sin Φ + F × cos Φ) ∆Φ∆ θ . (A6)If an intermediate waveform is not available to providethe functions A ( t ; θ ) and Φ( t ; θ ), they can be evaluatedby e.g. A ( t ; θ ) = ( A ( t ; θ − ) + A ( t ; θ )) / The dephasing of numerical waveforms is even more significantfor BNS inspiral. We believe that Ref. [12] which used thismethod underestimated the derivatives in some cases by a factorof ∼ We have found that this method works reasonably wellfor the inspiral waveform. If, however, the amplitude ofone of the numerical waveforms drops to zero, then thephase of the waveform becomes undefined. Because theamplitude of the numerical BHNS waveforms fall to zeroat different times for different EOS, as shown in Fig. 3,the derivative d Φ /dθ becomes meaningless towards theend when the average amplitude is still nonzero. It islikely one could work around this difficulty. However, wechoose instead to use another more robust method.
3. Finite differencing of Fourier transform
Because we will need to calculate the Fourier trans-form of the derivative dh/dθ to find the Fisher matrix,we first Fourier transform each waveform and then eval-uate the numerical derivative. Since the derivative d/dθ commutes with the Fourier transform, the Fisher matrixcan be written explicitly as (cid:18) ∂h∂θ A (cid:12)(cid:12)(cid:12)(cid:12) ∂h∂θ B (cid:19) = 4Re (cid:90) f f f i ∂ ˜ h∂θ A ∂ ˜ h ∗ ∂θ B S n ( f ) df, (A7)where the contribution to the integral below f i and above f f is negligable.As in the second method we break up each Fouriertransformed waveform into amplitude A ( f ; θ ) and accu-mulated phase Φ( f ; θ )˜ h ( f ; θ ) = A ( f ; θ ) e − i Φ( f ; θ ) , (A8)where the phase of each waveform at f i is on the n = 0branch cut. As demonstrated by Figs. 5 and 6, both theamplitude and phase are non-oscillatory functions of θ at a fixed frequency f , and can be well approximated bya low-order polynomial. In contrast to the accumulatedphase of the complex numerical waveform h + − ih × , theaccumulated phase of the Fourier transform of the strain˜ h is always well defined for numerical BHNS waveformsin the frequency range f i to f f .Finally, we find that one obtains better accuracy bydifferentiating ln A instead of A , decomposing ˜ h as˜ h ( f ; θ ) = e ln A ( f ; θ ) − i Φ( f ; θ ) . (A9)The derivative is now approximated by ∂ ˜ h∂θ = e ln A − i Φ (cid:18) ∆ ln A ∆ θ − i ∆Φ∆ θ (cid:19) . (A10)interpolating when needed to evaluate ln A and Φ at themidpoint.
4. Parameter spacing
Finally, we note that the EOS parameter spacing mustbe carefully chosen. If two waveforms are too close in8parameter space, the error in each waveform will domi-nate over the truncation error due to finite differencing.The most significant source of this error comes from thespurius oscillations in the amplitude of the Fourier trans-form in the frequency range ∼ ∼ Q = 2, we find that a spacing between waveforms of∆ log( p / (dyne cm − )) = 0 . . Q = 3, we havefound that a spacing of ∆ log( p / (dyne cm − )) ≥ . × pattern now form a+ pattern, and in the transformed parameter space thenew axes are not alligned with a degenerate contour. TheFisher matrix can be calculated in the transformed pa-rameter space then transformed back to the original pa-rameter space.We find that as long as these two requirements are met,the uncertainties in σ λ / and σ R have only an O (20%)fractional dependence on the EOS parameter spacing.However, the dependence of the orientation of the errorellipses on the EOS parameter spacing does not allow oneto distinguish between Λ and R as the best extracted pa-rameter. Appendix B: Prescription for calculating EOBwaveforms
In this appendix we compile the necessary ingredientsneeded to produce nonspinning BBH waveforms usingthe effective one body (EOB) formalism first introducedin Ref. [66]. The version used here is exactly that ofRef. [41], and is described in more detail in a review [40].The only ingredients not listed here are terms for the re-sumed waveform in Ref. [67] and coefficients to determinethe ringdown waveform found in Ref. [68].
1. Hamiltonian dynamics
In the EOB formalism the two-body dynamics are re-placed by a test particle of reduced mass µ = m m /M moving in a modified Schwarzschild metric of total mass M = m + m given by ds = − A ( r ) dt + B ( r ) dr + r ( dθ + sin θdφ ) . (B1)The metric potentials A and B can be calculated frompost-Newtonian theory. The first function is A ( u ) = P [1 − u +2 νu + (cid:18) − π (cid:19) νu + a νu + a νu ] , (B2)where u = 1 /r , ν = µ/M is the symmetric mass ratio,and P mn [ · ] denotes a Pad´e approximant of order m inthe numerator and n in the denominator. The 4 and5 PN coefficients, a and a , are fit to numerical BBHwaveforms. The values that give the optimal fit form adegenerate curve in the a – a parameter space, and thespecific values chosen here are ( a , a ) = (0 , − D ( r ) = B ( r ) A ( r ) , (B3)and has been calculated to 2PN order D ( u ) = P [1 − νu + 2(3 ν − νu ] . (B4)The motion of the EOB particle of mass µ is deter-mined by the Hamiltonianˆ H = 1 ν (cid:113) ν ( ˆ H eff − , (B5)whereˆ H eff = (cid:118)(cid:117)(cid:117)(cid:116) A (1 /r ) (cid:32) p φ r + p r B + 2 ν (4 − ν ) p r r (cid:33) (B6)is the effective Hamiltonian. The equations of motiongiven this conservative Hamiltonian ˆ H and a dissipative The expressions below will be written exclusively in terms ofrescaled dimensionless quantities. The coordinates (
T, R, φ ) andconjugate momenta ( P R , P φ ) have been rescaled to dimensionlesscoordinates ( t, r, φ ) and momenta ( p r , p φ ) given by: t = T/M and r = R/M for the coordinates, and p r = P R /µ , p φ = P φ /µM for the conjugate momenta. Other quantities are then rescaled inthe following way: ω = M Ω =
Mdφ/dT is the angular velocity,ˆ D = D/M is the distance to the source, ˆ H = H/µ and ˆ H eff = H eff /µ are the Hamiltonian and effective Hamiltonian, and ˆ F φ = F φ /µ is the radiation reaction force. F i are drdt = ∂ ˆ H∂p r (B7) dφdt = ∂ ˆ H∂p φ = ω (B8) dp r dt = − ∂ ˆ H∂r + ˆ F r (B9) dp φ dt = − ∂ ˆ H∂φ + ˆ F φ . (B10)Here, ∂ ˆ H∂φ = 0 because the EOB Hamiltonian does nothave an explicit φ dependence. In addition, for cir-cularized binary inspiral the radial component of theradiation-reaction force ˆ F r is of higher post-Newtonianorder than the tangential component, so it is set to zero.To increase resolution near the black hole, the radialcoordinate can be rewritten in terms of a tortoise coor-dinate [69] defined by dr ∗ dr = (cid:18) BA (cid:19) / . (B11)The new radial momentum is then p r ∗ = ( A/B ) / p r .Using this definition, the effective Hamiltonian becomesˆ H eff = (cid:118)(cid:117)(cid:117)(cid:116) p r ∗ + A (1 /r ) (cid:32) p φ r + 2 ν (4 − ν ) p r ∗ r (cid:33) (B12)where the parts that are 4PN and higher are neglected.(The 4 and 5 PN terms are however accounted for in thefree parameters a and a which were fit to numericalwaveforms). The equations of motion become drdt = A √ D ∂ ˆ H∂p r ∗ (B13) dφdt = ∂ ˆ H∂p φ = ω (B14) dp r ∗ dt = − A √ D ∂ ˆ H∂r (B15) dp φ dt = ˆ F φ . (B16)
2. Radiation reaction
For the radiation reaction term ˆ F φ , which is writtenin terms of the PN parameter x , we will need a way towrite x in terms of the dynamical variables. The usualmethod is to use the Newtonion potential 1 /r and veloc-ity squared ( ωr ) as PN counting parameters and thenrewrite them in terms of the gauge invariant angularvelocity ω using the Kepler law ω r = 1 which holdsin the Newtonian limit, and for circular orbits, in the Schwarzschild ( ν →
0) limit. The Kepler relation can beextended to circular orbits in the EOB metric by defininga new radial parameter, r ω = rψ / , where ψ ( r, p φ ) = 2 r (cid:18) dAdr (cid:19) − ν (cid:118)(cid:117)(cid:117)(cid:116) A ( r ) (cid:32) p φ r (cid:33) − , (B17)for which ω r ω = 1 holds for all circular orbits. In addi-tion, for noncircular orbits (in particular for the plunge),this relation also relaxes the quasicircular condition bynot requiring that the Kepler relation hold. The specificchoice of PN parameter used here is x = ( ωr ω ) . (B18)See Ref. [70] for an extensive discussion.The radiation reaction term ˆ F φ used in Ref. [41] takesthe form of a summation over all multipolesˆ F φ = − πνω (cid:88) (cid:96) =2 (cid:96) (cid:88) m =1 ( mω ) | ˆ Dh (cid:96)m | . (B19)Instead of the standard Taylor expanded version of h (cid:96)m which can be found in Ref. [71], Ref. [67] decomposed thewaveform into a product of terms: h = h Newt22 ˆ S eff T e iδ f ( x ) f NQC22 (B20)for (cid:96) = m = 2, and h (cid:96)m = h Newt (cid:96)m ˆ S eff T (cid:96)m e iδ (cid:96)m ρ (cid:96)(cid:96)m ( x ) (B21)for the other values of (cid:96) and m . The leading Newtonianpart h Newt (cid:96)m is given in the usual form as a function of xh Newt (cid:96)m = ν ˆ D n (cid:96)m c (cid:96) + (cid:15) ( ν ) x ( (cid:96) + (cid:15) ) / Y (cid:96) − (cid:15), − m (cid:16) π , φ (cid:17) (B22)where the coefficients n (cid:96)m and c (cid:96) + (cid:15) ( ν ) are defined byEqs. (5–7) of Ref. [67], and the parity (cid:15) is 0 for (cid:96) + m even and 1 for (cid:96) + m odd.The PN terms in the resummation which had beenwritten as functions of x in Ref. [67] are now written interms of the dynamical variables. The effective sourceterm ˆ S eff becomes [72]ˆ S eff ( r, p r ∗ , p φ ) = (cid:40) ˆ H eff ( r, p r ∗ , p φ ) (cid:15) = 0 p φ r ω ω (cid:15) = 1 . (B23)The tail term is T (cid:96)m ( r, p r ∗ , p φ ) = Γ( (cid:96) + 1 − i ˆˆ k )Γ( (cid:96) + 1) e π ˆˆ k e i ˆˆ k ln 2 kr , (B24)where ˆˆ k = νm ˆ H ( r, p r ∗ , p φ ) ω ( r, p r ∗ , p φ ), k = mω ( r, p r ∗ , p φ ), and r = 2. The phase of this tailterm is corrected with a term of the form e iδ (cid:96)m . The0first ten δ (cid:96)m are given in Eqs. (20–29) of Ref. [67]. Thefirst one is δ = 73 y / + 428 π y − ν ¯ y / , (B25)where y = ( ν ˆ H ( r, p r ∗ , p φ ) ω ( r, p r ∗ , p φ )) / and ¯ y , whichhas several possible forms, is chosen to be ¯ y = ω / [72]. Finally, the remainder term of the resummation f (cid:96)m isexpanded in powers of x . For (cid:96) = m = 2 this is thenre-summed with a Pad´e approximant f ( x ) = P [ f Taylor22 ( x )] , (B26)where f Taylor22 ( ν, x ) = 1 + 55 ν − x + 2047 ν − ν − x + (cid:18) ν − ν π ν − ν − (x) + 21428357727650 (cid:19) x + (cid:18) (x) − (cid:19) x + (cid:18) (x) − (cid:19) x , (B27)and the eulerln m (x) = γ E +ln 2+ ln x+ln m terms are treated as coefficients when calculating the Pad´e approximant.For the other values of (cid:96) and m , f (cid:96)m is re-summed in the form f (cid:96)m = ρ (cid:96)(cid:96)m . The quantity ρ (cid:96)m is given in Eqs. (C1–C35)of Ref. [67]. ρ is for example ρ = 1 + (cid:18) ν − (cid:19) x + (cid:18) ν − ν − (cid:19) x + (cid:18) − (x) (cid:19) x + (cid:18) (x) − (cid:19) x . (B28)The final product in the resummation of h is a nextto quasicircular (NQC) correction term that is used tocorrect the dynamics and waveform amplitude during theplunge f NQC22 ( a , a ) = 1 + a p r ∗ ( rω ) + a ¨ rrω . (B29)The free parameters a and a are determined by thefollowing conditions: (i) the time when the orbital fre-quency ω is a maximum (the EOB merger time t M ) co-incides with the time when the amplitude | h | is a max-imum, and (ii) the value of the maximum amplitude isequal to a fitting function that was fit to several BBHsimulations and is given by | h | max ( ν ) = 1 . ν (1 − . − ν )) . (B30)
3. Integrating the equations of motion
The equations of motion are solved by starting withinitial conditions { r , φ , p r ∗ , p φ } and numerically inte-grating the equations of motion. In this paper we areinterested in long, zero-eccentricity orbits. This can beachieved in the EOB framework by starting the integra-tion with large r , where radiation reaction effects aresmall, and using the quasicircular condition p r ∗ = 0. Eq. (B15) then becomes ∂H∂r ( r, p r ∗ = 0 , p φ ) = 0 (B31)and results in the condition p φ = − ddu A ( u ) ddu ( u A ( u )) (B32)for p φ . If this quasicircular initial condition is used forsmaller r , the radiation reaction term is no longer neg-ligable, and this initial condition will result in eccentricorbits. If desired, one can use an initial condition thatmore accurately approximates a zero eccentricity inspiralsuch as post-circular or post-post-circular initial condi-tions with nonzero p r ∗ [73].To numerically solve Eqs. (B13–B16), they must bewritten as a system of first order equations. However,the term ˆ F φ in Eq. (B16) contains the square of ¨ r fromthe NQC term f NQC22 in h . Since f NQC22 gives a smallcorrection of order 10% during the plunge, the easiestmethod, and that used in Ref. [41], is iteration [72]: (i)First solve the system of equations with f NQC22 set to one.(ii) Use the solution of Eqs. (B13–B16) to evaluate ¨ r andthe other quantities in f NQC22 . (iii) Re-solve the equationsof motion with the NQC coefficients no longer set to one.(iv) Repeat steps (ii) and (iii) until the solution converges1to the desired accuracy. In practice this iteration onlyneeds to be done roughly 2–5 times.A second method is to directly rewrite Eq. (B16) as afirst order equation. This can be done by replacing ¨ r inthe NQC term on the right hand side with an expressioncontaining ˙ p φ and then solving for ˙ p φ . The equations ofmotion (B13–B16) and the chain rule give¨ r = ddt (cid:32) A √ D ∂ ˆ H∂p r ∗ (cid:33) (B33)= L + M + N ˙ p φ (B34)where L = 12 ∂∂r A D (cid:32) ∂ ˆ H∂p r ∗ (cid:33) (B35) M = − A D ∂ ˆ H∂r ∂ ˆ H∂p r ∗ (B36) N = A √ D ∂ ˆ H∂p r ∗ ∂p φ . (B37)Plugging this expression into Eq. (B16) yields an equa-tion quadratic in ˙ p φ which can be solved exactly ifdesired. To first order in the NQC correction term,Eq. (B16) now becomes the first order equation dp φ dt = ˆ F φ, Higher + ˆ F QC φ, (cid:104) a p r ∗ ( rω ) + 2 a rω ( L + M ) (cid:105) − F QC φ, a rω N , (B38)whereˆ F φ, Higher = − πνω (cid:88) (cid:96) =2 (cid:96) (cid:88) m =1( (cid:96),m ) (cid:54) =(2 , ( mω ) | ˆ Dh (cid:96)m | (B39)includes just the higher order terms ( (cid:96), m ) (cid:54) = (2 , F QC φ, = − πνω (2 ω ) | ˆ Dh QC22 | . (B40)The QC in h QC22 means that the NQC term f NQC22 hasbeen factored out.The solution to the equations of motion { r ( t ) , φ ( t ) , p r ∗ ( t ) , p φ ( t ) } are then plugged back intoEqs. (B20–B21) to give the waveform h inspiral (cid:96)m ( t ).
4. Ringdown
In the EOB formalism the ringdown waveform of thefinal Kerr black hole is smoothly matched onto the in- spiral waveform at the EOB merger time t M . The massof the black hole remnant is given by the energy of theEOB particle at the merger time t M M BH ≡ µ ˆ H ( t M ) = M (cid:113) ν ( ˆ H eff ( t M ) − , (B41)and the Kerr parameter is given by the final angular mo-mentum of the EOB particle [74]ˆ a BH ≡ P φ ( t M ) M = νp φ ( t M )1 + 2 ν ( ˆ H eff ( t M ) − . (B42)The ringdown waveform is given by the first five posi-tive quasinormal modes (QNM) for a black hole of mass M BH and spin ˆ a BH : h ringdown22 ( t ) = 1ˆ D (cid:88) n =0 C +22 n e − σ +22 n ( t − t M ) , (B43)where σ +22 n = α n + iω n is the nth complex (cid:96) = m = 2QNM frequency for a Kerr BH with mass ˆ M BH and spinˆ a BH , and C +22 n are complex constants that determine themagnitude and phase of each QNM. The amplitude ofthe negative frequency modes is small [69]. The firstthree QNMs have been tabulated in Ref. [68], and fittingformuli are also provided. The QNM frequency ω n canbe approximated by M BH ω n = f + f (1 − ˆ a BH ) f , (B44)and the inverse damping time α n is given in terms ofthe quality factor approximated by12 ω n α n = q + q (1 − ˆ a BH ) q . (B45)The coefficients for n = 0–2 can be found in table VIIIof Ref. [68]. For n = 3–4, α n and ω n can be linearlyextrapolated from the values for n = 1 and 2 as was donein Ref. [73].The constants C +22 n are determined by requiring thatthe inspiral and ringdown waveforms be continuous on a“matching comb” centered on the EOB merger time t M .Specifically, at the times { t M − δ, t M − δ, t M , t M + δ, t M +2 δ } we require h inspiral22 ( t ) = h ringdown22 ( t ). In Ref. [41], δ was chosen to be equal to 2 . M BH /M . This gives 5 com-plex equations for the 5 unknown complex coefficients C +22 n .The full inspiral plus ringdown waveform is then givenby h ( t ) = (cid:40) h inspiral22 ( t ) t < t M h ringdown22 ( t ) t > t M ..